/[escript]/branches/trilinos_from_5897/paso/src/SystemMatrix.h
ViewVC logotype

Contents of /branches/trilinos_from_5897/paso/src/SystemMatrix.h

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6058 - (show annotations)
Thu Mar 10 06:51:55 2016 UTC (3 years, 1 month ago) by caltinay
File MIME type: text/plain
File size: 11155 byte(s)
added getPtr() to AbstractSystemMatrix so we can now use shared systemmatrix
pointers rather than circumventing them.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
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 class 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 } else {
73 mainBlock->saveMM(filename.c_str());
74 }
75 }
76
77 virtual inline void saveHB(const std::string& filename) const
78 {
79 if (mpi_info->size > 1) {
80 throw PasoException("SystemMatrix::saveHB: Only single rank supported.");
81 } else if (!(type & MATRIX_FORMAT_CSC)) {
82 throw PasoException("SystemMatrix::saveHB: Only CSC format supported.");
83 } else {
84 mainBlock->saveHB_CSC(filename.c_str());
85 }
86 }
87
88 virtual void resetValues();
89
90 /// Nullifies rows in the matrix.
91 /// The rows are marked by positive values in mask_row. Values on the
92 /// main diagonal which are marked to set to zero by mask_row are set
93 /// to main_diagonal_value.
94 void nullifyRows(double* mask_row, double main_diagonal_value);
95
96 void add(dim_t, index_t*, dim_t, dim_t, index_t*, dim_t, double*);
97
98 void makeZeroRowSums(double* left_over);
99
100 /// copies the col_coupleBlock into row_coupleBlock.
101 /// WARNING: this method uses mpi_requests of the coupler attached to the
102 /// matrix. No reordering on the received columns is performed.
103 /// In practice this means that components in
104 /// row_coupleBlock->pattern->index and
105 /// row_coupler->connector->recv->shared
106 /// are ordered by increasing value.
107 /// Note that send and receive row_coupler->connectors are swapping roles.
108 void copyColCoupleBlock();
109
110 void copyRemoteCoupleBlock(bool recreatePattern);
111
112 void fillWithGlobalCoordinates(double f1);
113
114 void print() const;
115
116 /// Merges the system matrix which is distributed on several MPI ranks
117 /// into a complete sparse matrix on rank 0. Used by the Merged Solver.
118 SparseMatrix_ptr mergeSystemMatrix() const;
119
120 void mergeMainAndCouple(index_t** p_ptr, index_t** p_idx, double** p_val) const;
121
122 void mergeMainAndCouple_CSR_OFFSET0(index_t** p_ptr, index_t** p_idx, double** p_val) const;
123 void mergeMainAndCouple_CSR_OFFSET0_Block(index_t** p_ptr, index_t** p_idx, double** p_val) const;
124
125 void mergeMainAndCouple_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val) const;
126
127 void copyMain_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val);
128
129 void extendedRowsForST(dim_t* degree_ST, index_t* offset_ST, index_t* ST);
130
131 void applyBalanceInPlace(double* x, bool RHS) const;
132
133 void applyBalance(double* x_out, const double* x, bool RHS) const;
134
135 void balance();
136
137 double getGlobalSize() const;
138
139 void setPreconditioner(Options* options);
140
141 /// Applies the preconditioner.
142 /// This method needs to be called within a parallel region.
143 /// Barrier synchronization is performed before the evaluation to make
144 /// sure that the input vector is available
145 void solvePreconditioner(double* x, double* b);
146
147 void freePreconditioner();
148
149 index_t* borrowMainDiagonalPointer() const;
150
151 inline void startCollect(const double* in) const
152 {
153 startColCollect(in);
154 }
155
156 inline double* finishCollect() const
157 {
158 return finishColCollect();
159 }
160
161 inline void startColCollect(const double* in) const
162 {
163 col_coupler->startCollect(in);
164 }
165
166 inline double* finishColCollect() const
167 {
168 return col_coupler->finishCollect();
169 }
170
171 inline void startRowCollect(const double* in)
172 {
173 row_coupler->startCollect(in);
174 }
175
176 inline double* finishRowCollect()
177 {
178 return row_coupler->finishCollect();
179 }
180
181 inline dim_t getNumRows() const
182 {
183 return mainBlock->numRows;
184 }
185
186 inline dim_t getNumCols() const
187 {
188 return mainBlock->numCols;
189 }
190
191 inline dim_t getTotalNumRows() const
192 {
193 return getNumRows() * row_block_size;
194 }
195
196 inline dim_t getTotalNumCols() const
197 {
198 return getNumCols() * col_block_size;
199 }
200
201 inline dim_t getRowOverlap() const
202 {
203 return row_coupler->getNumOverlapComponents();
204 }
205
206 inline dim_t getColOverlap() const
207 {
208 return col_coupler->getNumOverlapComponents();
209 }
210
211 inline dim_t getGlobalNumRows() const
212 {
213 if (type & MATRIX_FORMAT_CSC) {
214 return pattern->input_distribution->getGlobalNumComponents();
215 }
216 return pattern->output_distribution->getGlobalNumComponents();
217 }
218
219 inline dim_t getGlobalNumCols() const
220 {
221 if (type & MATRIX_FORMAT_CSC) {
222 return pattern->output_distribution->getGlobalNumComponents();
223 }
224 return pattern->input_distribution->getGlobalNumComponents();
225 }
226
227 inline dim_t getGlobalTotalNumRows() const
228 {
229 return getGlobalNumRows() * row_block_size;
230 }
231
232 inline dim_t getGlobalTotalNumCols() const
233 {
234 return getGlobalNumCols() * col_block_size;
235 }
236
237 inline double getSparsity() const
238 {
239 return getGlobalSize() /
240 ((double)getGlobalTotalNumRows()*getGlobalTotalNumCols());
241 }
242
243 inline dim_t getNumOutput() const
244 {
245 return pattern->getNumOutput();
246 }
247
248 inline void copyBlockFromMainDiagonal(double* out) const
249 {
250 mainBlock->copyBlockFromMainDiagonal(out);
251 }
252
253 inline void copyBlockToMainDiagonal(const double* in)
254 {
255 mainBlock->copyBlockToMainDiagonal(in);
256 }
257
258 inline void copyFromMainDiagonal(double* out) const
259 {
260 mainBlock->copyFromMainDiagonal(out);
261 }
262
263 inline void copyToMainDiagonal(const double* in)
264 {
265 mainBlock->copyToMainDiagonal(in);
266 }
267
268 inline void setValues(double value)
269 {
270 mainBlock->setValues(value);
271 col_coupleBlock->setValues(value);
272 row_coupleBlock->setValues(value);
273 is_balanced = false;
274 }
275
276 inline void rowSum(double* row_sum) const
277 {
278 if ((type & MATRIX_FORMAT_CSC) || (type & MATRIX_FORMAT_OFFSET1)) {
279 throw PasoException("SystemMatrix::rowSum: No normalization "
280 "available for compressed sparse column or index offset 1.");
281 } else {
282 const dim_t nrow = mainBlock->numRows*row_block_size;
283 #pragma omp parallel for
284 for (index_t irow=0; irow<nrow; ++irow) {
285 row_sum[irow]=0.;
286 }
287 mainBlock->addRow_CSR_OFFSET0(row_sum);
288 col_coupleBlock->addRow_CSR_OFFSET0(row_sum);
289 }
290 }
291
292 void MatrixVector(double alpha, const double* in, double beta,
293 double* out) const;
294
295 void MatrixVector_CSR_OFFSET0(double alpha, const double* in, double beta,
296 double* out) const;
297
298 static SystemMatrix_ptr loadMM_toCSR(const char* filename);
299
300 static SystemMatrix_ptr loadMM_toCSC(const char* filename);
301
302 static int getSystemMatrixTypeId(int solver, int preconditioner,
303 int package, bool symmetry,
304 const escript::JMPI& mpi_info);
305
306 SystemMatrixType type;
307 SystemMatrixPattern_ptr pattern;
308
309 dim_t logical_row_block_size;
310 dim_t logical_col_block_size;
311
312 dim_t row_block_size;
313 dim_t col_block_size;
314 dim_t block_size;
315
316 Distribution_ptr row_distribution;
317 Distribution_ptr col_distribution;
318 escript::JMPI mpi_info;
319
320 Coupler_ptr col_coupler;
321 Coupler_ptr row_coupler;
322
323 /// main block
324 SparseMatrix_ptr mainBlock;
325 /// coupling to neighbouring processors (row - col)
326 SparseMatrix_ptr col_coupleBlock;
327 /// coupling to neighbouring processors (col - row)
328 SparseMatrix_ptr row_coupleBlock;
329 /// coupling of rows-cols on neighbouring processors (may not be valid)
330 SparseMatrix_ptr remote_coupleBlock;
331
332 bool is_balanced;
333
334 /// matrix may be balanced by a diagonal matrix D=diagonal(balance_vector)
335 /// if is_balanced is true, the matrix stored is D*A*D where A is the
336 /// original matrix.
337 /// When the system of linear equations is solved we solve D*A*D*y=c.
338 /// So to solve A*x=b one needs to set c=D*b and x=D*y.
339 double* balance_vector;
340
341 /// stores the global ids for all cols in col_coupleBlock
342 mutable index_t* global_id;
343
344 /// package code controlling the solver pointer
345 mutable index_t solver_package;
346
347 /// pointer to data needed by a solver
348 void* solver_p;
349
350 private:
351 virtual void setToSolution(escript::Data& out, escript::Data& in,
352 boost::python::object& options) const;
353
354 virtual void ypAx(escript::Data& y, escript::Data& x) const;
355
356 void solve(double* out, double* in, Options* options) const;
357 };
358
359
360 void RHS_loadMM_toCSR(const char* filename, double* b, dim_t size);
361
362
363 } // namespace paso
364
365 #endif // __PASO_SYSTEMMATRIX_H__
366

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26