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