Flow123d  master-f44eb46
linsys.hh
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.hh
15  * @brief Wrappers for linear systems based on MPIAIJ and MATIS format.
16  * @author Jan Brezina
17  *
18  * Linear system only keep together matrix, RHS and the solution vector.
19  *
20  */
21 
22 //=============================================================================
23 //
24 // LINER SYSTEM ROUTINES - linear system use wrapers of PETSc assembling routins
25 // in order to allow counting of allocation and filling of matrix to be done by the same code
26 //
27 // we want to allow allocations of nonlocal rows, to this end we count on- and off-processor
28 // members in an parallel Vector
29 //
30 //=============================================================================
31 
32 #ifndef LA_LINSYS_HH_
33 #define LA_LINSYS_HH_
34 
35 /**
36  * @brief Abstract linear system class.
37  *
38  * Linear system consists of Matrix, RHS and solution.
39  * It provides general methods for:
40  * - matrix preallocation
41  * - assembling matrix and RHS
42  * - application of linear constrains (Dirichlet conditions) either during assembly or
43  * on assembled system
44  * - solution of the system
45  * - output in Matlab format
46  *
47  * Method operates on the system as single object. But some methods for matrix only manipulation
48  * can be provided until we have matrix as separate class.
49  *
50  * TODO:
51  * - simplify constructors
52  * - introduce set_solution_vector() and set_rhs() in order to solve multiple systems with same matrix
53  * - simplify constructor, make common interface, rather introduce particular seters for particular solvers
54  * - which parameters expose through Input::Record
55  * - why we need lsize in constructors if we have Distribution
56  */
57 
58 #include <mpi.h> // for MPI_Comm, MPI_...
59 #include <stdlib.h> // for NULL, abs
60 #include <string.h> // for memcpy
61 //
62 
63 #include <cmath> // for abs, fabs
64 #include <new> // for operator new[]
65 #include <string> // for basic_string
66 #include <typeinfo> // for type_info
67 #include <utility> // for swap, pair
68 #include <vector> // for vector, allocator
69 
70 #include <armadillo>
71 
72 // PETSc includes
73 #include "petsclog.h" // for MPI_Allreduce
74 #include "petscmat.h" // for Mat, MatStructure
75 #include "petscmath.h" // for PetscScalar
76 #include "petscsys.h" // for PetscErrorCode
77 #include "petscvec.h" // for Vec, VecCreate...
78 
79 #include "input/accessors.hh" // for Record
80 #include "la/local_system.hh"
81 #include "la/distribution.hh" // for Distribution
82 #include "system/exceptions.hh" // for ExcAssertMsg::...
83 #include "system/system.hh" // for chkerr
84 #include "system/asserts.hh" // for ASSERT, msg
85 
86 namespace Input { namespace Type { class Abstract; } }
87 
88 
89 class LinSys
90 {
91 friend class SchurComplement;
92 public:
93  // Abstract Input Record for LinSys initialization
95 
96  typedef enum {
97  INSERT=INSERT_VALUES,
98  ADD=ADD_VALUES,
101  NONE
103 
104  struct SolveInfo {
105  SolveInfo(int cr, int nits)
106  : converged_reason(cr), n_iterations(nits)
107  {}
110  };
111 
112 protected:
113  typedef std::pair<unsigned,double> Constraint_;
115 
116 public:
117  /**
118  * Constructor.
119  * Constructor of abstract class should not be called directly, but is used for initialization of member common
120  * to all particular solvers.
121  *
122  * @param comm - MPI communicator
123  *
124  * TODO: Vector solution_ is now initialized to NULL, but it should be rather allocated
125  * in the constructor instead of the method set_solution().
126  */
127  LinSys(const Distribution *rows_ds)
128  : comm_( rows_ds->get_comm() ), status_( NONE ), lsize_( rows_ds->lsize() ), rows_ds_(rows_ds),
129  symmetric_( false ), positive_definite_( false ), negative_definite_( false ),
130  spd_via_symmetric_general_( false ), solution_(NULL), v_solution_(NULL),
131  own_vec_(false), own_solution_(false)
132  {
133  int lsizeInt = static_cast<int>( rows_ds->lsize() );
134  int sizeInt;
135  MPI_Allreduce ( &lsizeInt, &sizeInt, 1, MPI_INT, MPI_SUM, comm_ );
136  size_ = static_cast<unsigned>( sizeInt );
137 
142  };
143 
144  /**
145  * Copy constructor.
146  */
147  LinSys(LinSys &other)
148  : r_tol_(other.r_tol_), a_tol_(other.a_tol_), d_tol_(other.d_tol_), max_it_(other.max_it_), comm_(other.comm_), status_(other.status_),
149  lsize_( other.rows_ds_->lsize() ), size_(other.size_), rows_ds_(other.rows_ds_), symmetric_(other.symmetric_),
154 
155  {
156  ASSERT_PERMANENT( false ).error("Using copy constructor of LinSys is not allowed!");
157  set_solution(other.v_solution_);
158  };
159 
160  /**
161  * Returns global system size.
162  */
163  inline unsigned int size()
164  {
165  return size_;
166  }
167 
168  /**
169  * Returns local system size. (for partitioning of solution vectors)
170  * for PETSC_MPIAIJ it is also partitioning of the matrix
171  */
172  inline unsigned int vec_lsize()
173  {
174  return lsize_;
175  }
176 
177  /**
178  * Returns PETSC matrix (only for PETSC solvers)
179  *
180  * If matrix is changed, method set_matrix_changed() must be called.
181  * Example:
182  * @CODE
183  * MatDiagonalSet(schur->get_matrix(), new_diagonal, ADD_VALUES);
184  * schur->set_matrix_changed();
185  * @ENDCODE
186  */
187  virtual const Mat *get_matrix()
188  {
189  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function get_matrix is not implemented for linsys type\n.");
190  return NULL;
191  }
192 
193  /**
194  * Returns RHS vector (only for PETSC solvers)
195  *
196  * If vector is changed, method set_rhs_changed() must be called.
197  * Example:
198  * @CODE
199  * VecScale(schur->get_rhs(), -1.0);
200  * schur->set_rhs_changed();
201  * @ENDCODE
202  */
203  virtual const Vec *get_rhs()
204  {
205  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function get_rhs is not implemented for linsys type\n.");
206  return NULL;
207  }
208 
209  /**
210  * Sets matrix changed flag.
211  */
213  { matrix_changed_ = true;}
214 
215  /**
216  * Sets rhs changed flag (only for PETSC solvers)
217  */
219  { rhs_changed_ = true; }
220 
221  /**
222  * Set relative tolerance, absolute tolerance, and maximum number of iterations of the linear solver.
223  *
224  * For each of these three parameters we first look for the value at user input
225  * if not set we use the value provided to this method and finally the default values are
226  * set by the call of this method in the constructor.
227  */
228  virtual void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) = 0;
229 
230  /**
231  * Returns true if the system matrix has changed since the last solve.
232  */
234  { return matrix_changed_;}
235 
236  /**
237  * Returns true if the system RHS has changed since the last solve.
238  */
240  { return rhs_changed_;}
241 
242 
243  /**
244  * Sets PETSC matrix (only for PETSC solvers)
245  */
246  virtual PetscErrorCode set_matrix(Mat&, MatStructure)
247  {
248  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function set_matrix is not implemented for linsys type \n.");
249  return 0;
250  }
251 
252  /**
253  * Sets RHS vector (only for PETSC solvers)
254  */
255  virtual PetscErrorCode set_rhs(Vec&)
256  {
257  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function set_rhs is not implemented for linsys type \n.");
258  return 0;
259  }
260 
261  /**
262  * Clears entries of the matrix
263  */
264  virtual PetscErrorCode mat_zero_entries()
265  {
266  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function mat_zero_entries is not implemented for linsys type \n.");
267  return 0;
268  }
269 
270  /**
271  * Clears entries of the right-hand side
272  */
273  virtual PetscErrorCode rhs_zero_entries()
274  {
275  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function vec_zero_entries is not implemented for linsys type \n.");
276  return 0;
277  }
278 
279  /**
280  * Returns PETSC vector with solution. Underlying array can be provided on construction.
281  */
282  const Vec &get_solution()
283  {
284  return solution_;
285  }
286 
287  /**
288  * Set PETSc solution
289  */
290  void set_solution(Vec sol_vec) {
291  solution_ = sol_vec;
292  own_vec_ = false;
293  own_solution_ = false;
294  double *out_array;
295  VecGetArray( solution_, &out_array );
296  v_solution_ = out_array;
297  VecRestoreArray( solution_, &out_array );
298  }
299 
300  /**
301  * Create PETSc solution
302  */
303  void set_solution(double *sol_array) {
304  v_solution_ = sol_array;
305  own_vec_ = true;
306  own_solution_ = false;
307  PetscErrorCode ierr;
308  ierr = VecCreateMPIWithArray( comm_,1, rows_ds_->lsize(), PETSC_DECIDE, v_solution_, &solution_ ); CHKERRV( ierr );
309  }
310 
311  /**
312  * Create PETSc solution
313  */
314  void set_solution() {
315  v_solution_ = new double[ rows_ds_->lsize() + 1 ];
316  own_vec_ = true;
317  own_solution_ = true;
318  PetscErrorCode ierr;
319  ierr = VecCreateMPIWithArray( comm_,1, rows_ds_->lsize(), PETSC_DECIDE, v_solution_, &solution_ ); CHKERRV( ierr );
320  }
321 
322  /**
323  * Returns PETSC subarray with solution. Underlying array can be provided on construction.
324  */
326  {
327  return v_solution_;
328  }
329 
330  /**
331  * Switch linear system into allocating assembly. (only for PETSC_MPIAIJ_preallocate_by_assembly)
332  */
333  virtual void start_allocation()
334  {
335  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function start_allocation is not implemented for linsys type \n.");
336  }
337 
338  /**
339  * Switch linear system into adding assembly. (the only one supported by triplets ??)
340  */
341  virtual void start_add_assembly()
342  {
343  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function start_add_assembly is not implemented for linsys type \n.");
344  }
345 
346  /**
347  * Switch linear system into insert assembly. (not currently used)
348  */
349  virtual void start_insert_assembly()
350  {
351  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function start_insert_assembly is not implemented for linsys type \n.");
352  }
353 
354  /**
355  * Finish assembly of the whole system. For PETSC this should call MatEndAssembly with MAT_FINAL_ASSEMBLY
356  */
357  virtual void finish_assembly( )=0;
358 
359  /**
360  * Assembly full rectangular submatrix into the system matrix.
361  * Should be virtual, implemented differently in particular solvers.
362  */
363  virtual void mat_set_values(int nrow,int *rows,int ncol,int *cols,double *vals)=0;
364 
365  /**
366  * Shortcut for assembling just one element into the matrix.
367  * Similarly we can provide method accepting armadillo matrices.
368  */
369  void mat_set_value(int row,int col,double val)
370  { mat_set_values(1,&row,1,&col,&val); }
371 
372  /**
373  * Set values of the system right-hand side.
374  * Should be virtual, implemented differently in particular solvers.
375  */
376  virtual void rhs_set_values(int nrow,int *rows,double *vals)=0;
377 
378  /**
379  * Shorcut for assembling just one element into RHS vector.
380  */
381  void rhs_set_value(int row,double val)
382  { rhs_set_values(1,&row,&val); }
383 
384 
385  /// Set values in the system matrix and values in the right-hand side vector on corresponding rows.
386  inline void set_values(int nrow,int *rows,int ncol,int *cols,PetscScalar *mat_vals, PetscScalar *rhs_vals)
387  {
388  mat_set_values(nrow, rows, ncol, cols, mat_vals);
389  rhs_set_values(nrow, rows, rhs_vals);
390  }
391 
393  // local.eliminate_solution();
394  arma::mat tmp = local.matrix.t();
395 // DBGCOUT(<< "\n" << tmp);
396 // DBGCOUT(<< "row dofs: \n");
397 // for(unsigned int i=0; i< local.row_dofs.size(); i++)
398 // cout << local.row_dofs[i] << " ";
399 // cout <<endl;
400 
401  // This is always done only once, see implementation.
402  mat_set_values(local.matrix.n_rows, (int *)(local.row_dofs.memptr()),
403  local.matrix.n_cols, (int *)(local.col_dofs.memptr()),
404  tmp.memptr());
405 
406  rhs_set_values(local.matrix.n_rows, (int *)(local.row_dofs.memptr()),
407  local.rhs.memptr());
408  }
409 
410  /// Sets local system, considering that the dof indices are locallized on current processor.
411  /// @param local_to_global_map - maps the local dof indices to global ones
412  void set_local_system(LocalSystem & local, const std::vector<LongIdx> & local_to_global_map){
413  // transpose the matrix to provide proper column format
414  arma::mat tmp = local.matrix.t();
415 
416  // map local dof indices to global
417  uint m = local.row_dofs.n_elem,
418  n = local.col_dofs.n_elem;
419  int row_dofs[m];
420  int col_dofs[n];
421  for (uint i=0; i<m; i++)
422  row_dofs[i]= local_to_global_map[local.row_dofs[i]];
423  for (uint i=0; i<n; i++)
424  col_dofs[i]= local_to_global_map[local.col_dofs[i]];
425 
426  mat_set_values(m, row_dofs, n, col_dofs, tmp.memptr());
427  rhs_set_values(m, row_dofs, local.rhs.memptr());
428  }
429 
430  /**
431  * Add given dense matrix to preallocation.
432  */
433  //virtual void allocate(std::vector<int> rows, std::vector<int> cols);
434 
435  /**
436  * Shortcut to assembly into matrix and RHS in one call, possibly apply Dirichlet boundary conditions.
437  * @p row_dofs - are global indices of rows of dense @p matrix and rows of dense vector @rhs in global system
438  * @p col_dofs - are global indices of columns of the matrix, and possibly
439  *
440  * Application of Dirichlet conditions:
441  * 1) Rows with negative dofs are set to zero.
442  * 2) Cols with negative dofs are eliminated.
443  * 3) If there are entries on global diagonal. We determine value K either from diagonal of local matrix, or (if it is zero) from
444  * diagonal average.
445  *
446  * Caveats:
447  * - can not set dirichlet condition on zero dof
448  * - Armadillo stores matrix in column first form (Fortran like) which makes it not well suited
449  * for passing local matrices.
450  *
451  */
452  void set_values(std::vector<int> &row_dofs, std::vector<int> &col_dofs,
453  const arma::mat &matrix, const arma::vec &rhs,
454  const arma::vec &row_solution, const arma::vec &col_solution)
455 
456  {
457  arma::mat tmp = matrix.t();
458  arma::vec tmp_rhs = rhs;
459  bool negative_row = false;
460  bool negative_col = false;
461 
462  for(unsigned int l_row = 0; l_row < row_dofs.size(); l_row++)
463  if (row_dofs[l_row] < 0) {
464  tmp_rhs(l_row)=0.0;
465  tmp.col(l_row).zeros();
466  negative_row=true;
467  }
468 
469  for(unsigned int l_col = 0; l_col < col_dofs.size(); l_col++)
470  if (col_dofs[l_col] < 0) {
471  tmp_rhs -= matrix.col(l_col) * col_solution[l_col];
472  tmp.row(l_col).zeros();
473  negative_col=true;
474  }
475 
476 
477  if (negative_row && negative_col) {
478  // look for diagonal entry
479  for(unsigned int l_row = 0; l_row < row_dofs.size(); l_row++)
480  if (row_dofs[l_row] < 0)
481  for(unsigned int l_col = 0; l_col < col_dofs.size(); l_col++)
482  if (col_dofs[l_col] < 0 && row_dofs[l_row] == col_dofs[l_col]) {
483  double new_diagonal = fabs(matrix.at(l_row,l_col));
484  if (new_diagonal == 0.0) {
485  if (matrix.is_square()) {
486  new_diagonal = arma::sum( abs(matrix.diag())) / matrix.n_rows;
487  } else {
488  new_diagonal = arma::accu( abs(matrix) ) / matrix.n_elem;
489  }
490  }
491  tmp.at(l_col, l_row) = new_diagonal;
492  tmp_rhs(l_row) = new_diagonal * row_solution[l_row];
493  }
494 
495  }
496 
497  if (negative_row)
498  for(int &row : row_dofs) row=abs(row);
499 
500  if (negative_col)
501  for(int &col : col_dofs) col=abs(col);
502 
503 
504  mat_set_values(row_dofs.size(), const_cast<int *>(&(row_dofs[0])),
505  col_dofs.size(), const_cast<int *>(&(col_dofs[0])), tmp.memptr() );
506  rhs_set_values(row_dofs.size(), const_cast<int *>(&(row_dofs[0])), tmp_rhs.memptr() );
507  }
508 
509  /**
510  * Adds Dirichlet constrain.
511  * @param row - global number of row that should be eliminated.
512  * @param value - solution value at the given row
513  */
514  void add_constraint(int row, double value) {
515 
516  constraints_.push_back( Constraint_( static_cast<unsigned>( row ), value ) );
517  }
518 
519  /**
520  * Apply constrains to assembled matrix. Constrains are given by pairs: global row index, value.
521  * i.e. typedef pair<unsigned int, double> Constrain;
522  *
523  * What is th meaning of ( const double factor ) form Cambridge code?
524  */
525  virtual void apply_constrains( double scalar )=0;
526 
527  /**
528  * Solve the system and return convergence reason.
529  */
530  virtual SolveInfo solve()=0;
531 
532  /**
533  * Returns norm of the initial right hand side
534  */
536  return residual_norm_;
537  };
538 
539  /**
540  * Explicitly compute residual and its norm for current solution.
541  */
542  virtual double compute_residual() =0;
543 
544  /**
545  * Returns information on relative solver accuracy
546  */
548  return r_tol_;
549  };
550 
551  /**
552  * Returns information on absolute solver accuracy
553  */
554  virtual double get_absolute_accuracy(){
555  return 0.0;
556  };
557 
558  /**
559  * Provides user knowledge about symmetry.
560  */
561  inline void set_symmetric(bool flag = true)
562  {
563  symmetric_ = flag;
564  if (!flag) {
565  set_positive_definite(false);
566  set_negative_definite(false);
567  }
568  }
569 
570  inline bool is_symmetric()
571  { return symmetric_; }
572 
573  /**
574  * Provides user knowledge about positive definiteness.
575  */
576  inline void set_positive_definite(bool flag = true)
577  {
578  positive_definite_ = flag;
579  if (flag) {
580  set_symmetric();
581  set_negative_definite(false);
582  }
583  }
584 
585  /**
586  * Provides user knowledge about negative definiteness.
587  */
588  inline void set_negative_definite(bool flag = true)
589  {
590  negative_definite_ = flag;
591  if (flag) {
592  set_symmetric();
593  set_positive_definite(false);
594  }
595  }
596 
597  inline bool is_positive_definite()
598  { return positive_definite_; }
599 
600  inline bool is_negative_definite()
601  { return negative_definite_; }
602 
603  /// TODO: In fact we want to know if the matrix is already preallocated
604  /// However to do this we need explicit finalisation of preallocating cycle.
605  inline bool is_new() {
606  return ( status_ == NONE );
607  }
608 
609  inline bool is_preallocated() {
610  return ( status_ == INSERT || status_ == ADD);
611  }
612 
613  /**
614  * Provides user knowledge about positive definiteness via symmetric general approach.
615  * This is useful for solving Darcy flow by mixed hybrid method, where blocks on subdomains are saddle point but
616  * interface among subdomains is only at the block of Lagrange multipliers and is symmetric positive definite.
617  * Problem condensed to interface can thus be solved by PCG method, although original problem is saddle point.
618  */
619  inline void set_spd_via_symmetric_general(bool flag = true)
620  {
622  if (flag) set_symmetric();
623  }
624 
626  { return spd_via_symmetric_general_; }
627 
628 
629  /**
630  * Output the system in the Matlab format possibly with given ordering.
631  * Rather we shoud provide output operator <<, since it is more flexible.
632  */
633  virtual void view(string)
634  {
635  ASSERT_PERMANENT( false )(typeid(*this).name()).error("Function view is not implemented for linsys type %s \n.");
636  }
637 
638  /**
639  * Sets basic parameters of LinSys defined by user in input file and used to calculate
640  */
641  virtual void set_from_input(const Input::Record in_rec)
642  {
643  in_rec_ = in_rec;
645  }
646 
647  /**
648  * Get precision of solving
649  */
650  virtual double get_solution_precision() = 0;
651 
652  virtual ~LinSys()
653  {
654  if ( own_vec_ && solution_ ) { chkerr(VecDestroy(&solution_)); }
655  if ( own_solution_ ) delete[] v_solution_;
656  }
657 
658 protected:
659  // Default values initialized in constructor
660  static constexpr double default_r_tol_ = 1e-7;
661  static constexpr double default_a_tol_ = 1e-11;
662  static constexpr double default_d_tol_ = 10000;
663  static constexpr unsigned int default_max_it_ = 1000;
664 
665  double r_tol_; ///< relative tolerance of linear solver
666  double a_tol_; ///< absolute tolerance of linear solver
667  double d_tol_; ///< tolerance for divergence of linear solver
668  unsigned int max_it_; ///< maximum number of iterations of linear solver
669 
671  SetValuesMode status_; //!< Set value status of the linear system.
672 
673  const unsigned lsize_; //!< local number of matrix rows (non-overlapping division of rows)
674  unsigned size_; //!< global number of matrix rows, i.e. problem size
675 
676  const Distribution * rows_ds_; //!< final distribution of rows of MH matrix
677 
682 
683  bool matrix_changed_; //!< true if the matrix was changed since the last solve
684  bool rhs_changed_; //!< true if the right hand side was changed since the last solve
685 
686  Vec solution_; //!< PETSc vector constructed with vb array.
687  double *v_solution_; //!< local solution array pointing into Vec solution_
688  bool own_vec_; //!< Indicates if the solution vector has been allocated by this class
689  bool own_solution_; //!< Indicates if the solution array has been allocated by this class
690 
691  double residual_norm_; //!< local solution array pointing into Vec solution_
692 
694 
695  std::vector<double> globalSolution_; //!< global solution in numbering for linear system
696 
697  Input::Record in_rec_; // structure contained parameters of LinSys defined in input file
698 
699 };
700 
701 #endif /* LA_LINSYS_HH_ */
Definitions of ASSERTS.
#define ASSERT_PERMANENT(expr)
Allow use shorter versions of macro names if these names is not used with external library.
Definition: asserts.hh:348
unsigned int lsize(int proc) const
get local size
Accessor to the data with type Type::Record.
Definition: accessors.hh:291
Class for declaration of polymorphic Record.
const Vec & get_solution()
Definition: linsys.hh:282
static constexpr double default_r_tol_
Definition: linsys.hh:660
unsigned int max_it_
maximum number of iterations of linear solver
Definition: linsys.hh:668
virtual const Vec * get_rhs()
Definition: linsys.hh:203
void set_negative_definite(bool flag=true)
Definition: linsys.hh:588
void set_local_system(LocalSystem &local, const std::vector< LongIdx > &local_to_global_map)
Definition: linsys.hh:412
void set_solution(Vec sol_vec)
Definition: linsys.hh:290
void set_rhs_changed()
Definition: linsys.hh:218
void set_values(std::vector< int > &row_dofs, std::vector< int > &col_dofs, const arma::mat &matrix, const arma::vec &rhs, const arma::vec &row_solution, const arma::vec &col_solution)
Definition: linsys.hh:452
static constexpr double default_a_tol_
Definition: linsys.hh:661
double get_residual_norm()
Definition: linsys.hh:535
bool spd_via_symmetric_general_
Definition: linsys.hh:681
double get_relative_accuracy()
Definition: linsys.hh:547
unsigned int size()
Definition: linsys.hh:163
double d_tol_
tolerance for divergence of linear solver
Definition: linsys.hh:667
virtual ~LinSys()
Definition: linsys.hh:652
void set_matrix_changed()
Definition: linsys.hh:212
MPI_Comm comm_
Definition: linsys.hh:670
std::vector< double > globalSolution_
global solution in numbering for linear system
Definition: linsys.hh:695
bool is_spd_via_symmetric_general()
Definition: linsys.hh:625
void set_solution()
Definition: linsys.hh:314
double * get_solution_array()
Definition: linsys.hh:325
bool is_negative_definite()
Definition: linsys.hh:600
virtual void set_from_input(const Input::Record in_rec)
Definition: linsys.hh:641
virtual void rhs_set_values(int nrow, int *rows, double *vals)=0
bool is_positive_definite()
Definition: linsys.hh:597
bool symmetric_
Definition: linsys.hh:678
Input::Record in_rec_
Definition: linsys.hh:697
virtual double get_solution_precision()=0
virtual void start_add_assembly()
Definition: linsys.hh:341
virtual void finish_assembly()=0
static constexpr double default_d_tol_
Definition: linsys.hh:662
bool is_new()
Definition: linsys.hh:605
bool is_matrix_changed()
Definition: linsys.hh:233
void rhs_set_value(int row, double val)
Definition: linsys.hh:381
virtual const Mat * get_matrix()
Definition: linsys.hh:187
bool is_rhs_changed()
Definition: linsys.hh:239
virtual void apply_constrains(double scalar)=0
void set_solution(double *sol_array)
Definition: linsys.hh:303
virtual void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it)=0
ConstraintVec_ constraints_
Definition: linsys.hh:693
unsigned size_
global number of matrix rows, i.e. problem size
Definition: linsys.hh:674
unsigned int vec_lsize()
Definition: linsys.hh:172
virtual void mat_set_values(int nrow, int *rows, int ncol, int *cols, double *vals)=0
static constexpr unsigned int default_max_it_
Definition: linsys.hh:663
bool positive_definite_
Definition: linsys.hh:679
virtual double compute_residual()=0
bool is_preallocated()
Definition: linsys.hh:609
bool own_solution_
Indicates if the solution array has been allocated by this class.
Definition: linsys.hh:689
void set_local_system(LocalSystem &local)
Definition: linsys.hh:392
void set_symmetric(bool flag=true)
Definition: linsys.hh:561
virtual PetscErrorCode rhs_zero_entries()
Definition: linsys.hh:273
double a_tol_
absolute tolerance of linear solver
Definition: linsys.hh:666
void mat_set_value(int row, int col, double val)
Definition: linsys.hh:369
SetValuesMode
Definition: linsys.hh:96
@ NONE
Definition: linsys.hh:101
@ INSERT
Definition: linsys.hh:97
@ DONE
Definition: linsys.hh:100
@ ADD
Definition: linsys.hh:98
@ ALLOCATE
Definition: linsys.hh:99
void add_constraint(int row, double value)
Definition: linsys.hh:514
std::pair< unsigned, double > Constraint_
Definition: linsys.hh:113
SetValuesMode status_
Set value status of the linear system.
Definition: linsys.hh:671
bool own_vec_
Indicates if the solution vector has been allocated by this class.
Definition: linsys.hh:688
virtual void start_insert_assembly()
Definition: linsys.hh:349
virtual void start_allocation()
Definition: linsys.hh:333
virtual PetscErrorCode set_matrix(Mat &, MatStructure)
Definition: linsys.hh:246
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
virtual SolveInfo solve()=0
virtual void view(string)
Definition: linsys.hh:633
void set_values(int nrow, int *rows, int ncol, int *cols, PetscScalar *mat_vals, PetscScalar *rhs_vals)
Set values in the system matrix and values in the right-hand side vector on corresponding rows.
Definition: linsys.hh:386
double * v_solution_
local solution array pointing into Vec solution_
Definition: linsys.hh:687
double residual_norm_
local solution array pointing into Vec solution_
Definition: linsys.hh:691
virtual PetscErrorCode set_rhs(Vec &)
Definition: linsys.hh:255
void set_spd_via_symmetric_general(bool flag=true)
Definition: linsys.hh:619
bool is_symmetric()
Definition: linsys.hh:570
LinSys(const Distribution *rows_ds)
Definition: linsys.hh:127
void set_positive_definite(bool flag=true)
Definition: linsys.hh:576
double r_tol_
relative tolerance of linear solver
Definition: linsys.hh:665
virtual double get_absolute_accuracy()
Definition: linsys.hh:554
bool matrix_changed_
true if the matrix was changed since the last solve
Definition: linsys.hh:683
const unsigned lsize_
local number of matrix rows (non-overlapping division of rows)
Definition: linsys.hh:673
virtual PetscErrorCode mat_zero_entries()
Definition: linsys.hh:264
bool negative_definite_
Definition: linsys.hh:680
std::vector< Constraint_ > ConstraintVec_
Definition: linsys.hh:114
LinSys(LinSys &other)
Definition: linsys.hh:147
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
LocDofVec col_dofs
Definition: local_system.hh:54
arma::vec rhs
local system RHS
LocDofVec row_dofs
Definition: local_system.hh:54
arma::mat matrix
local system matrix
Support classes for parallel programing.
static constexpr bool value
Definition: json.hpp:87
unsigned int uint
#define MPI_INT
Definition: mpi.h:160
#define MPI_Allreduce(sendbuf, recvbuf, count, datatype, op, comm)
Definition: mpi.h:612
int MPI_Comm
Definition: mpi.h:141
#define MPI_SUM
Definition: mpi.h:196
ArmaMat< double, N, M > mat
Definition: armor.hh:888
ArmaVec< double, N > vec
Definition: armor.hh:885
Abstract linear system class.
Definition: balance.hh:40
int converged_reason
Definition: linsys.hh:108
SolveInfo(int cr, int nits)
Definition: linsys.hh:105
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.
Definition: system.hh:142