LinearSystem.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #ifndef _LINEARSYSTEM_HPP_
00031 #define _LINEARSYSTEM_HPP_
00032 
00033 #include "ChasteSerialization.hpp"
00034 #include "UblasCustomFunctions.hpp" // needs to be 'first'
00035 
00036 #include "PetscTools.hpp"
00037 #include "PetscVecTools.hpp"
00038 #include "PetscMatTools.hpp"
00039 #include "OutputFileHandler.hpp"
00040 #include "PCBlockDiagonal.hpp"
00041 #include "PCLDUFactorisation.hpp"
00042 #include "PCTwoLevelsBlockDiagonal.hpp"
00043 #include "ArchiveLocationInfo.hpp"
00044 #include <boost/serialization/shared_ptr.hpp>
00045 
00046 #include <petscvec.h>
00047 #include <petscmat.h>
00048 #include <petscksp.h>
00049 #include <petscviewer.h>
00050 
00051 #include <string>
00052 #include <cassert>
00053 
00059 class LinearSystem
00060 {
00061     friend class TestLinearSystem;
00062     friend class TestPCBlockDiagonal;
00063     friend class TestPCTwoLevelsBlockDiagonal;
00064     friend class TestPCLDUFactorisation;
00065     friend class TestChebyshevIteration;
00066 
00067 private:
00068 
00069     Mat mLhsMatrix;  
00070     Mat mPrecondMatrix; 
00071     Vec mRhsVector;  
00072     PetscInt mSize;  
00078     PetscInt mOwnershipRangeLo; 
00079     PetscInt mOwnershipRangeHi; 
00081     MatNullSpace mMatNullSpace; 
00084     bool mDestroyMatAndVec;
00085 
00086     KSP mKspSolver;   
00087     bool mKspIsSetup; 
00088     double mNonZerosUsed;  
00089     bool mMatrixIsConstant; 
00090     double mTolerance; 
00095     bool mUseAbsoluteTolerance;
00096     std::string mKspType;
00097     std::string mPcType;
00099     Vec mDirichletBoundaryConditionsVector; 
00101 
00102 
00103     PCBlockDiagonal* mpBlockDiagonalPC;
00105     PCLDUFactorisation* mpLDUFactorisationPC;
00107     PCTwoLevelsBlockDiagonal* mpTwoLevelsBlockDiagonalPC;
00108     
00110     boost::shared_ptr<std::vector<PetscInt> > mpBathNodes;
00111 
00113     bool mPrecondMatrixIsNotLhs;
00114 
00116     unsigned mRowPreallocation;
00117     
00119     bool mUseFixedNumberIterations;
00120         
00126     unsigned mEvaluateNumItsEveryNSolves;
00127 
00129     void* mpConvergenceTestContext;
00130 
00132     unsigned mNumSolves;
00133 
00135     PetscReal mEigMin;
00136 
00138     PetscReal mEigMax;
00139     
00141     bool mForceSpectrumReevaluation;
00142 
00143 #ifdef TRACE_KSP
00144     unsigned mTotalNumIterations;
00145     unsigned mMaxNumIterations;
00146 #endif
00147 
00148     friend class boost::serialization::access;
00155     template<class Archive>
00156     void serialize(Archive & archive, const unsigned int version)
00157     {
00158         //MatNullSpace mMatNullSpace; // Gets re-created by calling code on load
00159 
00160         archive & mNonZerosUsed;
00161         archive & mMatrixIsConstant;
00162         archive & mTolerance;
00163         archive & mUseAbsoluteTolerance;
00164         archive & mKspType;
00165         archive & mPcType;
00166 
00167         //Vec mDirichletBoundaryConditionsVector; // Gets re-created by calling code on load
00168     }
00169 
00170 public:
00171 
00182     LinearSystem(PetscInt lhsVectorSize, unsigned rowPreallocation=UINT_MAX);
00183 
00197     LinearSystem(Vec templateVector, unsigned rowPreallocation);
00198 
00211     LinearSystem(Vec residualVector, Mat jacobianMatrix);
00212 
00220     LinearSystem(PetscInt lhsVectorSize, Mat lhsMatrix, Vec rhsVector);
00221 
00225     ~LinearSystem();
00226 
00227 
00228 //    bool IsMatrixEqualTo(Mat testMatrix);
00229 //    bool IsRhsVectorEqualTo(Vec testVector);
00230 
00238     void SetMatrixElement(PetscInt row, PetscInt col, double value);
00239 
00247     void AddToMatrixElement(PetscInt row, PetscInt col, double value);
00248 
00254     void AssembleFinalLinearSystem();
00255 
00261     void AssembleIntermediateLinearSystem();
00262 
00266     void AssembleFinalLhsMatrix();
00267 
00271     void AssembleFinalPrecondMatrix();
00272 
00276     void AssembleIntermediateLhsMatrix();
00277 
00281     void AssembleRhsVector();
00282 
00288     void SetMatrixIsSymmetric(bool isSymmetric=true);
00289 
00295     bool IsMatrixSymmetric();
00296 
00302     void SetMatrixIsConstant(bool matrixIsConstant);
00303 
00309     void SetRelativeTolerance(double relativeTolerance);
00310 
00316     void SetAbsoluteTolerance(double absoluteTolerance);
00317 
00323     void SetKspType(const char* kspType);
00324 
00331 
00332     void SetPcType(const char* pcType, boost::shared_ptr<std::vector<PetscInt> > pBathNodes=boost::shared_ptr<std::vector<PetscInt> >() );
00333 
00337     void DisplayMatrix();
00338 
00342     void DisplayRhs();
00343 
00352     void SetMatrixRow(PetscInt row, double value);
00353 
00360     Vec GetMatrixRowDistributed(unsigned rowIndex);
00361 
00370     void ZeroMatrixRowsWithValueOnDiagonal(std::vector<unsigned>& rRows, double diagonalValue);
00371 
00372 
00379     void ZeroMatrixRowsAndColumnsWithValueOnDiagonal(std::vector<unsigned>& rRowColIndices, double diagonalValue);
00380 
00381 
00391     void ZeroMatrixColumn(PetscInt col);
00392 
00396     void ZeroLhsMatrix();
00397 
00401     void ZeroRhsVector();
00402 
00406     void ZeroLinearSystem();
00407 
00413     Vec Solve(Vec lhsGuess=NULL);
00414 
00421     void SetRhsVectorElement(PetscInt row, double value);
00422 
00429     void AddToRhsVectorElement(PetscInt row, double value);
00430 
00434     unsigned GetSize() const;
00435 
00448     void SetNullBasis(Vec nullbasis[], unsigned numberOfBases);
00449 
00456     void RemoveNullSpace();
00457 
00461     Vec& rGetRhsVector();
00462 
00466     Vec GetRhsVector() const;
00467 
00471     Mat& rGetLhsMatrix();
00472 
00476     Mat GetLhsMatrix() const;
00477 
00481     Mat& rGetPrecondMatrix();
00482 
00488     Vec& rGetDirichletBoundaryConditionsVector();
00489 
00490     // DEBUGGING CODE:
00497     void GetOwnershipRange(PetscInt& lo, PetscInt& hi);
00498 
00506     double GetMatrixElement(PetscInt row, PetscInt col);
00507 
00514     double GetRhsVectorElement(PetscInt row);
00515 
00519     unsigned GetNumIterations() const;
00520 
00530     template<size_t MATRIX_SIZE>
00531     void AddLhsMultipleValues(unsigned* matrixRowAndColIndices, c_matrix<double, MATRIX_SIZE, MATRIX_SIZE>& rSmallMatrix)
00532     {
00533         PetscMatTools::AddMultipleValues(mLhsMatrix, matrixRowAndColIndices, rSmallMatrix);
00534     }
00535 
00545     template<size_t VECTOR_SIZE>
00546     void AddRhsMultipleValues(unsigned* vectorIndices, c_vector<double, VECTOR_SIZE>& smallVector)
00547     {
00548         PetscVecTools::AddMultipleValues(mRhsVector, vectorIndices, smallVector);
00549     }
00550 
00555     void SetPrecondMatrixIsDifferentFromLhs(bool precondIsDifferent = true);
00556     
00562     void SetUseFixedNumberIterations(bool useFixedNumberIterations = true, unsigned evaluateNumItsEveryNSolves = UINT_MAX);
00563     
00568     void ResetKspSolver();
00569 };
00570 
00571 #include "SerializationExportWrapper.hpp"
00572 // Declare identifier for the serializer
00573 CHASTE_CLASS_EXPORT(LinearSystem)
00574 
00575 namespace boost
00576 {
00577 namespace serialization
00578 {
00579 
00580 template<class Archive>
00581 inline void save_construct_data(
00582     Archive & ar, const LinearSystem * t, const unsigned int file_version)
00583 {
00584 
00585     std::string archive_filename_lhs = ArchiveLocationInfo::GetArchiveDirectory() + "lhs.mat";
00586     std::string archive_filename_rhs = ArchiveLocationInfo::GetArchiveDirectory() + "rhs.vec";
00587     const unsigned size = t->GetSize();
00588     ar << size;
00589 
00590     PetscTools::DumpPetscObject(t->GetRhsVector(), archive_filename_rhs);
00591     PetscTools::DumpPetscObject(t->GetLhsMatrix(), archive_filename_lhs);
00592 
00593     //Is the matrix structurally symmetric?
00594     PetscTruth symm_set, is_symmetric;
00595     is_symmetric = PETSC_FALSE;
00596     //Note that the following call only changes is_symmetric when symm_set is true
00597     MatIsSymmetricKnown(t->GetLhsMatrix(), &symm_set, &is_symmetric);
00598     assert(symm_set == is_symmetric);
00599     ar << symm_set;
00600 }
00601 
00606 template<class Archive>
00607 inline void load_construct_data(
00608     Archive & ar, LinearSystem * t, const unsigned int file_version)
00609 {
00610      std::string archive_filename_lhs = ArchiveLocationInfo::GetArchiveDirectory() + "lhs.mat";
00611      std::string archive_filename_rhs = ArchiveLocationInfo::GetArchiveDirectory() + "rhs.vec";
00612 
00613      PetscInt size;
00614      ar >> size;
00615 
00616      Vec new_vec;
00617      PetscTools::ReadPetscObject(new_vec, archive_filename_rhs);
00618 
00619      Mat new_mat;
00620      PetscTools::ReadPetscObject(new_mat, archive_filename_lhs);
00621 
00622      //This has to occur after the call to MatLoad as the matrix does not exist until MatLoad is called.
00623      //The property will be copied & set correctly in the LinearSystem constructor.
00624      PetscTruth symm_set;
00625      ar >> symm_set;
00626      if (symm_set == PETSC_TRUE)
00627      {
00628         PetscMatTools::SetOption(new_mat, MAT_SYMMETRIC);
00629         PetscMatTools::SetOption(new_mat, MAT_SYMMETRY_ETERNAL);
00630      }
00631 
00632      ::new(t)LinearSystem(size, new_mat, new_vec);
00633 }
00634 }
00635 } // namespace ...
00636 
00637 #endif //_LINEARSYSTEM_HPP_

Generated on Tue May 31 14:31:42 2011 for Chaste by  doxygen 1.5.5