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

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

Parent Directory Parent Directory | Revision Log Revision Log


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


1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2018 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Apache License, version 2.0
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************/
19
20 /* Paso: SystemMatrix */
21
22 /****************************************************************************/
23
24 /* Copyrights by ACcESS Australia 2003,2004,2005,2006 */
25 /* Author: Lutz Gross, l.gross@uq.edu.au */
26
27 /****************************************************************************/
28
29 #ifndef __PASO_SYSTEMMATRIX_H__
30 #define __PASO_SYSTEMMATRIX_H__
31
32 #include "SparseMatrix.h"
33 #include "SystemMatrixPattern.h"
34
35 #include <escript/AbstractSystemMatrix.h>
36
37 namespace paso {
38
39 struct Options;
40 class SystemMatrix;
41 typedef boost::shared_ptr<SystemMatrix> SystemMatrix_ptr;
42 typedef boost::shared_ptr<const SystemMatrix> const_SystemMatrix_ptr;
43
44 typedef int SystemMatrixType;
45
46 /// this class holds a (distributed) stiffness matrix
47 class SystemMatrix : public escript::AbstractSystemMatrix
48 {
49 public:
50 /// default constructor - throws exception.
51 SystemMatrix();
52
53 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 ~SystemMatrix();
59
60 /// 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 virtual void nullifyRowsAndCols(escript::Data& mask_row,
65 escript::Data& mask_col,
66 double main_diagonal_value);
67
68 virtual inline void saveMM(const std::string& filename) const
69 {
70 if (mpi_info->size > 1) {
71 //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 } 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 throw PasoException("SystemMatrix::saveHB: Only single rank supported.");
84 } else if (!(type & MATRIX_FORMAT_CSC)) {
85 throw PasoException("SystemMatrix::saveHB: Only CSC format supported.");
86 } else {
87 mainBlock->saveHB_CSC(filename.c_str());
88 }
89 }
90
91 virtual void resetValues(bool preserveSolverData = false);
92
93 /// 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
99 void add(dim_t, index_t*, dim_t, dim_t, index_t*, dim_t, double*);
100
101 void makeZeroRowSums(double* left_over);
102
103 /// copies the col_coupleBlock into row_coupleBlock.
104 /// 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
113 void copyRemoteCoupleBlock(bool recreatePattern);
114
115 void fillWithGlobalCoordinates(double f1);
116
117 void print() const;
118
119 /// 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
123 void mergeMainAndCouple(index_t** p_ptr, index_t** p_idx, double** p_val) const;
124
125 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
128 void mergeMainAndCouple_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val) const;
129
130 void copyMain_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val);
131
132 void extendedRowsForST(dim_t* degree_ST, index_t* offset_ST, index_t* ST);
133
134 void applyBalanceInPlace(double* x, bool RHS) const;
135
136 void applyBalance(double* x_out, const double* x, bool RHS) const;
137
138 void balance();
139
140 double getGlobalSize() const;
141
142 void setPreconditioner(Options* options);
143
144 /// Applies the preconditioner.
145 /// This method needs to be called within a parallel region.
146 /// 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
150 void freePreconditioner();
151
152 index_t* borrowMainDiagonalPointer() const;
153
154 inline void startCollect(const double* in) const
155 {
156 startColCollect(in);
157 }
158
159 inline double* finishCollect() const
160 {
161 return finishColCollect();
162 }
163
164 inline void startColCollect(const double* in) const
165 {
166 col_coupler->startCollect(in);
167 }
168
169 inline double* finishColCollect() const
170 {
171 return col_coupler->finishCollect();
172 }
173
174 inline void startRowCollect(const double* in)
175 {
176 row_coupler->startCollect(in);
177 }
178
179 inline double* finishRowCollect()
180 {
181 return row_coupler->finishCollect();
182 }
183
184 inline dim_t getNumRows() const
185 {
186 return mainBlock->numRows;
187 }
188
189 inline dim_t getNumCols() const
190 {
191 return mainBlock->numCols;
192 }
193
194 inline dim_t getTotalNumRows() const
195 {
196 return getNumRows() * row_block_size;
197 }
198
199 inline dim_t getTotalNumCols() const
200 {
201 return getNumCols() * col_block_size;
202 }
203
204 inline dim_t getRowOverlap() const
205 {
206 return row_coupler->getNumOverlapComponents();
207 }
208
209 inline dim_t getColOverlap() const
210 {
211 return col_coupler->getNumOverlapComponents();
212 }
213
214 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
222 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
230 inline dim_t getGlobalTotalNumRows() const
231 {
232 return getGlobalNumRows() * row_block_size;
233 }
234
235 inline dim_t getGlobalTotalNumCols() const
236 {
237 return getGlobalNumCols() * col_block_size;
238 }
239
240 inline double getSparsity() const
241 {
242 return getGlobalSize() /
243 ((double)getGlobalTotalNumRows()*getGlobalTotalNumCols());
244 }
245
246 inline dim_t getNumOutput() const
247 {
248 return pattern->getNumOutput();
249 }
250
251 inline void copyBlockFromMainDiagonal(double* out) const
252 {
253 mainBlock->copyBlockFromMainDiagonal(out);
254 }
255
256 inline void copyBlockToMainDiagonal(const double* in)
257 {
258 mainBlock->copyBlockToMainDiagonal(in);
259 }
260
261 inline void copyFromMainDiagonal(double* out) const
262 {
263 mainBlock->copyFromMainDiagonal(out);
264 }
265
266 inline void copyToMainDiagonal(const double* in)
267 {
268 mainBlock->copyToMainDiagonal(in);
269 }
270
271 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
279 inline void rowSum(double* row_sum) const
280 {
281 if ((type & MATRIX_FORMAT_CSC) || (type & MATRIX_FORMAT_OFFSET1)) {
282 throw PasoException("SystemMatrix::rowSum: No normalization "
283 "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
295 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 static SystemMatrix_ptr loadMM_toCSR(const char* filename);
302
303 static SystemMatrix_ptr loadMM_toCSC(const char* filename);
304
305 static int getSystemMatrixTypeId(int solver, int preconditioner,
306 int package, bool symmetry,
307 const escript::JMPI& mpi_info);
308
309 SystemMatrixType type;
310 SystemMatrixPattern_ptr pattern;
311
312 dim_t logical_row_block_size;
313 dim_t logical_col_block_size;
314
315 dim_t row_block_size;
316 dim_t col_block_size;
317 dim_t block_size;
318
319 escript::Distribution_ptr row_distribution;
320 escript::Distribution_ptr col_distribution;
321 escript::JMPI mpi_info;
322
323 Coupler_ptr<real_t> col_coupler;
324 Coupler_ptr<real_t> row_coupler;
325
326 /// 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
335 bool is_balanced;
336
337 /// 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
344 /// stores the global ids for all cols in col_coupleBlock
345 mutable index_t* global_id;
346
347 /// package code controlling the solver pointer
348 mutable index_t solver_package;
349
350 /// pointer to data needed by a solver
351 void* solver_p;
352
353 private:
354 virtual void setToSolution(escript::Data& out, escript::Data& in,
355 boost::python::object& options) const;
356
357 virtual void ypAx(escript::Data& y, escript::Data& x) const;
358
359 void solve(double* out, double* in, Options* options) const;
360 };
361
362
363 void RHS_loadMM_toCSR(const char* filename, double* b, dim_t size);
364
365
366 } // namespace paso
367
368 #endif // __PASO_SYSTEMMATRIX_H__
369

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26