/[escript]/trunk/paso/src/SystemMatrix.h
ViewVC logotype

Annotation of /trunk/paso/src/SystemMatrix.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6939 - (hide annotations)
Mon Jan 20 03:37:18 2020 UTC (5 weeks, 3 days ago) by uqaeller
File MIME type: text/plain
File size: 11434 byte(s)
Updated the copyright header.


1 ksteube 1312
2 jfenwick 3981 /*****************************************************************************
3 ksteube 1811 *
4 uqaeller 6939 * Copyright (c) 2003-2020 by The University of Queensland
5 jfenwick 3981 * http://www.uq.edu.au
6 ksteube 1811 *
7     * Primary Business: Queensland, Australia
8 jfenwick 6112 * Licensed under the Apache License, version 2.0
9     * http://www.apache.org/licenses/LICENSE-2.0
10 ksteube 1811 *
11 jfenwick 3981 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 * Development 2012-2013 by School of Earth Sciences
13 uqaeller 6939 * Development from 2014-2017 by Centre for Geoscience Computing (GeoComp)
14     * Development from 2019 by School of Earth and Environmental Sciences
15     **
16 jfenwick 3981 *****************************************************************************/
17 dhawcroft 631
18 ksteube 1811
19 caltinay 4814 /****************************************************************************/
20 jgs 150
21 caltinay 4836 /* Paso: SystemMatrix */
22 jgs 150
23 caltinay 4814 /****************************************************************************/
24 jgs 150
25 gross 412 /* Copyrights by ACcESS Australia 2003,2004,2005,2006 */
26 jfenwick 2608 /* Author: Lutz Gross, l.gross@uq.edu.au */
27 jgs 150
28 caltinay 4814 /****************************************************************************/
29 jgs 150
30 caltinay 4836 #ifndef __PASO_SYSTEMMATRIX_H__
31     #define __PASO_SYSTEMMATRIX_H__
32 jgs 150
33 ksteube 1312 #include "SparseMatrix.h"
34 jgs 150 #include "SystemMatrixPattern.h"
35    
36 caltinay 5929 #include <escript/AbstractSystemMatrix.h>
37    
38 caltinay 4836 namespace paso {
39 gross 1552
40 aellery 6726 struct Options;
41 caltinay 5929 class SystemMatrix;
42 caltinay 4836 typedef boost::shared_ptr<SystemMatrix> SystemMatrix_ptr;
43     typedef boost::shared_ptr<const SystemMatrix> const_SystemMatrix_ptr;
44 jgs 150
45 caltinay 4836 typedef int SystemMatrixType;
46    
47 caltinay 5929 /// this class holds a (distributed) stiffness matrix
48 caltinay 6119 class SystemMatrix : public escript::AbstractSystemMatrix
49 caltinay 4814 {
50 caltinay 5929 public:
51     /// default constructor - throws exception.
52     SystemMatrix();
53 ksteube 1312
54 caltinay 5929 SystemMatrix(SystemMatrixType type, SystemMatrixPattern_ptr pattern,
55     dim_t rowBlockSize, dim_t columnBlockSize,
56     bool patternIsUnrolled, const escript::FunctionSpace& rowFS,
57     const escript::FunctionSpace& colFS);
58    
59 caltinay 4836 ~SystemMatrix();
60 jgs 150
61 caltinay 4836 /// Nullifies rows and columns in the matrix.
62     /// The rows and columns are marked by positive values in mask_row and
63     /// mask_col. Values on the main diagonal which are marked to set to
64     /// zero by both mask_row and mask_col are set to main_diagonal_value.
65 caltinay 5929 virtual void nullifyRowsAndCols(escript::Data& mask_row,
66     escript::Data& mask_col,
67     double main_diagonal_value);
68 jgs 150
69 caltinay 5929 virtual inline void saveMM(const std::string& filename) const
70     {
71     if (mpi_info->size > 1) {
72 caltinay 6200 //throw PasoException("SystemMatrix::saveMM: Only single rank supported.");
73     SparseMatrix_ptr merged(mergeSystemMatrix());
74     if (mpi_info->rank == 0)
75     merged->saveMM(filename.c_str());
76 caltinay 5929 } else {
77     mainBlock->saveMM(filename.c_str());
78     }
79     }
80    
81     virtual inline void saveHB(const std::string& filename) const
82     {
83     if (mpi_info->size > 1) {
84 caltinay 5996 throw PasoException("SystemMatrix::saveHB: Only single rank supported.");
85 caltinay 5929 } else if (!(type & MATRIX_FORMAT_CSC)) {
86 caltinay 5996 throw PasoException("SystemMatrix::saveHB: Only CSC format supported.");
87 caltinay 5929 } else {
88     mainBlock->saveHB_CSC(filename.c_str());
89     }
90     }
91    
92 caltinay 6393 virtual void resetValues(bool preserveSolverData = false);
93 caltinay 5929
94 caltinay 4836 /// Nullifies rows in the matrix.
95     /// The rows are marked by positive values in mask_row. Values on the
96     /// main diagonal which are marked to set to zero by mask_row are set
97     /// to main_diagonal_value.
98     void nullifyRows(double* mask_row, double main_diagonal_value);
99 jgs 150
100 caltinay 4836 void add(dim_t, index_t*, dim_t, dim_t, index_t*, dim_t, double*);
101 jgs 150
102 caltinay 4873 void makeZeroRowSums(double* left_over);
103 gross 3094
104 caltinay 4873 /// copies the col_coupleBlock into row_coupleBlock.
105 caltinay 4836 /// WARNING: this method uses mpi_requests of the coupler attached to the
106     /// matrix. No reordering on the received columns is performed.
107     /// In practice this means that components in
108     /// row_coupleBlock->pattern->index and
109     /// row_coupler->connector->recv->shared
110     /// are ordered by increasing value.
111     /// Note that send and receive row_coupler->connectors are swapping roles.
112     void copyColCoupleBlock();
113 gross 3369
114 caltinay 4836 void copyRemoteCoupleBlock(bool recreatePattern);
115 gross 3369
116 caltinay 4836 void fillWithGlobalCoordinates(double f1);
117 lgao 3827
118 caltinay 4836 void print() const;
119 jgs 150
120 caltinay 4849 /// Merges the system matrix which is distributed on several MPI ranks
121     /// into a complete sparse matrix on rank 0. Used by the Merged Solver.
122     SparseMatrix_ptr mergeSystemMatrix() const;
123 ksteube 1312
124 caltinay 4849 void mergeMainAndCouple(index_t** p_ptr, index_t** p_idx, double** p_val) const;
125 jgs 150
126 caltinay 4849 void mergeMainAndCouple_CSR_OFFSET0(index_t** p_ptr, index_t** p_idx, double** p_val) const;
127     void mergeMainAndCouple_CSR_OFFSET0_Block(index_t** p_ptr, index_t** p_idx, double** p_val) const;
128 caltinay 3594
129 caltinay 4849 void mergeMainAndCouple_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val) const;
130    
131 caltinay 4836 void copyMain_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val);
132 caltinay 3594
133 caltinay 4836 void extendedRowsForST(dim_t* degree_ST, index_t* offset_ST, index_t* ST);
134 jgs 150
135 caltinay 4836 void applyBalanceInPlace(double* x, bool RHS) const;
136 caltinay 3594
137 caltinay 4836 void applyBalance(double* x_out, const double* x, bool RHS) const;
138 caltinay 3594
139 caltinay 4836 void balance();
140 caltinay 3594
141 caltinay 4836 double getGlobalSize() const;
142 caltinay 3594
143 caltinay 4846 void setPreconditioner(Options* options);
144 caltinay 3594
145 caltinay 4836 /// Applies the preconditioner.
146 caltinay 4873 /// This method needs to be called within a parallel region.
147 caltinay 4836 /// Barrier synchronization is performed before the evaluation to make
148     /// sure that the input vector is available
149     void solvePreconditioner(double* x, double* b);
150 caltinay 3594
151 caltinay 4836 void freePreconditioner();
152 gross 3369
153 caltinay 4836 index_t* borrowMainDiagonalPointer() const;
154 caltinay 3594
155 caltinay 5929 inline void startCollect(const double* in) const
156 caltinay 4836 {
157     startColCollect(in);
158     }
159 caltinay 3594
160 caltinay 5929 inline double* finishCollect() const
161 caltinay 4836 {
162     return finishColCollect();
163     }
164 caltinay 3594
165 caltinay 5929 inline void startColCollect(const double* in) const
166 caltinay 4836 {
167     col_coupler->startCollect(in);
168     }
169 caltinay 3594
170 caltinay 5929 inline double* finishColCollect() const
171 caltinay 4836 {
172     return col_coupler->finishCollect();
173     }
174 caltinay 3594
175 caltinay 4836 inline void startRowCollect(const double* in)
176     {
177     row_coupler->startCollect(in);
178     }
179 caltinay 3594
180 caltinay 4836 inline double* finishRowCollect()
181     {
182     return row_coupler->finishCollect();
183     }
184 caltinay 3594
185 caltinay 4836 inline dim_t getNumRows() const
186     {
187     return mainBlock->numRows;
188     }
189 caltinay 3594
190 caltinay 4836 inline dim_t getNumCols() const
191     {
192     return mainBlock->numCols;
193     }
194 gross 3369
195 caltinay 4836 inline dim_t getTotalNumRows() const
196     {
197     return getNumRows() * row_block_size;
198     }
199 caltinay 3594
200 caltinay 4836 inline dim_t getTotalNumCols() const
201     {
202     return getNumCols() * col_block_size;
203     }
204 caltinay 3594
205 caltinay 4836 inline dim_t getRowOverlap() const
206     {
207     return row_coupler->getNumOverlapComponents();
208     }
209 caltinay 3594
210 caltinay 4836 inline dim_t getColOverlap() const
211     {
212     return col_coupler->getNumOverlapComponents();
213     }
214 caltinay 3594
215 caltinay 4836 inline dim_t getGlobalNumRows() const
216     {
217     if (type & MATRIX_FORMAT_CSC) {
218     return pattern->input_distribution->getGlobalNumComponents();
219     }
220     return pattern->output_distribution->getGlobalNumComponents();
221     }
222 caltinay 3594
223 caltinay 4836 inline dim_t getGlobalNumCols() const
224     {
225     if (type & MATRIX_FORMAT_CSC) {
226     return pattern->output_distribution->getGlobalNumComponents();
227     }
228     return pattern->input_distribution->getGlobalNumComponents();
229     }
230 caltinay 3594
231 caltinay 4836 inline dim_t getGlobalTotalNumRows() const
232     {
233     return getGlobalNumRows() * row_block_size;
234     }
235 caltinay 3594
236 caltinay 4836 inline dim_t getGlobalTotalNumCols() const
237     {
238     return getGlobalNumCols() * col_block_size;
239     }
240 caltinay 3594
241 caltinay 4836 inline double getSparsity() const
242     {
243     return getGlobalSize() /
244 caltinay 4869 ((double)getGlobalTotalNumRows()*getGlobalTotalNumCols());
245 caltinay 4836 }
246 caltinay 3594
247 caltinay 4836 inline dim_t getNumOutput() const
248     {
249     return pattern->getNumOutput();
250     }
251 caltinay 3594
252 caltinay 4836 inline void copyBlockFromMainDiagonal(double* out) const
253     {
254     mainBlock->copyBlockFromMainDiagonal(out);
255     }
256 caltinay 3594
257 caltinay 4836 inline void copyBlockToMainDiagonal(const double* in)
258     {
259     mainBlock->copyBlockToMainDiagonal(in);
260     }
261 jgs 150
262 caltinay 4836 inline void copyFromMainDiagonal(double* out) const
263     {
264     mainBlock->copyFromMainDiagonal(out);
265     }
266 gross 3445
267 caltinay 4836 inline void copyToMainDiagonal(const double* in)
268     {
269     mainBlock->copyToMainDiagonal(in);
270     }
271 gross 3445
272 caltinay 4836 inline void setValues(double value)
273     {
274     mainBlock->setValues(value);
275     col_coupleBlock->setValues(value);
276     row_coupleBlock->setValues(value);
277     is_balanced = false;
278     }
279 caltinay 3594
280 caltinay 4836 inline void rowSum(double* row_sum) const
281     {
282     if ((type & MATRIX_FORMAT_CSC) || (type & MATRIX_FORMAT_OFFSET1)) {
283 caltinay 5996 throw PasoException("SystemMatrix::rowSum: No normalization "
284 caltinay 4836 "available for compressed sparse column or index offset 1.");
285     } else {
286     const dim_t nrow = mainBlock->numRows*row_block_size;
287     #pragma omp parallel for
288     for (index_t irow=0; irow<nrow; ++irow) {
289     row_sum[irow]=0.;
290     }
291     mainBlock->addRow_CSR_OFFSET0(row_sum);
292     col_coupleBlock->addRow_CSR_OFFSET0(row_sum);
293     }
294     }
295 caltinay 3594
296 caltinay 5929 void MatrixVector(double alpha, const double* in, double beta,
297     double* out) const;
298    
299     void MatrixVector_CSR_OFFSET0(double alpha, const double* in, double beta,
300     double* out) const;
301    
302 caltinay 4836 static SystemMatrix_ptr loadMM_toCSR(const char* filename);
303 caltinay 3594
304 caltinay 4836 static SystemMatrix_ptr loadMM_toCSC(const char* filename);
305 gross 3369
306 caltinay 5929 static int getSystemMatrixTypeId(int solver, int preconditioner,
307     int package, bool symmetry,
308 caltinay 5997 const escript::JMPI& mpi_info);
309 caltinay 3594
310 caltinay 4836 SystemMatrixType type;
311     SystemMatrixPattern_ptr pattern;
312 caltinay 3594
313 caltinay 4836 dim_t logical_row_block_size;
314     dim_t logical_col_block_size;
315 gross 3369
316 caltinay 4836 dim_t row_block_size;
317     dim_t col_block_size;
318     dim_t block_size;
319 caltinay 3594
320 caltinay 6197 escript::Distribution_ptr row_distribution;
321     escript::Distribution_ptr col_distribution;
322 caltinay 5997 escript::JMPI mpi_info;
323 caltinay 3594
324 caltinay 6581 Coupler_ptr<real_t> col_coupler;
325     Coupler_ptr<real_t> row_coupler;
326 caltinay 3594
327 caltinay 4836 /// main block
328     SparseMatrix_ptr mainBlock;
329     /// coupling to neighbouring processors (row - col)
330     SparseMatrix_ptr col_coupleBlock;
331     /// coupling to neighbouring processors (col - row)
332     SparseMatrix_ptr row_coupleBlock;
333     /// coupling of rows-cols on neighbouring processors (may not be valid)
334     SparseMatrix_ptr remote_coupleBlock;
335 caltinay 3594
336 caltinay 4836 bool is_balanced;
337 jgs 150
338 caltinay 4836 /// matrix may be balanced by a diagonal matrix D=diagonal(balance_vector)
339     /// if is_balanced is true, the matrix stored is D*A*D where A is the
340     /// original matrix.
341     /// When the system of linear equations is solved we solve D*A*D*y=c.
342     /// So to solve A*x=b one needs to set c=D*b and x=D*y.
343     double* balance_vector;
344 caltinay 3594
345 caltinay 4836 /// stores the global ids for all cols in col_coupleBlock
346 caltinay 4849 mutable index_t* global_id;
347 gross 3369
348 caltinay 4836 /// package code controlling the solver pointer
349 caltinay 5929 mutable index_t solver_package;
350 caltinay 3594
351 caltinay 4836 /// pointer to data needed by a solver
352     void* solver_p;
353 caltinay 3594
354 caltinay 5929 private:
355     virtual void setToSolution(escript::Data& out, escript::Data& in,
356     boost::python::object& options) const;
357 caltinay 3594
358 caltinay 5929 virtual void ypAx(escript::Data& y, escript::Data& x) const;
359 caltinay 3594
360 caltinay 5929 void solve(double* out, double* in, Options* options) const;
361     };
362 gross 1562
363 caltinay 3594
364 caltinay 4869 void RHS_loadMM_toCSR(const char* filename, double* b, dim_t size);
365 caltinay 3594
366    
367 caltinay 4836 } // namespace paso
368 caltinay 4873
369 caltinay 4836 #endif // __PASO_SYSTEMMATRIX_H__
370 ksteube 1312

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26