Flow123d  master-f44eb46
linsys_PETSC.cc
Go to the documentation of this file.
1 /*!
2  *
3  * Copyright (C) 2015 Technical University of Liberec. All rights reserved.
4  *
5  * This program is free software; you can redistribute it and/or modify it under
6  * the terms of the GNU General Public License version 3 as published by the
7  * Free Software Foundation. (http://www.gnu.org/licenses/gpl-3.0.en.html)
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12  *
13  *
14  * @file linsys_PETSC.cc
15  * @brief Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construction
16  * @author Jakub Sistek
17  */
18 
19 // derived from base linsys
20 #include "la/linsys_PETSC.hh"
21 #include "petscvec.h"
22 #include "petscksp.h"
23 #include "petscmat.h"
24 #include "system/sys_profiler.hh"
25 #include "system/system.hh"
26 
27 
28 //#include <boost/bind.hpp>
29 
30 namespace it = Input::Type;
31 
33  return it::Record("Petsc", "PETSc solver settings.\n It provides interface to various PETSc solvers. The convergence criteria is:\n"
34  "```\n"
35  "norm( res_i ) < max( norm( res_0 ) * r_tol, a_tol )\n"
36  "```\n"
37  "where ```res_i``` is the residuum vector after i-th iteration of the solver and ```res_0``` is the estimate of the norm of the initial residual. "
38  "If the initial guess of the solution is provided (usually only for transient equations) the residual of this estimate is used, "
39  "otherwise the norm of preconditioned RHS is used. "
40  "The default norm is (($L_2$)) norm of preconditioned residual: (($ P^{-1}(Ax-b)$)), usage of other norm may be prescribed using the 'option' key. "
41  "See also PETSc documentation for KSPSetNormType.")
43  .declare_key("r_tol", it::Double(0.0, 1.0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
44  "If not, we use the value 1.0e-7."),
45  "Residual tolerance relative to the initial error.")
46  .declare_key("a_tol", it::Double(0.0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
47  "If not, we use the value 1.0e-11."),
48  "Absolute residual tolerance.")
49  .declare_key("d_tol", it::Double(0.0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
50  "If not, we use the value 10000."),
51  "Tolerance for divergence.")
52  .declare_key("max_it", it::Integer(0), it::Default::read_time("Default value is set by the nonlinear solver or the equation. "
53  "If not, we use the value 1000."),
54  "Maximum number of outer iterations of the linear solver.")
55  .declare_key("options", it::String(), it::Default("\"\""), "This options is passed to PETSC to create a particular KSP (Krylov space method).\n"
56  "If the string is left empty (by default), the internal default options is used.")
57  .close();
58 }
59 
60 
62 
63 
64 LinSys_PETSC::LinSys_PETSC( const Distribution * rows_ds, const std::string &params)
65  : LinSys( rows_ds ),
66  params_(params),
67  init_guess_nonzero(false),
68  matrix_(0)
69 {
70  // create PETSC vectors:
71  PetscErrorCode ierr;
72  // rhs
73  v_rhs_= new double[ rows_ds_->lsize() + 1 ];
74  ierr = VecCreateMPIWithArray( comm_, 1, rows_ds_->lsize(), PETSC_DECIDE, v_rhs_, &rhs_ ); CHKERRV( ierr );
75  ierr = VecZeroEntries( rhs_ ); CHKERRV( ierr );
76  VecDuplicate(rhs_, &residual_);
77 
78  matrix_ = NULL;
79  solution_precision_ = std::numeric_limits<double>::infinity();
80  matrix_changed_ = true;
81  rhs_changed_ = true;
82 }
83 
85  : LinSys(other), params_(other.params_), v_rhs_(NULL), solution_precision_(other.solution_precision_)
86 {
87  MatCopy(other.matrix_, matrix_, DIFFERENT_NONZERO_PATTERN);
88  VecCopy(other.rhs_, rhs_);
89  VecCopy(other.on_vec_, on_vec_);
90  VecCopy(other.off_vec_, off_vec_);
91 }
92 
93 void LinSys_PETSC::set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it)
94 {
95  if (! in_rec_.is_empty()) {
96  // input record is set
97  r_tol_ = in_rec_.val<double>("r_tol", r_tol);
98  a_tol_ = in_rec_.val<double>("a_tol", a_tol);
99  d_tol_ = in_rec_.val<double>("d_tol", d_tol);
100  max_it_ = in_rec_.val<unsigned int>("max_it", max_it);
101  } else {
102  r_tol_ = r_tol;
103  a_tol_ = a_tol;
104  d_tol_ = d_tol;
105  max_it_ = max_it;
106 
107  }
108 }
109 
110 
112 {
113  PetscErrorCode ierr;
114 
115  ierr = VecCreateMPI( comm_, rows_ds_->lsize(), PETSC_DECIDE, &(on_vec_) ); CHKERRV( ierr );
116  ierr = VecDuplicate( on_vec_, &(off_vec_) ); CHKERRV( ierr );
117  status_ = ALLOCATE;
118 }
119 
121 {
122  switch ( status_ ) {
123  case ALLOCATE:
124  this->preallocate_matrix( );
125  break;
126  case INSERT:
127  this->finish_assembly( MAT_FLUSH_ASSEMBLY );
128  break;
129  case ADD:
130  case DONE:
131  break;
132  default:
133  ASSERT_PERMANENT(false).error("Can not set values. Matrix is not preallocated.\n");
134  }
135  status_ = ADD;
136 }
137 
139 {
140  switch ( status_ ) {
141  case ALLOCATE:
142  this->preallocate_matrix();
143  break;
144  case ADD:
145  this->finish_assembly( MAT_FLUSH_ASSEMBLY );
146  break;
147  case INSERT:
148  case DONE:
149  break;
150  default:
151  ASSERT_PERMANENT(false).error("Can not set values. Matrix is not preallocated.\n");
152  }
153  status_ = INSERT;
154 }
155 
156 void LinSys_PETSC::mat_set_values( int nrow, int *rows, int ncol, int *cols, double *vals )
157 {
158  // here vals would need to be converted from double to PetscScalar if it was ever something else than double :-)
159  switch (status_) {
160  case INSERT:
161  case ADD:
162  chkerr(MatSetValues(matrix_,nrow,rows,ncol,cols,vals,(InsertMode)status_));
163  break;
164  case ALLOCATE:
165  this->preallocate_values(nrow,rows,ncol,cols);
166  break;
167  default: DebugOut() << "LS SetValues with non allowed insert mode.\n";
168  }
169 
170  matrix_changed_ = true;
171 }
172 
173 void LinSys_PETSC::rhs_set_values( int nrow, int *rows, double *vals )
174 {
175  PetscErrorCode ierr;
176 
177  switch (status_) {
178  case INSERT:
179  case ADD:
180  ierr = VecSetValues(rhs_,nrow,rows,vals,(InsertMode)status_); CHKERRV( ierr );
181  break;
182  case ALLOCATE:
183  break;
184  default: ASSERT_PERMANENT(false).error("LinSys's status disallow set values.\n");
185  }
186 
187  rhs_changed_ = true;
188 }
189 
190 void LinSys_PETSC::preallocate_values(int nrow,int *rows,int ncol,int *cols)
191 {
192  int i,j;
193  int col;
194  PetscInt row;
195 
196  for (i=0; i<nrow; i++) {
197  row=rows[i];
198  for(j=0; j<ncol; j++) {
199  col = cols[j];
200  if (rows_ds_->get_proc(row) == rows_ds_->get_proc(col))
201  VecSetValue(on_vec_,row,1.0,ADD_VALUES);
202  else
203  VecSetValue(off_vec_,row,1.0,ADD_VALUES);
204  }
205  }
206 }
207 
209 {
210  ASSERT_EQ(status_, ALLOCATE).error("Linear system has to be in ALLOCATE status.");
211 
212  PetscErrorCode ierr;
213  PetscInt *on_nz, *off_nz;
214  PetscScalar *on_array, *off_array;
215 
216  // assembly and get values from counting vectors, destroy them
217  VecAssemblyBegin(on_vec_);
218  VecAssemblyBegin(off_vec_);
219 
220  on_nz = new PetscInt[ rows_ds_->lsize() ];
221  off_nz = new PetscInt[ rows_ds_->lsize() ];
222 
223  VecAssemblyEnd(on_vec_);
224  VecAssemblyEnd(off_vec_);
225 
226  VecGetArray( on_vec_, &on_array );
227  VecGetArray( off_vec_, &off_array );
228 
229  for ( unsigned int i=0; i<rows_ds_->lsize(); i++ ) {
230  on_nz[i] = std::min( rows_ds_->lsize(), static_cast<uint>( on_array[i]+0.1 ) ); // small fraction to ensure correct rounding
231  off_nz[i] = std::min( rows_ds_->size() - rows_ds_->lsize(), static_cast<uint>( off_array[i]+0.1 ) );
232  }
233 
234  VecRestoreArray(on_vec_,&on_array);
235  VecRestoreArray(off_vec_,&off_array);
236  VecDestroy(&on_vec_);
237  VecDestroy(&off_vec_);
238 
239  // create PETSC matrix with preallocation
240  if (matrix_ != NULL)
241  {
242  chkerr(MatDestroy(&matrix_));
243  }
244  ierr = MatCreateAIJ(PETSC_COMM_WORLD, rows_ds_->lsize(), rows_ds_->lsize(), PETSC_DETERMINE, PETSC_DETERMINE,
245  0, on_nz, 0, off_nz, &matrix_); CHKERRV( ierr );
246 
247  if (symmetric_) MatSetOption(matrix_, MAT_SYMMETRIC, PETSC_TRUE);
248  MatSetOption(matrix_, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
249 
250  // This option is used in order to assembly larger local matrices with own non-zero structure.
251  // Zero entries are ignored so we must prevent adding exact zeroes.
252  // Add LocalSystem::almost_zero for entries that should not be eliminated.
253  MatSetOption(matrix_, MAT_IGNORE_ZERO_ENTRIES, PETSC_TRUE);
254 
255 
256 
257  delete[] on_nz;
258  delete[] off_nz;
259 }
260 
262 {
263  MatAssemblyType assemblyType = MAT_FINAL_ASSEMBLY;
264  this->finish_assembly( assemblyType );
265 }
266 
267 void LinSys_PETSC::finish_assembly( MatAssemblyType assembly_type )
268 {
269  PetscErrorCode ierr;
270 
271  if (status_ == ALLOCATE) {
272  WarningOut() << "Finalizing linear system without setting values.\n";
273  this->preallocate_matrix();
274  }
275  ierr = MatAssemblyBegin(matrix_, assembly_type); CHKERRV( ierr );
276  ierr = VecAssemblyBegin(rhs_); CHKERRV( ierr );
277  ierr = MatAssemblyEnd(matrix_, assembly_type); CHKERRV( ierr );
278  ierr = VecAssemblyEnd(rhs_); CHKERRV( ierr );
279 
280  if (assembly_type == MAT_FINAL_ASSEMBLY) status_ = DONE;
281 
282  //PetscViewerPushFormat(PETSC_VIEWER_STDOUT_SELF, PETSC_VIEWER_ASCII_INDEX);
283  //MatView(matrix_, PETSC_VIEWER_STDOUT_SELF);
284  //VecView(rhs_, PETSC_VIEWER_STDOUT_SELF);
285  //this->view();
286 
287  matrix_changed_ = true;
288  rhs_changed_ = true;
289 }
290 
291 void LinSys_PETSC::apply_constrains( double scalar )
292 {
293  PetscErrorCode ierr;
294 
295  // check that system matrix is assembled
296  ASSERT_EQ(status_, DONE).error("System matrix and right-hand side are not assembled when applying constraints." );
297 
298  // number of constraints
299  PetscInt numConstraints = static_cast<PetscInt>( constraints_.size() );
300 
301  // Additional multiplier for numerical reasons (criterion to be established)
302  const PetscScalar diagScalar = static_cast<PetscScalar>( scalar );
303 
304  std::vector<PetscInt> globalDofs;
306 
307  // Constraint iterators
308  ConstraintVec_::const_iterator cIter = constraints_.begin( );
309  ConstraintVec_::const_iterator cEnd = constraints_.end( );
310  // collect global dof indices and the correpsonding values
311  for ( ; cIter != cEnd; ++cIter ) {
312  globalDofs.push_back( static_cast<PetscInt>( cIter -> first ) );
313  values.push_back( static_cast<PetscScalar>( cIter -> second ) * diagScalar );
314  }
315 
316  // prepare pointers to be passed to PETSc
317  PetscInt * globalDofPtr = this->makePetscPointer_( globalDofs );
318  PetscScalar * valuePtr = this->makePetscPointer_( values );
319 
320  // set matrix rows to zero
321  ierr = MatZeroRows( matrix_, numConstraints, globalDofPtr, diagScalar, PETSC_NULL, PETSC_NULL ); CHKERRV( ierr );
322  matrix_changed_ = true;
323 
324  // set RHS entries to values (crashes if called with NULL pointers)
325  if ( numConstraints ) {
326  ierr = VecSetValues( rhs_, numConstraints, globalDofPtr, valuePtr, INSERT_VALUES ); CHKERRV( ierr );
327  rhs_changed_ = true;
328  }
329 
330  // perform communication in the rhs vector
331  ierr = VecAssemblyBegin( rhs_ ); CHKERRV( ierr );
332  ierr = VecAssemblyEnd( rhs_ ); CHKERRV( ierr );
333 }
334 
335 
337 {
338  init_guess_nonzero = set_nonzero;
339 }
340 
341 
343 {
344 
345  const char *petsc_dflt_opt;
346  int nits;
347 
348  // -mat_no_inode ... inodes are usefull only for
349  // vector problems e.g. MH without Schur complement reduction
350 
351  /* Comment to PETSc options:
352  *
353  * -ksp_diagonal_scale scales the matrix before solution, while -ksp_diagonal_scale_fix just fixes the scaling after solution
354  * -pc_asm_type basic enforces classical Schwartz method, which seems more stable for positive definite systems.
355  * The default 'restricted' probably violates s.p.d. structure, many tests fail.
356  */
357  if (rows_ds_->np() > 1) {
358  // parallel setting
359  if (this->is_positive_definite())
360  petsc_dflt_opt="-ksp_type cg -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_type asm -pc_asm_type basic -pc_asm_overlap 4 -sub_pc_type icc -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
361  //petsc_dflt_opt="-ksp_type bcgs -ksp_diagonal_scale_fix -pc_type asm -pc_asm_overlap 4 -sub_pc_type ilu -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
362  else
363  petsc_dflt_opt="-ksp_type bcgs -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_type asm -pc_asm_overlap 4 -sub_pc_type ilu -sub_pc_factor_levels 3 -sub_pc_factor_fill 6.0";
364 
365  }
366  else {
367  // serial setting
368  if (this->is_positive_definite())
369  petsc_dflt_opt="-ksp_type cg -pc_type icc -pc_factor_levels 3 -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
370  //petsc_dflt_opt="-ksp_type bcgs -pc_type ilu -pc_factor_levels 5 -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
371  else
372  petsc_dflt_opt="-ksp_type bcgs -pc_type ilu -pc_factor_levels 5 -ksp_diagonal_scale -ksp_diagonal_scale_fix -pc_factor_fill 6.0";
373  }
374 
375  if (params_ == "") params_ = petsc_dflt_opt;
376  LogOut().fmt("inserting petsc options: {}\n",params_.c_str());
377 
378  // now takes an optional PetscOptions object as the first argument
379  // value NULL will preserve previous behaviour previous behavior.
380  PetscOptionsInsertString(NULL, params_.c_str()); // overwrites previous options values
381 
382  MatSetOption( matrix_, MAT_USE_INODES, PETSC_FALSE );
383 
384  chkerr(KSPCreate( comm_, &system ));
385  chkerr(KSPSetOperators(system, matrix_, matrix_));
386 
387 
388  // TODO take care of tolerances - shall we support both input file and command line petsc setting
389  chkerr(KSPSetTolerances(system, r_tol_, a_tol_, d_tol_,PETSC_DEFAULT));
390  chkerr(KSPSetTolerances(system, r_tol_, a_tol_, d_tol_, max_it_));
391  KSPSetFromOptions(system);
392  // We set the KSP flag set_initial_guess_nonzero
393  // unless KSP type is preonly.
394  // In such case PETSc fails (version 3.4.1)
395  if (init_guess_nonzero)
396  {
397  KSPType type;
398  KSPGetType(system, &type);
399  if (strcmp(type, KSPPREONLY) != 0)
400  KSPSetInitialGuessNonzero(system, PETSC_TRUE);
401  }
402 
403  {
404  START_TIMER("PETSC linear solver");
405  START_TIMER("PETSC linear iteration");
406  chkerr(KSPSolve(system, rhs_, solution_ ));
407  KSPGetConvergedReason(system,&reason);
408  KSPGetIterationNumber(system,&nits);
409  ADD_CALLS(nits);
410  }
411  // substitute by PETSc call for residual
412  VecNorm(rhs_, NORM_2, &residual_norm_);
413 
414  LogOut().fmt("convergence reason {}, number of iterations is {}\n", reason, nits);
415 
416  // get residual norm
417  KSPGetResidualNorm(system, &solution_precision_);
418 
419  // TODO: I do not understand this
420  //Profiler::instance()->set_timer_subframes("SOLVING MH SYSTEM", nits);
421 
422  chkerr(KSPDestroy(&system));
423 
424  return LinSys::SolveInfo(static_cast<int>(reason), static_cast<int>(nits));
425 
426 }
427 
428 void LinSys_PETSC::view(string text )
429 {
430  FilePath matFileName(text + "_flow123d_matrix.m",FilePath::FileType::output_file);
431  FilePath rhsFileName(text + "_flow123d_rhs.m",FilePath::FileType::output_file);
432  FilePath solFileName(text + "_flow123d_sol.m",FilePath::FileType::output_file);
433 
434  PetscViewer myViewer;
435 
436  if ( matrix_ != NULL ) {
437  PetscViewerASCIIOpen(comm_,((string)matFileName).c_str(),&myViewer);
438  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
439  MatView( matrix_, myViewer );
440  PetscViewerDestroy(&myViewer);
441  }
442  else
443  WarningOut() << "PetscViewer: the matrix of LinSys is not set.\n";
444 
445  if ( rhs_ != NULL ) {
446  PetscViewerASCIIOpen(comm_,((string)rhsFileName).c_str(),&myViewer);
447  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
448  VecView( rhs_, myViewer );
449  PetscViewerDestroy(&myViewer);
450  }
451  else
452  WarningOut() << "PetscViewer: the rhs of LinSys is not set.\n";
453 
454  if ( solution_ != NULL ) {
455  PetscViewerASCIIOpen(comm_,((string)solFileName).c_str(),&myViewer);
456  PetscViewerSetFormat(myViewer,PETSC_VIEWER_ASCII_MATLAB);
457  VecView( solution_, myViewer );
458  PetscViewerDestroy(&myViewer);
459  }
460  else
461  WarningOut() << "PetscViewer: the solution of LinSys is not set.\n";
462 }
463 
465 {
466  if (matrix_ != NULL) { chkerr(MatDestroy(&matrix_)); }
467  chkerr(VecDestroy(&rhs_));
468 
469  if (residual_ != NULL) chkerr(VecDestroy(&residual_));
470  if (v_rhs_ != NULL) delete[] v_rhs_;
471 }
472 
473 
474 
476 {
477  LinSys::set_from_input( in_rec );
478 
479  // PETSC specific parameters
480  // If parameters are specified in input file, they are used,
481  // otherwise keep settings provided in constructor of LinSys_PETSC.
482  std::string user_params = in_rec.val<string>("options");
483  if (user_params != "") params_ = user_params;
484 }
485 
486 
488 {
489  return solution_precision_;
490 }
491 
492 
494 {
495  MatMult(matrix_, solution_, residual_);
496  VecAXPY(residual_,-1.0, rhs_);
497  double residual_norm;
498  VecNorm(residual_, NORM_2, &residual_norm);
499  return residual_norm;
500 }
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
#define ASSERT_EQ(a, b)
Definition of comparative assert macro (EQual) only for debug mode.
Definition: asserts.hh:333
unsigned int lsize(int proc) const
get local size
unsigned int get_proc(unsigned int idx) const
get processor of the given index
unsigned int size() const
get global size
unsigned int np() const
get num of processors
Dedicated class for storing path to input and output files.
Definition: file_path.hh:54
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
bool is_empty() const
Definition: accessors.hh:365
const Ret val(const string &key) const
Class Input::Type::Default specifies default value of keys of a Input::Type::Record.
Definition: type_record.hh:61
static Default read_time(const std::string &description)
The factory function to make an default value that will be specified at the time when a key will be r...
Definition: type_record.hh:97
Class for declaration of the input data that are floating point numbers.
Definition: type_base.hh:534
Class for declaration of the integral input data.
Definition: type_base.hh:483
Record type proxy class.
Definition: type_record.hh:182
unsigned int size() const
Returns number of keys in the Record.
Definition: type_record.hh:602
virtual Record & derive_from(Abstract &parent)
Method to derive new Record from an AbstractRecord parent.
Definition: type_record.cc:196
Record & close() const
Close the Record for further declarations of keys.
Definition: type_record.cc:304
Record & declare_key(const string &key, std::shared_ptr< TypeBase > type, const Default &default_value, const string &description, TypeBase::attribute_map key_attributes=TypeBase::attribute_map())
Declares a new key of the Record.
Definition: type_record.cc:503
Class for declaration of the input data that are in string format.
Definition: type_base.hh:582
void finish_assembly() override
void start_add_assembly() override
void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals) override
void apply_constrains(double scalar=1.) override
double compute_residual() override
Mat matrix_
Petsc matrix of the problem.
void preallocate_matrix()
LinSys::SolveInfo solve() override
double solution_precision_
void preallocate_values(int nrow, int *rows, int ncol, int *cols)
bool init_guess_nonzero
flag for starting from nonzero guess
void start_insert_assembly() override
void set_from_input(const Input::Record in_rec) override
static const Input::Type::Record & get_input_type()
Definition: linsys_PETSC.cc:32
void set_initial_guess_nonzero(bool set_nonzero=true)
std::string params_
command-line-like options for the PETSc solver
LinSys_PETSC(const Distribution *rows_ds, const std::string &params="")
Definition: linsys_PETSC.cc:64
KSPConvergedReason reason
static const int registrar
Registrar of class to factory.
Vec rhs_
PETSc vector constructed with vx array.
Vec off_vec_
Vectors for counting non-zero entries in off-diagonal block.
Vec on_vec_
Vectors for counting non-zero entries in diagonal block.
T * makePetscPointer_(std::vector< T > &array)
void view(string text="") override
void rhs_set_values(int nrow, int *rows, double *vals) override
double * v_rhs_
local RHS array pointing to Vec rhs_
void start_allocation() override
double get_solution_precision() override
void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override
Definition: linsys_PETSC.cc:93
unsigned int max_it_
maximum number of iterations of linear solver
Definition: linsys.hh:668
double d_tol_
tolerance for divergence of linear solver
Definition: linsys.hh:667
MPI_Comm comm_
Definition: linsys.hh:670
virtual void set_from_input(const Input::Record in_rec)
Definition: linsys.hh:641
bool is_positive_definite()
Definition: linsys.hh:597
bool symmetric_
Definition: linsys.hh:678
Input::Record in_rec_
Definition: linsys.hh:697
ConstraintVec_ constraints_
Definition: linsys.hh:693
double a_tol_
absolute tolerance of linear solver
Definition: linsys.hh:666
@ INSERT
Definition: linsys.hh:97
@ DONE
Definition: linsys.hh:100
@ ADD
Definition: linsys.hh:98
@ ALLOCATE
Definition: linsys.hh:99
SetValuesMode status_
Set value status of the linear system.
Definition: linsys.hh:671
const Distribution * rows_ds_
final distribution of rows of MH matrix
Definition: linsys.hh:676
Vec solution_
PETSc vector constructed with vb array.
Definition: linsys.hh:686
double residual_norm_
local solution array pointing into Vec solution_
Definition: linsys.hh:691
double r_tol_
relative tolerance of linear solver
Definition: linsys.hh:665
bool matrix_changed_
true if the matrix was changed since the last solve
Definition: linsys.hh:683
static Input::Type::Abstract & get_input_type()
Definition: linsys.cc:29
bool rhs_changed_
true if the right hand side was changed since the last solve
Definition: linsys.hh:684
Solver based on the original PETSc solver using MPIAIJ matrix and succesive Schur complement construc...
#define WarningOut()
Macro defining 'warning' record of log.
Definition: logger.hh:278
#define LogOut()
Macro defining 'log' record of log.
Definition: logger.hh:281
#define DebugOut()
Macro defining 'debug' record of log.
Definition: logger.hh:284
unsigned int uint
#define ADD_CALLS(n_calls)
Increase number of calls in actual timer.
#define START_TIMER(tag)
Starts a timer with specified tag.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:142