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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6726 - (hide annotations)
Mon Oct 8 04:44:29 2018 UTC (12 months, 1 week ago) by aellery
File MIME type: text/plain
File size: 11358 byte(s)
Fixed bug #435 (compile warnings under clang++).


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26