00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00263
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
00289
00290
00291
00292
00293
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
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
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
00488
00489
00490
00491
00492 MatInfo mat_info;
00493 MatGetInfo(mLhsMatrix, MAT_GLOBAL_SUM, &mat_info);
00494
00495 if (!mKspIsSetup)
00496 {
00497 mNonZerosUsed=mat_info.nz_used;
00498
00499 PC prec;
00500
00501 KSPCreate(PETSC_COMM_WORLD, &mKspSolver);
00502
00503
00504
00505
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
00516
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
00527 KSPSetType(mKspSolver, mKspType);
00528 KSPGetPC(mKspSolver, &prec);
00529
00530
00531
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
00549
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
00566
00567
00568
00569
00570
00571 #undef COVERAGE_IGNORE
00572 }
00573
00574
00576 Vec lhs_vector;
00577 VecDuplicate(mRhsVector, &lhs_vector);
00578 if (lhsGuess)
00579 {
00580 VecCopy(lhsGuess, lhs_vector);
00581
00582 }
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621
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
00631 KSPConvergedReason reason;
00632 KSPGetConvergedReason(mKspSolver, &reason);
00633 KSPEXCEPT(reason);
00634 }
00635 catch (const Exception& e)
00636 {
00637
00638 VecDestroy(lhs_vector);
00639 throw e;
00640 }
00641
00642 return lhs_vector;
00643 }