Flow123d  master-f44eb46
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,
70  NONE
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;
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 
LinSys Vector RHS
LinSys * block(i, j) MatrixArray *matrix() virtual solve(Vector solution
virtual SolveInfo solve()=0
std::vector< std::vector< MatrixSimple * > > mat_array
MatrixArray vice konstruktoru pro ruzne slepovani MatrixSimple a jiz existujicich MatrixArray MatrixSimple * block(i_row, i_col) multiply_vec(Vector x
virtual void preallocate_values(int nrow, int *rows, int ncol, int *cols)=0
Distribution col_ds
Distribution of multiplied vector.
void mat_set_value(int row, int col, PetscScalar val)
Set one element of the system matrix.
virtual void start_allocation()=0
unsigned int size()
Get global system size.
SetValuesMode
possible states of the matrix
void set_symmetric(bool flag=true)
void multiply(Vec x, Vec y)
multiply PETSC vector: y=Ax
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.
void finalize(MatAssemblyType assembly_type=MAT_FINAL_ASSEMBLY)
bool positive_definite
Flag for positive definite system.
bool is_symmetric()
void multiply_add(Vec x, Vec b, Vec y)
multiply vector and ad onother one: y=Ax+b
Distribution & get_row_ds()
void mat_set_values(int nrow, int *rows, int ncol, int *cols, PetscScalar *vals)
Set full rectangular submatrix of the system matrix.
const Mat & get_matrix()
Get matrix. SHOULD NOT BE USED !!!
virtual void view_local_matrix()=0
SetValuesMode status
Set value status of the linear system.
void multiply_mat(...)
multiply matrix .. has to deal with preallocation of result and reuse of structure
virtual void preallocate_matrix()=0
void aypx(...)
bool is_positive_definite()
virtual ~LinSys()
void set_positive_definite(bool flag=true)
MatrixSimple(unsigned int loc_row_size, unsigned int col_size)
Construct a parallel system with given local size.
void axpy(...)
void start_add_assembly()
Mat matrix
Petsc matrix of the problem.
Distribution & get_col_ds()
Distribution row_ds
Distribution of continuous blocks of system rows among the processors.
bool symmetric
Flag for the symmetric system.
void start_insert_assembly()
class MatrixArray VecArray