Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
LinearSystem.hpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF 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
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
73private:
74
87 MatNullSpace mMatNullSpace;
91
96 double mTolerance;
102 std::string mKspType;
103 std::string mPcType;
108
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
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
176public:
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
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
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
341 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
471
475 Vec GetRhsVector() const;
476
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
584namespace boost
585{
586namespace serialization
587{
588template<class Archive>
589inline 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
614template<class Archive>
615inline 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_
gcov doesn't like this file...
#define CHASTE_CLASS_EXPORT(T)
static std::string GetArchiveDirectory()
Vec & rGetDirichletBoundaryConditionsVector()
void SetNullBasis(Vec nullbasis[], unsigned numberOfBases)
void serialize(Archive &archive, const unsigned int version)
void SetPrecondMatrixIsDifferentFromLhs(bool precondIsDifferent=true)
boost::shared_ptr< std::vector< PetscInt > > mpBathNodes
Vec Solve(Vec lhsGuess=nullptr)
void FinaliseLhsMatrix()
Vec mDirichletBoundaryConditionsVector
PetscReal mEigMin
void SetMatrixIsSymmetric(bool isSymmetric=true)
void AssembleFinalLinearSystem()
Mat GetLhsMatrix() const
bool mUseAbsoluteTolerance
void ZeroMatrixColumn(PetscInt col)
bool mForceSpectrumReevaluation
void SetMatrixIsConstant(bool matrixIsConstant)
unsigned mRowPreallocation
void FinalisePrecondMatrix()
void GetOwnershipRange(PetscInt &lo, PetscInt &hi)
Mat & rGetLhsMatrix()
PCLDUFactorisation * mpLDUFactorisationPC
void SwitchWriteModeLhsMatrix()
void AddToMatrixElement(PetscInt row, PetscInt col, double value)
double GetMatrixElement(PetscInt row, PetscInt col)
PetscInt mOwnershipRangeLo
void SetKspType(const char *kspType)
void SetMatrixElement(PetscInt row, PetscInt col, double value)
unsigned GetNumIterations() const
std::string mKspType
bool IsMatrixSymmetric()
PetscInt mOwnershipRangeHi
Vec GetRhsVector() const
void SetUseFixedNumberIterations(bool useFixedNumberIterations=true, unsigned evaluateNumItsEveryNSolves=UINT_MAX)
void FinaliseRhsVector()
std::string mPcType
MatNullSpace mMatNullSpace
void ZeroMatrixRowsWithValueOnDiagonal(std::vector< unsigned > &rRows, double diagonalValue)
unsigned GetSize() const
void ZeroMatrixRowsAndColumnsWithValueOnDiagonal(std::vector< unsigned > &rRowColIndices, double diagonalValue)
void SetAbsoluteTolerance(double absoluteTolerance)
void AssembleIntermediateLinearSystem()
void SetRhsVectorElement(PetscInt row, double value)
void SetMatrixRow(PetscInt row, double value)
void * mpConvergenceTestContext
PCTwoLevelsBlockDiagonal * mpTwoLevelsBlockDiagonalPC
Mat & rGetPrecondMatrix()
void ZeroLinearSystem()
Vec GetMatrixRowDistributed(unsigned rowIndex)
PetscReal mEigMax
friend class boost::serialization::access
PCBlockDiagonal * mpBlockDiagonalPC
void AddLhsMultipleValues(unsigned *matrixRowAndColIndices, c_matrix< double, MATRIX_SIZE, MATRIX_SIZE > &rSmallMatrix)
void AddToRhsVectorElement(PetscInt row, double value)
unsigned mEvaluateNumItsEveryNSolves
Vec & rGetRhsVector()
bool mPrecondMatrixIsNotLhs
void SetRelativeTolerance(double relativeTolerance)
void RemoveNullSpace()
bool mUseFixedNumberIterations
double GetRhsVectorElement(PetscInt row)
void SetPcType(const char *pcType, boost::shared_ptr< std::vector< PetscInt > > pBathNodes=boost::shared_ptr< std::vector< PetscInt > >())
void AddRhsMultipleValues(unsigned *vectorIndices, c_vector< double, VECTOR_SIZE > &smallVector)
unsigned mNumSolves
double mNonZerosUsed
static void SetOption(Mat matrix, MatOption option)
static void AddMultipleValues(Mat matrix, unsigned *matrixRowAndColIndices, c_matrix< double, MATRIX_SIZE, MATRIX_SIZE > &rSmallMatrix)
static void DumpPetscObject(const Mat &rMat, const std::string &rOutputFileFullPath)
static void ReadPetscObject(Mat &rMat, const std::string &rOutputFileFullPath, Vec rParallelLayout=nullptr)
static void AddMultipleValues(Vec vector, unsigned *vectorIndices, c_vector< double, VECTOR_SIZE > &smallVector)