Flow123d  JS_before_hm-1804-gf2ad740aa
la_linsys_new.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 la_linsys_new.hh
15  * @brief
16  */
17 
18 /**
19  * Classes for block linear systems. Namely it should provide following functionality:
20  * 1) interface to various PETSC (or possibly other) matrices, with preallocation
21  * through the set_values method
22  * 2) For coupling problems creation of the global linear system by connectiong sub systems of individual
23  * equations. Working cycle could be:
24  * 1. creation of sub problems up to the global problem (just matrix dimensions and types)
25  * 2. pre assembly - > preallocation
26  * 3. assembly - by the same functions
27  * 4. matrix modifications (schur complements, assembly of one global matrix, preconditioners ...
28  * 5. iterative solution in parallel
29  * 6. for nonlinear/time problems: jump to 3. and reuse the structures
30  *
31  * 3) Solution through explicit or implicit schur complements
32  * 4) Preconditioning for coupled problems .
33  * 5) In assembly, the access to the matrices in LinSys hierarchy is through successive call of block(i,j) method, but
34  * there also should be a class or another mechanism how to represent a coordinate in the hierarchy. Then one can access a MatriSimple
35  * directly by two coordinates and constantly access correct part of rhs vector
36  *
37  * Virtual test cases:
38  * - Only water flow solved by 1, 2, 3, or 4 (Newton BC) Schur complements, explict and implicit
39  * Ideal: switch between Implicit and explicit only changing the particular class or some thing similarly simple
40  * -
41  *
42  * Unsolved design problems:
43  * - PETSC neumoznuje blokove vektory, tj. vektory by mely byt globalni a z nich by se mely pouzivat podvektory
44  * mozne reseni: blokovy globalni vektor se hierarchicky vytvari (pouze velikosti) a pak slouzi jako
45  * "fabrika" na globalni vektory i subvektory a nebo pro pristup k subvektorum globalniho vektoru
46  * .. detaily
47  *
48  */
49 
50 
51 
52 /**
53  * Abstract class for frontend for PETSc matrices.
54  * purpose:
55  * - allow preallocation through double assembly
56  * - make interface to PETSC since PETSC evolve and it is simpler to make changes on one place
57  * - make documented interface to PETSC and document used PETSC functionality this way
58  *
59  * Questions:
60  * - how to setup column distribution. It is necessary only for some of PETSC formats, isn't it?
61  */
63 {public:
64  /// possible states of the matrix
65  typedef enum {
66  INSERT=INSERT_VALUES,
67  ADD=ADD_VALUES,
71  } SetValuesMode;
72 
73  /// Construct a parallel system with given local size.
74  MatrixSimple(unsigned int loc_row_size, unsigned int col_size);
75 
76  /// multiply PETSC vector: y=Ax
77  void multiply(Vec x, Vec y);
78  /// multiply vector and ad onother one: y=Ax+b
79  void multiply_add(Vec x, Vec b, Vec y);
80  /// multiply matrix .. has to deal with preallocation of result and reuse of structure
81  void multiply_mat(...);
82  /// add and scale sparse matrices ... has to deal with different nonzero structures,
83  /// in such a case petsc rutines are inceribly slow
84  void axpy(...);
85  void aypx(...);
86 
87  /// @name access members @{
88  /// Access to distribution of rows, thhrou that you can get also matix size and local size. But for convenience
89  /// we provide also direct functions
92  /// Get global system size.
93  inline unsigned int size();
94  /// Get matrix. SHOULD NOT BE USED !!!
95  inline const Mat &get_matrix()
96  { return matrix; }
97  /// @}
98 
99  virtual void start_allocation()=0;
100  void start_add_assembly();
101  void start_insert_assembly();
102  void finalize(MatAssemblyType assembly_type=MAT_FINAL_ASSEMBLY);
103 
104  virtual void preallocate_matrix()=0;
105  virtual void preallocate_values(int nrow,int *rows,int ncol,int *cols)=0;
106 
107  virtual void view_local_matrix()=0;
108 
109  /// Set full rectangular submatrix of the system matrix.
110  void mat_set_values(int nrow,int *rows,int ncol,int *cols,PetscScalar *vals);
111 
112  /// Set one element of the system matrix.
113  inline void mat_set_value(int row,int col,PetscScalar val);
114 
115  inline void set_symmetric(bool flag = true);
116  inline bool is_symmetric();
117 
118  inline void set_positive_definite(bool flag = true);
119  inline bool is_positive_definite();
120 
121  /// Output the matrix in the Matlab format possibly with given renumbering of rows/cols
122  void view(std::ostream output_stream, int * row_mapping = NULL, int * col_mapping = NULL);
123 
124  virtual ~LinSys();
125 
126 protected:
127  Distribution row_ds; ///< Distribution of continuous blocks of system rows among the processors.
128  Distribution col_ds; ///< Distribution of multiplied vector.
129  bool symmetric; ///< Flag for the symmetric system.
130  bool positive_definite; ///< Flag for positive definite system.
131  SetValuesMode status; ///< Set value status of the linear system.
132 
133  Mat matrix; ///< Petsc matrix of the problem.
134 
135 };
136 
137 /**
138  * Array of MatrixSimple
139  *
140  */
142 {
143 private:
145 
146 public:
147  MatrixArray ... vice konstruktoru pro ruzne slepovani MatrixSimple a jiz existujicich MatrixArray
148 
149  MatrixSimple *block(i_row, i_col)
150  multiply_vec( Vector x, Vector y)
151 }
152 
153 /**
154  * Array or hierarchy over the global vector.
155  */
156 class VecArray
157 {
158 
159 };
160 
161 
162 /**
163  * LinSys - matrix and a particular way how to compute solution for given RHS, i.e. this class can perform
164  * action of the matrix inverse
165  *
166  * This is abstract class for members of possible hierarchical tree of the whole system.
167  *
168  */
169 class LinSys {
170  LinSys *block(i,j)
171  MatrixArray *matrix()
172 
173  virtual solve( Vector solution, Vector RHS);
174 }
175 
176 /**
177  * The LinSys for direct factorization, its matrix has to be simple. Factorization will be made through PETSc
178  * so the used backend library library can be influenced by on-line parameters.
179  */
180 class LinSysDirect : public LinSys
181 {
183 }
184 
185 /**
186  * Array of four submatrices A00 A10 A01 A11 where A00 has to be LinSys and A11 has to be Simple.
187  * Solves the system by recursive solution of schur complement system. Can be implemented as
188  * interface to PETSC functionality see PCFIELDSPLIT
189  */
190 class LinSysImplicitSchur : public LinSys
191 {
192  LinSys *block(i,j)
193  MatrixArray *matrix()
194 }
195 /**
196  * This LinSys should after assembly construct an explicit Schur complemnt system. User has to provide
197  * an LinSysDirect for one of the diagonal blocks.
198  *
199  * There is a problem how to make possible several successive schur complements.
200  * This can not be organised bottom up like ImplicitSchur.
201  */
202 class LinSysExplicitSchur : public LinSys
203 {
204 }
205 
206 /**
207  * Array of four matrices 00 01 10 11, where 00 and 11 are LinSys and their application is used as preconditioner for the whole system.
208  * again see PCFIELDSPLIT
209  */
210 class LinSysJacobiCoupling : public LinSys
211 {
212  LinSys *block11()
213  MatrixArray *block12()
214  MatrixArray *block21()
215  LinSys *block22()
216 }
217 
218 
219 
220 
221 
MatrixSimple::positive_definite
bool positive_definite
Flag for positive definite system.
Definition: la_linsys_new.hh:130
MatrixSimple::matrix
Mat matrix
Petsc matrix of the problem.
Definition: la_linsys_new.hh:133
MatrixSimple::col_ds
Distribution col_ds
Distribution of multiplied vector.
Definition: la_linsys_new.hh:128
MatrixSimple::MatrixSimple
MatrixSimple(unsigned int loc_row_size, unsigned int col_size)
Construct a parallel system with given local size.
MatrixSimple::finalize
void finalize(MatAssemblyType assembly_type=MAT_FINAL_ASSEMBLY)
LinSys
Definition: la_linsys_new.hh:169
MatrixSimple
Definition: la_linsys_new.hh:62
MatrixArray
Definition: la_linsys_new.hh:141
LinSys::RHS
LinSys Vector RHS
Definition: la_linsys_new.hh:173
std::vector
Definition: doxy_dummy_defs.hh:7
MatrixSimple::~LinSys
virtual ~LinSys()
MatrixSimple::view
void view(std::ostream output_stream, int *row_mapping=NULL, int *col_mapping=NULL)
Output the matrix in the Matlab format possibly with given renumbering of rows/cols.
LinSys::solve
virtual SolveInfo solve()=0
MatrixSimple::start_allocation
virtual void start_allocation()=0
MatrixSimple::start_insert_assembly
void start_insert_assembly()
MatrixSimple::is_symmetric
bool is_symmetric()
MatrixSimple::symmetric
bool symmetric
Flag for the symmetric system.
Definition: la_linsys_new.hh:129
MatrixSimple::INSERT
@ INSERT
Definition: la_linsys_new.hh:66
VecArray
class MatrixArray VecArray
Distribution
Definition: distribution.hh:50
MatrixSimple::ALLOCATE
@ ALLOCATE
Definition: la_linsys_new.hh:68
MatrixSimple::start_add_assembly
void start_add_assembly()
MatrixSimple::row_ds
Distribution row_ds
Distribution of continuous blocks of system rows among the processors.
Definition: la_linsys_new.hh:127
MatrixSimple::SetValuesMode
SetValuesMode
possible states of the matrix
Definition: la_linsys_new.hh:65
MatrixSimple::multiply_add
void multiply_add(Vec x, Vec b, Vec y)
multiply vector and ad onother one: y=Ax+b
MatrixSimple::preallocate_matrix
virtual void preallocate_matrix()=0
MatrixArray::block
MatrixArray vice konstruktoru pro ruzne slepovani MatrixSimple a jiz existujicich MatrixArray MatrixSimple * block(i_row, i_col) multiply_vec(Vector x
MatrixArray::mat_array
std::vector< std::vector< MatrixSimple * > > mat_array
Definition: la_linsys_new.hh:144
MatrixSimple::NONE
@ NONE
Definition: la_linsys_new.hh:70
MatrixSimple::view_local_matrix
virtual void view_local_matrix()=0
MatrixSimple::aypx
void aypx(...)
MatrixSimple::multiply_mat
void multiply_mat(...)
multiply matrix .. has to deal with preallocation of result and reuse of structure
MatrixSimple::get_row_ds
Distribution & get_row_ds()
MatrixSimple::preallocate_values
virtual void preallocate_values(int nrow, int *rows, int ncol, int *cols)=0
MatrixSimple::ADD
@ ADD
Definition: la_linsys_new.hh:67
MatrixSimple::get_col_ds
Distribution & get_col_ds()
MatrixSimple::DONE
@ DONE
Definition: la_linsys_new.hh:69
MatrixSimple::size
unsigned int size()
Get global system size.
MatrixSimple::set_positive_definite
void set_positive_definite(bool flag=true)
MatrixSimple::multiply
void multiply(Vec x, Vec y)
multiply PETSC vector: y=Ax
MatrixSimple::set_symmetric
void set_symmetric(bool flag=true)
MatrixSimple::status
SetValuesMode status
Set value status of the linear system.
Definition: la_linsys_new.hh:131
MatrixSimple::axpy
void axpy(...)
MatrixSimple::mat_set_values
void mat_set_values(int nrow, int *rows, int ncol, int *cols, PetscScalar *vals)
Set full rectangular submatrix of the system matrix.
MatrixSimple::get_matrix
const Mat & get_matrix()
Get matrix. SHOULD NOT BE USED !!!
Definition: la_linsys_new.hh:95
MatrixSimple::mat_set_value
void mat_set_value(int row, int col, PetscScalar val)
Set one element of the system matrix.
MatrixSimple::is_positive_definite
bool is_positive_definite()
LinSys::block
LinSys * block(i, j) MatrixArray *matrix() virtual solve(Vector solution