Chaste  Release::2018.1
LinearSystem.hpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #ifndef _LINEARSYSTEM_HPP_
37 #define _LINEARSYSTEM_HPP_
38 
39 #include "ChasteSerialization.hpp"
40 #include "UblasCustomFunctions.hpp" // needs to be 'first'
41 
42 #include "PetscTools.hpp"
43 #include "PetscVecTools.hpp"
44 #include "PetscMatTools.hpp"
45 #include "OutputFileHandler.hpp"
46 #include "PCBlockDiagonal.hpp"
47 #include "PCLDUFactorisation.hpp"
48 #include "PCTwoLevelsBlockDiagonal.hpp"
49 #include "ArchiveLocationInfo.hpp"
50 //#include <boost/serialization/shared_ptr.hpp>
51 
52 #include <petscvec.h>
53 #include <petscmat.h>
54 #include <petscksp.h>
55 #include <petscviewer.h>
56 
57 #include <string>
58 #include <cassert>
59 
66 {
67  friend class TestLinearSystem;
68  friend class TestPCBlockDiagonal;
69  friend class TestPCTwoLevelsBlockDiagonal;
70  friend class TestPCLDUFactorisation;
71  friend class TestChebyshevIteration;
72 
73 private:
74 
87  MatNullSpace mMatNullSpace;
91 
92  KSP mKspSolver;
93  bool mKspIsSetup;
94  double mNonZerosUsed;
96  double mTolerance;
102  std::string mKspType;
103  std::string mPcType;
107 
114 
116  boost::shared_ptr<std::vector<PetscInt> > mpBathNodes;
117 
120 
123 
126 
133 
136 
138  unsigned mNumSolves;
139 
141  PetscReal mEigMin;
142 
144  PetscReal mEigMax;
145 
148 
149 #ifdef TRACE_KSP
150  unsigned mTotalNumIterations;
151  unsigned mMaxNumIterations;
152 #endif
153 
161  template<class Archive>
162  void serialize(Archive & archive, const unsigned int version)
163  {
164  //MatNullSpace mMatNullSpace; // Gets re-created by calling code on load
165 
166  archive & mNonZerosUsed;
167  archive & mMatrixIsConstant;
168  archive & mTolerance;
169  archive & mUseAbsoluteTolerance;
170  archive & mKspType;
171  archive & mPcType;
172 
173  //Vec mDirichletBoundaryConditionsVector; // Gets re-created by calling code on load
174  }
175 
176 public:
177 
188  LinearSystem(PetscInt lhsVectorSize, unsigned rowPreallocation=UINT_MAX);
189 
207  LinearSystem(Vec templateVector, unsigned rowPreallocation, bool newAllocationError=true);
208 
221  LinearSystem(Vec residualVector, Mat jacobianMatrix);
222 
230  LinearSystem(PetscInt lhsVectorSize, Mat lhsMatrix, Vec rhsVector);
231 
235  ~LinearSystem();
236 
237 // bool IsMatrixEqualTo(Mat testMatrix);
238 // bool IsRhsVectorEqualTo(Vec testVector);
239 
247  void SetMatrixElement(PetscInt row, PetscInt col, double value);
248 
256  void AddToMatrixElement(PetscInt row, PetscInt col, double value);
257 
264 
271 
275  void FinaliseLhsMatrix();
276 
280  void FinalisePrecondMatrix();
281 
286 
290  void FinaliseRhsVector();
291 
297  void SetMatrixIsSymmetric(bool isSymmetric=true);
298 
304  bool IsMatrixSymmetric();
305 
311  void SetMatrixIsConstant(bool matrixIsConstant);
312 
318  void SetRelativeTolerance(double relativeTolerance);
319 
325  void SetAbsoluteTolerance(double absoluteTolerance);
326 
332  void SetKspType(const char* kspType);
333 
340  void SetPcType(const char* pcType, boost::shared_ptr<std::vector<PetscInt> > pBathNodes=boost::shared_ptr<std::vector<PetscInt> >() );
342 
346  void DisplayMatrix();
347 
351  void DisplayRhs();
352 
361  void SetMatrixRow(PetscInt row, double value);
362 
369  Vec GetMatrixRowDistributed(unsigned rowIndex);
370 
379  void ZeroMatrixRowsWithValueOnDiagonal(std::vector<unsigned>& rRows, double diagonalValue);
380 
381 
389  void ZeroMatrixRowsAndColumnsWithValueOnDiagonal(std::vector<unsigned>& rRowColIndices, double diagonalValue);
390 
391 
401  void ZeroMatrixColumn(PetscInt col);
402 
406  void ZeroLhsMatrix();
407 
411  void ZeroRhsVector();
412 
416  void ZeroLinearSystem();
417 
423  Vec Solve(Vec lhsGuess=nullptr);
424 
431  void SetRhsVectorElement(PetscInt row, double value);
432 
439  void AddToRhsVectorElement(PetscInt row, double value);
440 
444  unsigned GetSize() const;
445 
457  void SetNullBasis(Vec nullbasis[], unsigned numberOfBases);
458 
465  void RemoveNullSpace();
466 
470  Vec& rGetRhsVector();
471 
475  Vec GetRhsVector() const;
476 
480  Mat& rGetLhsMatrix();
481 
485  Mat GetLhsMatrix() const;
486 
491 
498 
499  // DEBUGGING CODE:
506  void GetOwnershipRange(PetscInt& lo, PetscInt& hi);
507 
515  double GetMatrixElement(PetscInt row, PetscInt col);
516 
523  double GetRhsVectorElement(PetscInt row);
524 
528  unsigned GetNumIterations() const;
529 
539  template<size_t MATRIX_SIZE>
540  void AddLhsMultipleValues(unsigned* matrixRowAndColIndices, c_matrix<double, MATRIX_SIZE, MATRIX_SIZE>& rSmallMatrix)
541  {
542  PetscMatTools::AddMultipleValues(mLhsMatrix, matrixRowAndColIndices, rSmallMatrix);
543  }
544 
554  template<size_t VECTOR_SIZE>
555  void AddRhsMultipleValues(unsigned* vectorIndices, c_vector<double, VECTOR_SIZE>& smallVector)
556  {
557  PetscVecTools::AddMultipleValues(mRhsVector, vectorIndices, smallVector);
558  }
559 
564  void SetPrecondMatrixIsDifferentFromLhs(bool precondIsDifferent = true);
565 
571  void SetUseFixedNumberIterations(bool useFixedNumberIterations = true, unsigned evaluateNumItsEveryNSolves = UINT_MAX);
572 
577  void ResetKspSolver();
578 };
579 
581 // Declare identifier for the serializer
583 
584 namespace boost
585 {
586 namespace serialization
587 {
588 template<class Archive>
589 inline void save_construct_data(
590  Archive & ar, const LinearSystem * t, const unsigned int file_version)
591 {
592 
593  std::string archive_filename_lhs = ArchiveLocationInfo::GetArchiveDirectory() + "lhs.mat";
594  std::string archive_filename_rhs = ArchiveLocationInfo::GetArchiveDirectory() + "rhs.vec";
595  const unsigned size = t->GetSize();
596  ar << size;
597 
598  PetscTools::DumpPetscObject(t->GetRhsVector(), archive_filename_rhs);
599  PetscTools::DumpPetscObject(t->GetLhsMatrix(), archive_filename_lhs);
600 
601  //Is the matrix structurally symmetric?
602  PetscBool symm_set, is_symmetric;
603  is_symmetric = PETSC_FALSE;
604  //Note that the following call only changes is_symmetric when symm_set is true
605  MatIsSymmetricKnown(t->GetLhsMatrix(), &symm_set, &is_symmetric);
606  assert(symm_set == is_symmetric);
607  ar << symm_set;
608 }
609 
614 template<class Archive>
615 inline void load_construct_data(
616  Archive & ar, LinearSystem * t, const unsigned int file_version)
617 {
618  std::string archive_filename_lhs = ArchiveLocationInfo::GetArchiveDirectory() + "lhs.mat";
619  std::string archive_filename_rhs = ArchiveLocationInfo::GetArchiveDirectory() + "rhs.vec";
620 
621  PetscInt size;
622  ar >> size;
623 
624  Vec new_vec;
625  PetscTools::ReadPetscObject(new_vec, archive_filename_rhs);
626 
627  Mat new_mat;
628  PetscTools::ReadPetscObject(new_mat, archive_filename_lhs);
629 
630  //This has to occur after the call to MatLoad as the matrix does not exist until MatLoad is called.
631  //The property will be copied & set correctly in the LinearSystem constructor.
632  PetscBool symm_set;
633 
634  ar >> symm_set;
635  if (symm_set == PETSC_TRUE)
636  {
637  PetscMatTools::SetOption(new_mat, MAT_SYMMETRIC);
638  PetscMatTools::SetOption(new_mat, MAT_SYMMETRY_ETERNAL);
639  }
640 
641  ::new(t)LinearSystem(size, new_mat, new_vec);
642 }
643 }
644 } // namespace ...
645 
646 #endif //_LINEARSYSTEM_HPP_
void AssembleFinalLinearSystem()
Mat & rGetPrecondMatrix()
static void AddMultipleValues(Mat matrix, unsigned *matrixRowAndColIndices, c_matrix< double, MATRIX_SIZE, MATRIX_SIZE > &rSmallMatrix)
void FinalisePrecondMatrix()
unsigned mRowPreallocation
void FinaliseLhsMatrix()
friend class boost::serialization::access
void AddToRhsVectorElement(PetscInt row, double value)
void GetOwnershipRange(PetscInt &lo, PetscInt &hi)
static void DumpPetscObject(const Mat &rMat, const std::string &rOutputFileFullPath)
Definition: PetscTools.cpp:333
bool mUseAbsoluteTolerance
void * mpConvergenceTestContext
void serialize(Archive &archive, const unsigned int version)
void ZeroLinearSystem()
void AddLhsMultipleValues(unsigned *matrixRowAndColIndices, c_matrix< double, MATRIX_SIZE, MATRIX_SIZE > &rSmallMatrix)
void SetKspType(const char *kspType)
Vec Solve(Vec lhsGuess=nullptr)
void ZeroMatrixRowsWithValueOnDiagonal(std::vector< unsigned > &rRows, double diagonalValue)
unsigned GetSize() const
void ZeroRhsVector()
unsigned mEvaluateNumItsEveryNSolves
void SetRelativeTolerance(double relativeTolerance)
void SetMatrixElement(PetscInt row, PetscInt col, double value)
void AddToMatrixElement(PetscInt row, PetscInt col, double value)
double mTolerance
Vec mDirichletBoundaryConditionsVector
Vec & rGetDirichletBoundaryConditionsVector()
PCTwoLevelsBlockDiagonal * mpTwoLevelsBlockDiagonalPC
bool mMatrixIsConstant
void ZeroMatrixRowsAndColumnsWithValueOnDiagonal(std::vector< unsigned > &rRowColIndices, double diagonalValue)
double GetRhsVectorElement(PetscInt row)
boost::shared_ptr< std::vector< PetscInt > > mpBathNodes
Mat & rGetLhsMatrix()
double GetMatrixElement(PetscInt row, PetscInt col)
unsigned mNumSolves
void RemoveNullSpace()
void ZeroMatrixColumn(PetscInt col)
void SetPrecondMatrixIsDifferentFromLhs(bool precondIsDifferent=true)
PetscInt mOwnershipRangeHi
std::string mKspType
Vec & rGetRhsVector()
void SetMatrixIsSymmetric(bool isSymmetric=true)
LinearSystem(PetscInt lhsVectorSize, unsigned rowPreallocation=UINT_MAX)
PCBlockDiagonal * mpBlockDiagonalPC
void SetNullBasis(Vec nullbasis[], unsigned numberOfBases)
MatNullSpace mMatNullSpace
PCLDUFactorisation * mpLDUFactorisationPC
void SetPcType(const char *pcType, boost::shared_ptr< std::vector< PetscInt > > pBathNodes=boost::shared_ptr< std::vector< PetscInt > >())
bool mDestroyMatAndVec
void SetAbsoluteTolerance(double absoluteTolerance)
void ZeroLhsMatrix()
Vec GetRhsVector() const
static void AddMultipleValues(Vec vector, unsigned *vectorIndices, c_vector< double, VECTOR_SIZE > &smallVector)
void FinaliseRhsVector()
bool mForceSpectrumReevaluation
void SetRhsVectorElement(PetscInt row, double value)
Mat GetLhsMatrix() const
std::string mPcType
PetscReal mEigMax
void AssembleIntermediateLinearSystem()
void ResetKspSolver()
double mNonZerosUsed
static void SetOption(Mat matrix, MatOption option)
bool IsMatrixSymmetric()
void AddRhsMultipleValues(unsigned *vectorIndices, c_vector< double, VECTOR_SIZE > &smallVector)
#define CHASTE_CLASS_EXPORT(T)
void SwitchWriteModeLhsMatrix()
PetscInt mSize
unsigned GetNumIterations() const
static std::string GetArchiveDirectory()
static void ReadPetscObject(Mat &rMat, const std::string &rOutputFileFullPath, Vec rParallelLayout=nullptr)
Definition: PetscTools.cpp:361
PetscInt mOwnershipRangeLo
gcov doesn't like this file...
void SetUseFixedNumberIterations(bool useFixedNumberIterations=true, unsigned evaluateNumItsEveryNSolves=UINT_MAX)
bool mPrecondMatrixIsNotLhs
void SetMatrixRow(PetscInt row, double value)
PetscReal mEigMin
Vec GetMatrixRowDistributed(unsigned rowIndex)
void SetMatrixIsConstant(bool matrixIsConstant)
bool mUseFixedNumberIterations
void DisplayMatrix()