LinearSystem.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
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 
00034 #include "LinearSystem.hpp"
00035 #include "PetscException.hpp"
00036 #include <iostream>
00037 #include "OutputFileHandler.hpp"
00038 #include "PetscTools.hpp"
00039 #include <cassert>
00040 #include "HeartEventHandler.hpp"
00041 
00042 
00043 
00044 LinearSystem::LinearSystem(PetscInt lhsVectorSize, MatType matType)
00045    :mMatNullSpace(NULL),
00046     mDestroyMatAndVec(true),
00047     mKspIsSetup(false),
00048     mMatrixIsConstant(false),
00049     mTolerance(1e-6),
00050     mUseAbsoluteTolerance(false),
00051     mDirichletBoundaryConditionsVector(NULL)
00052 {
00053     VecCreate(PETSC_COMM_WORLD, &mRhsVector);
00054     VecSetSizes(mRhsVector, PETSC_DECIDE, lhsVectorSize);
00055     VecSetFromOptions(mRhsVector);
00056 
00057     PetscTools::SetupMat(mLhsMatrix, lhsVectorSize, lhsVectorSize, matType);
00058 
00059     mSize = lhsVectorSize;
00060 
00061     VecGetOwnershipRange(mRhsVector, &mOwnershipRangeLo, &mOwnershipRangeHi);
00062     
00065     strcpy(mKspType, "gmres");
00066     strcpy(mPcType, "jacobi");
00067 }
00068 
00076 LinearSystem::LinearSystem(Vec templateVector)
00077    :mMatNullSpace(NULL),
00078     mDestroyMatAndVec(true),
00079     mKspIsSetup(false),
00080     mMatrixIsConstant(false),
00081     mTolerance(1e-6),
00082     mUseAbsoluteTolerance(false),
00083     mDirichletBoundaryConditionsVector(NULL)
00084 {
00085     VecDuplicate(templateVector, &mRhsVector);
00086     VecGetSize(mRhsVector, &mSize);
00087     VecGetOwnershipRange(mRhsVector, &mOwnershipRangeLo, &mOwnershipRangeHi);
00088     PetscInt local_size = mOwnershipRangeHi - mOwnershipRangeLo;
00089 
00090     PetscTools::SetupMat(mLhsMatrix, mSize, mSize, (MatType) MATMPIAIJ, local_size, local_size);
00091 
00094     strcpy(mKspType, "gmres");
00095     strcpy(mPcType, "jacobi");
00096     
00097 }
00098 
00106 LinearSystem::LinearSystem(Vec residualVector, Mat jacobianMatrix)
00107     :mMatNullSpace(NULL),
00108     mDestroyMatAndVec(false),
00109     mKspIsSetup(false),
00110     mMatrixIsConstant(false),
00111     mTolerance(1e-6),
00112     mUseAbsoluteTolerance(false),
00113     mDirichletBoundaryConditionsVector(NULL)
00114 {
00115     assert(residualVector || jacobianMatrix);
00116     mRhsVector = residualVector;
00117     mLhsMatrix = jacobianMatrix;
00118 
00119     PetscInt mat_size=0, vec_size=0;
00120     if (mRhsVector)
00121     {
00122         VecGetSize(mRhsVector, &vec_size);
00123         mSize = (unsigned)vec_size;
00124         VecGetOwnershipRange(mRhsVector, &mOwnershipRangeLo, &mOwnershipRangeHi);
00125     }
00126     if (mLhsMatrix)
00127     {
00128         PetscInt mat_cols;
00129         MatGetSize(mLhsMatrix, &mat_size, &mat_cols);
00130         assert(mat_size == mat_cols);
00131         mSize = (unsigned)mat_size;
00132         MatGetOwnershipRange(mLhsMatrix, &mOwnershipRangeLo, &mOwnershipRangeHi);
00133     }
00134     assert(!mRhsVector || !mLhsMatrix || vec_size == mat_size);
00135 
00138     strcpy(mKspType, "gmres");
00139     strcpy(mPcType, "jacobi");
00140     
00141 }
00142 
00143 LinearSystem::~LinearSystem()
00144 {
00145     if (mDestroyMatAndVec)
00146     {
00147         VecDestroy(mRhsVector);
00148         MatDestroy(mLhsMatrix);
00149     }
00150     
00151     if (mMatNullSpace)
00152     {
00153         MatNullSpaceDestroy(mMatNullSpace);
00154     }
00155     
00156     if (mKspIsSetup)
00157     {
00158         KSPDestroy(mKspSolver);
00159     }
00160     
00161     if(mDirichletBoundaryConditionsVector)
00162     {
00163         VecDestroy(mDirichletBoundaryConditionsVector);
00164     }    
00165 }
00166 
00167 void LinearSystem::SetMatrixElement(PetscInt row, PetscInt col, double value)
00168 {
00169     if (row >= mOwnershipRangeLo && row < mOwnershipRangeHi)
00170     {
00171         MatSetValue(mLhsMatrix, row, col, value, INSERT_VALUES);
00172     }
00173 }
00174 
00175 void LinearSystem::AddToMatrixElement(PetscInt row, PetscInt col, double value)
00176 {
00177     if (row >= mOwnershipRangeLo && row < mOwnershipRangeHi)
00178     {
00179         MatSetValue(mLhsMatrix, row, col, value, ADD_VALUES);
00180     }
00181 }
00182 
00183 
00184 void LinearSystem::AssembleFinalLinearSystem()
00185 {
00186     AssembleFinalLhsMatrix();
00187     AssembleRhsVector();
00188 }
00189 
00190 
00191 void LinearSystem::AssembleIntermediateLinearSystem()
00192 {
00193     AssembleIntermediateLhsMatrix();
00194     AssembleRhsVector();
00195 }
00196 
00197 void LinearSystem::AssembleFinalLhsMatrix()
00198 {
00199     MatAssemblyBegin(mLhsMatrix, MAT_FINAL_ASSEMBLY);
00200     MatAssemblyEnd(mLhsMatrix, MAT_FINAL_ASSEMBLY);
00201 }
00202 
00203 
00204 void LinearSystem::AssembleIntermediateLhsMatrix()
00205 {
00206     MatAssemblyBegin(mLhsMatrix, MAT_FLUSH_ASSEMBLY);
00207     MatAssemblyEnd(mLhsMatrix, MAT_FLUSH_ASSEMBLY);
00208 }
00209 
00210 void LinearSystem::AssembleRhsVector()
00211 {
00212     VecAssemblyBegin(mRhsVector);
00213     VecAssemblyEnd(mRhsVector);
00214 }
00215 
00216 
00217 
00218 void LinearSystem::SetRhsVectorElement(PetscInt row, double value)
00219 {
00220     if (row >= mOwnershipRangeLo && row < mOwnershipRangeHi)
00221     {
00222         VecSetValues(mRhsVector, 1, &row, &value, INSERT_VALUES);
00223     }
00224 }
00225 
00226 void LinearSystem::AddToRhsVectorElement(PetscInt row, double value)
00227 {
00228     if (row >= mOwnershipRangeLo && row < mOwnershipRangeHi)
00229     {
00230         VecSetValues(mRhsVector, 1, &row, &value, ADD_VALUES);
00231     }
00232 }
00233 
00234 void LinearSystem::DisplayMatrix()
00235 {
00236     MatView(mLhsMatrix,PETSC_VIEWER_STDOUT_WORLD);
00237 }
00238 
00239 void LinearSystem::DisplayRhs()
00240 {
00241     VecView(mRhsVector,PETSC_VIEWER_STDOUT_WORLD);
00242 }
00243 
00244 void LinearSystem::SetMatrixRow(PetscInt row, double value)
00245 {
00246     if (row >= mOwnershipRangeLo && row < mOwnershipRangeHi)
00247     {
00248         PetscInt rows, cols;
00249         MatGetSize(mLhsMatrix, &rows, &cols);
00250         for (PetscInt i=0; i<cols; i++)
00251         {
00252             this->SetMatrixElement(row, i, value);
00253         }
00254     }
00255 }
00256 
00257 void LinearSystem::ZeroMatrixRow(PetscInt row)
00258 {
00259     MatAssemblyBegin(mLhsMatrix, MAT_FINAL_ASSEMBLY);
00260     MatAssemblyEnd(mLhsMatrix, MAT_FINAL_ASSEMBLY);
00261     double diag_zero=0.0;
00262     // MatZeroRows allows a non-zero value to be placed on the diagonal
00263     // diag_zero is the value to put in the diagonal
00264 
00265 #if (PETSC_VERSION_MINOR == 2) //Old API
00266     IS is;
00267     ISCreateGeneral(PETSC_COMM_WORLD,1,&row,&is);
00268     MatZeroRows(mLhsMatrix, is, &diag_zero);
00269     ISDestroy(is);
00270 #else
00271 
00272     MatZeroRows(mLhsMatrix, 1, &row, diag_zero);
00273 #endif
00274 
00275 }
00276 
00283 void LinearSystem::ZeroMatrixColumn(PetscInt col)
00284 {
00285     MatAssemblyBegin(mLhsMatrix, MAT_FINAL_ASSEMBLY);
00286     MatAssemblyEnd(mLhsMatrix, MAT_FINAL_ASSEMBLY);
00287 
00288 //#if (PETSC_VERSION_MINOR == 2) //Old API
00289    // hello Joe.
00290    // Hello
00291 //#else
00292     // determine which rows in this column are non-zero (and
00293     // therefore need to be zeroed)
00294     std::vector<unsigned> nonzero_rows;
00295     for (PetscInt row = mOwnershipRangeLo; row < mOwnershipRangeHi; row++)
00296     {
00297         if (GetMatrixElement(row, col) != 0.0)
00298         {
00299             nonzero_rows.push_back(row);
00300         }    
00301     }
00302 
00303     // set those rows to be zero by calling MatSetValues    
00304     unsigned size = nonzero_rows.size();
00305     PetscInt* rows = new PetscInt[size];
00306     PetscInt cols[1];
00307     double* zeros = new double[size];
00308     
00309     cols[0] = col;
00310     
00311     for (unsigned i=0; i<size; i++)
00312     {
00313         rows[i] = nonzero_rows[i];
00314         zeros[i] = 0.0;        
00315     }
00316     
00317     MatSetValues(mLhsMatrix, size, rows, 1, cols, zeros, INSERT_VALUES);
00318     delete [] rows;
00319     delete [] zeros;
00320 //#endif
00321 }
00322 
00323 void LinearSystem::ZeroRhsVector()
00324 {
00325     double *p_rhs_vector_array;
00326     VecGetArray(mRhsVector, &p_rhs_vector_array);
00327     for (PetscInt local_index=0; local_index<mOwnershipRangeHi - mOwnershipRangeLo; local_index++)
00328     {
00329         p_rhs_vector_array[local_index]=0.0;
00330     }
00331     VecRestoreArray(mRhsVector, &p_rhs_vector_array);
00332 }
00333 
00334 
00335 void LinearSystem::ZeroLhsMatrix()
00336 {
00337     MatZeroEntries(mLhsMatrix);
00338 }
00339 
00340 void LinearSystem::ZeroLinearSystem()
00341 {
00342     ZeroRhsVector();
00343     ZeroLhsMatrix();
00344 }
00345 
00346 
00347 
00348 unsigned LinearSystem::GetSize()
00349 {
00350     return (unsigned) mSize;
00351 }
00352 
00353 void LinearSystem::SetNullBasis(Vec nullBasis[], unsigned numberOfBases)
00354 {
00355     PETSCEXCEPT( MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_FALSE, numberOfBases, nullBasis, &mMatNullSpace) );
00356 }
00357 
00361 void LinearSystem::GetOwnershipRange(PetscInt &lo, PetscInt &hi)
00362 {
00363     lo = mOwnershipRangeLo;
00364     hi = mOwnershipRangeHi;
00365 }
00366 
00371 double LinearSystem::GetMatrixElement(PetscInt row, PetscInt col)
00372 {
00373     assert(mOwnershipRangeLo <= row && row < mOwnershipRangeHi);
00374     PetscInt row_as_array[1];
00375     row_as_array[0] = row;
00376     PetscInt col_as_array[1];
00377     col_as_array[0] = col;
00378 
00379     double ret_array[1];
00380 
00381     MatGetValues(mLhsMatrix, 1, row_as_array, 1, col_as_array, ret_array);
00382 
00383     return ret_array[0];
00384 }
00385 
00390 double LinearSystem::GetRhsVectorElement(PetscInt row)
00391 {
00392     assert(mOwnershipRangeLo <= row && row < mOwnershipRangeHi);
00393 
00394     double *p_rhs_vector;
00395     PetscInt local_index=row-mOwnershipRangeLo;
00396     VecGetArray(mRhsVector, &p_rhs_vector);
00397     double answer=p_rhs_vector[local_index];
00398     VecRestoreArray(mRhsVector, &p_rhs_vector);
00399 
00400     return answer;
00401 }
00402 
00406 Vec& LinearSystem::rGetRhsVector()
00407 {
00408     return mRhsVector;
00409 }
00410 
00414 Mat& LinearSystem::rGetLhsMatrix()
00415 {
00416     return mLhsMatrix;
00417 }
00418 
00424 Vec& LinearSystem::rGetDirichletBoundaryConditionsVector()
00425 {
00426     return mDirichletBoundaryConditionsVector;   
00427 }
00428 
00431 void LinearSystem::SetMatrixIsSymmetric()
00432 {
00433     MatSetOption(mLhsMatrix, MAT_SYMMETRIC);
00434     MatSetOption(mLhsMatrix, MAT_SYMMETRY_ETERNAL);
00435 }
00436 void LinearSystem::SetMatrixIsConstant(bool matrixIsConstant)
00437 {
00438     mMatrixIsConstant=matrixIsConstant;
00439 }
00440 
00441 void LinearSystem::SetRelativeTolerance(double relativeTolerance)
00442 {
00443     mTolerance=relativeTolerance;
00444     mUseAbsoluteTolerance=false;
00445     if (mKspIsSetup)
00446     {
00447         KSPSetTolerances(mKspSolver, mTolerance, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
00448     }
00449 
00450 }
00451 
00452 void LinearSystem::SetAbsoluteTolerance(double absoluteTolerance)
00453 {
00454     mTolerance=absoluteTolerance;
00455     mUseAbsoluteTolerance=true;
00456     if (mKspIsSetup)
00457     {
00458         KSPSetTolerances(mKspSolver, DBL_EPSILON, mTolerance, PETSC_DEFAULT, PETSC_DEFAULT);
00459     }
00460 
00461 }
00462 
00463 void LinearSystem::SetKspType(const char *kspType)
00464 {
00465     strcpy(mKspType, kspType);
00466     if (mKspIsSetup)
00467     {
00468         KSPSetType(mKspSolver, kspType);
00469         KSPSetFromOptions(mKspSolver);
00470     }
00471 }
00472 
00473 void LinearSystem::SetPcType(const char *pcType)
00474 {
00475     strcpy(mPcType, pcType);
00476     if (mKspIsSetup)
00477     {
00478         PC prec;
00479         KSPGetPC(mKspSolver, &prec);
00480         PCSetType(prec, pcType);
00481         KSPSetFromOptions(mKspSolver);
00482     }
00483 }
00484 
00485 Vec LinearSystem::Solve(Vec lhsGuess)
00486 {
00487     /* The following lines are very useful for debugging
00488      *    MatView(mLhsMatrix,    PETSC_VIEWER_STDOUT_WORLD);
00489      *    VecView(mRhsVector,    PETSC_VIEWER_STDOUT_WORLD);
00490      */
00491     //Double check that the non-zero pattern hasn't changed
00492     MatInfo mat_info;
00493     MatGetInfo(mLhsMatrix, MAT_GLOBAL_SUM, &mat_info);
00494 
00495     if (!mKspIsSetup)
00496     {
00497         mNonZerosUsed=mat_info.nz_used;
00498         //MatNorm(mLhsMatrix, NORM_FROBENIUS, &mMatrixNorm);
00499         PC prec; //Type of pre-conditioner
00500 
00501         KSPCreate(PETSC_COMM_WORLD, &mKspSolver);
00502         //See
00503         //http://www-unix.mcs.anl.gov/petsc/petsc-2/snapshots/petsc-current/docs/manualpages/KSP/KSPSetOperators.html
00504         //The preconditioner flag (last argument) in the following calls says
00505         //how to reuse the preconditioner on subsequent iterations
00506         if (mMatrixIsConstant)
00507         {
00508             KSPSetOperators(mKspSolver, mLhsMatrix, mLhsMatrix, SAME_PRECONDITIONER);
00509         }
00510         else
00511         {
00512             KSPSetOperators(mKspSolver, mLhsMatrix, mLhsMatrix, SAME_NONZERO_PATTERN);
00513         }
00514 
00515         // Set either absolute or relative tolerance of the KSP solver.
00516         // The default is to use relative tolerance (1e-6)
00517         if (mUseAbsoluteTolerance)
00518         {
00519             KSPSetTolerances(mKspSolver, DBL_EPSILON, mTolerance, PETSC_DEFAULT, PETSC_DEFAULT);
00520         }
00521         else
00522         {
00523             KSPSetTolerances(mKspSolver, mTolerance, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
00524         }
00525 
00526         // set ksp and pc types
00527         KSPSetType(mKspSolver, mKspType);
00528         KSPGetPC(mKspSolver, &prec);
00529  
00530         
00531         // Turn off pre-conditioning if the system size is very small
00532        if (mSize <= 4)
00533         {
00534             PCSetType(prec, PCNONE);
00535         }
00536         else
00537         {
00538             PCSetType(prec, mPcType);
00539         }
00540 
00541         if (mMatNullSpace)
00542         {
00543             PETSCEXCEPT( KSPSetNullSpace(mKspSolver, mMatNullSpace) );
00544         }
00545 
00546         if (lhsGuess)
00547         {
00548             //Assume that the user of this method will always be kind enough
00549             //to give us a reasonable guess.
00550             KSPSetInitialGuessNonzero(mKspSolver,PETSC_TRUE);
00551         }
00552 
00553         KSPSetFromOptions(mKspSolver);
00554         KSPSetUp(mKspSolver);
00555 
00556         mKspIsSetup = true;
00557     }
00558     else
00559     {
00560         #define COVERAGE_IGNORE
00561         if (mNonZerosUsed!=mat_info.nz_used)
00562         {
00563             EXCEPTION("LinearSystem doesn't allow the non-zero pattern of a matrix to change. (I think you changed it).");
00564         }
00565 //        PetscScalar norm;
00566 //        MatNorm(mLhsMatrix, NORM_FROBENIUS, &norm);
00567 //        if (fabs(norm - mMatrixNorm) > 0)
00568 //        {
00569 //            EXCEPTION("LinearSystem doesn't allow the matrix norm to change");
00570 //        }
00571         #undef COVERAGE_IGNORE
00572     }
00573 
00574     // Create solution vector
00576     Vec lhs_vector;
00577     VecDuplicate(mRhsVector, &lhs_vector);//Sets the same size (doesn't copy)
00578     if (lhsGuess)
00579     {
00580         VecCopy(lhsGuess, lhs_vector);
00581         //VecZeroEntries(lhs_vector);
00582     }
00583 //    //Double check that the mRhsVector contains sensible values
00584 //    double *p_rhs, *p_guess;
00585 //    VecGetArray(mRhsVector, &p_rhs);
00586 //    VecGetArray(lhsGuess, &p_guess);
00587 //    for (int global_index=mOwnershipRangeLo; global_index<mOwnershipRangeHi; global_index++)
00588 //    {
00589 //        int local_index = global_index - mOwnershipRangeLo;
00590 //        assert(!isnan(p_rhs[local_index]));
00591 //        assert(!isnan(p_guess[local_index]));
00592 //        if (p_rhs[local_index] != p_rhs[local_index])
00593 //        {
00594 //            std::cout << "********* PETSc nan\n";
00595 //            assert(0);
00596 //        }
00597 //    }
00598 //    std::cout << "b[0] = " << p_rhs[0] << ", guess[0] = " << p_guess[0] << "\n";
00599 //    VecRestoreArray(mRhsVector, &p_rhs);
00600 //    VecRestoreArray(lhsGuess, &p_guess);
00601 //    // Test A*guess
00602 //    Vec temp;
00603 //    VecDuplicate(mRhsVector, &temp);
00604 //    MatMult(mLhsMatrix, lhs_vector, temp);
00605 //    double *p_temp;
00606 //    VecGetArray(temp, &p_temp);
00607 //    std::cout << "temp[0] = " << p_temp[0] << "\n";
00608 //    VecRestoreArray(temp, &p_temp);
00609 //    VecDestroy(temp);
00610 //    // Dump the matrix to file
00611 //    PetscViewer viewer;
00612 //    PetscViewerASCIIOpen(PETSC_COMM_WORLD,"mat.output",&viewer);
00613 //    MatView(mLhsMatrix, viewer);
00614 //    PetscViewerFlush(viewer);
00615 //    PetscViewerDestroy(viewer);
00616 //    // Dump the rhs vector to file
00617 //    PetscViewerASCIIOpen(PETSC_COMM_WORLD,"vec.output",&viewer);
00618 //    PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
00619 //    VecView(mRhsVector, viewer);
00620 //    PetscViewerFlush(viewer);
00621 //    PetscViewerDestroy(viewer);
00622     
00623     
00624     try
00625     {
00626         HeartEventHandler::BeginEvent(HeartEventHandler::SOLVE_LINEAR_SYSTEM);
00627         PETSCEXCEPT(KSPSolve(mKspSolver, mRhsVector, lhs_vector));
00628         HeartEventHandler::EndEvent(HeartEventHandler::SOLVE_LINEAR_SYSTEM);
00629 
00630         // Check that solver converged and throw if not
00631         KSPConvergedReason reason;
00632         KSPGetConvergedReason(mKspSolver, &reason);
00633         KSPEXCEPT(reason);
00634     }
00635     catch (const Exception& e)
00636     {
00637         // Destroy solution vector on error to avoid memory leaks
00638         VecDestroy(lhs_vector);
00639         throw e;
00640     }
00641 
00642     return lhs_vector;
00643 }

Generated on Wed Mar 18 12:51:52 2009 for Chaste by  doxygen 1.5.5