escript  Revision_Unversioneddirectory
SystemMatrix.h
Go to the documentation of this file.
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 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 #include "Options.h"
35 
36 namespace paso {
37 
38 struct SystemMatrix;
39 typedef boost::shared_ptr<SystemMatrix> SystemMatrix_ptr;
40 typedef boost::shared_ptr<const SystemMatrix> const_SystemMatrix_ptr;
41 
42 typedef int SystemMatrixType;
43 
44 // this struct holds a (distributed) stiffness matrix
46 struct SystemMatrix : boost::enable_shared_from_this<SystemMatrix>
47 {
48  SystemMatrix(SystemMatrixType, SystemMatrixPattern_ptr, dim_t, dim_t,
49  bool patternIsUnrolled);
50 
51  ~SystemMatrix();
52 
57  void nullifyRowsAndCols(double* mask_row, double* mask_col,
58  double main_diagonal_value);
59 
64  void nullifyRows(double* mask_row, double main_diagonal_value);
65 
66  void add(dim_t, index_t*, dim_t, dim_t, index_t*, dim_t, double*);
67 
68  void makeZeroRowSums(double* left_over);
69 
78  void copyColCoupleBlock();
79 
80  void copyRemoteCoupleBlock(bool recreatePattern);
81 
82  void fillWithGlobalCoordinates(double f1);
83 
84  void print() const;
85 
89 
90  void mergeMainAndCouple(index_t** p_ptr, index_t** p_idx, double** p_val) const;
91 
92  void mergeMainAndCouple_CSR_OFFSET0(index_t** p_ptr, index_t** p_idx, double** p_val) const;
93  void mergeMainAndCouple_CSR_OFFSET0_Block(index_t** p_ptr, index_t** p_idx, double** p_val) const;
94 
95  void mergeMainAndCouple_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val) const;
96 
97  void copyMain_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val);
98 
99  void extendedRowsForST(dim_t* degree_ST, index_t* offset_ST, index_t* ST);
100 
101  void applyBalanceInPlace(double* x, bool RHS) const;
102 
103  void applyBalance(double* x_out, const double* x, bool RHS) const;
104 
105  void balance();
106 
107  double getGlobalSize() const;
108 
109  void setPreconditioner(Options* options);
110 
115  void solvePreconditioner(double* x, double* b);
116 
117  void freePreconditioner();
118 
120 
121  inline void startCollect(const double* in)
122  {
123  startColCollect(in);
124  }
125 
126  inline double* finishCollect()
127  {
128  return finishColCollect();
129  }
130 
131  inline void startColCollect(const double* in)
132  {
133  col_coupler->startCollect(in);
134  }
135 
136  inline double* finishColCollect()
137  {
138  return col_coupler->finishCollect();
139  }
140 
141  inline void startRowCollect(const double* in)
142  {
143  row_coupler->startCollect(in);
144  }
145 
146  inline double* finishRowCollect()
147  {
148  return row_coupler->finishCollect();
149  }
150 
151  inline dim_t getNumRows() const
152  {
153  return mainBlock->numRows;
154  }
155 
156  inline dim_t getNumCols() const
157  {
158  return mainBlock->numCols;
159  }
160 
161  inline dim_t getTotalNumRows() const
162  {
163  return getNumRows() * row_block_size;
164  }
165 
166  inline dim_t getTotalNumCols() const
167  {
168  return getNumCols() * col_block_size;
169  }
170 
171  inline dim_t getRowOverlap() const
172  {
173  return row_coupler->getNumOverlapComponents();
174  }
175 
176  inline dim_t getColOverlap() const
177  {
178  return col_coupler->getNumOverlapComponents();
179  }
180 
181  inline dim_t getGlobalNumRows() const
182  {
183  if (type & MATRIX_FORMAT_CSC) {
184  return pattern->input_distribution->getGlobalNumComponents();
185  }
186  return pattern->output_distribution->getGlobalNumComponents();
187  }
188 
189  inline dim_t getGlobalNumCols() const
190  {
191  if (type & MATRIX_FORMAT_CSC) {
192  return pattern->output_distribution->getGlobalNumComponents();
193  }
194  return pattern->input_distribution->getGlobalNumComponents();
195  }
196 
198  {
199  return getGlobalNumRows() * row_block_size;
200  }
201 
203  {
204  return getGlobalNumCols() * col_block_size;
205  }
206 
207  inline double getSparsity() const
208  {
209  return getGlobalSize() /
211  }
212 
213  inline dim_t getNumOutput() const
214  {
215  return pattern->getNumOutput();
216  }
217 
218  inline void copyBlockFromMainDiagonal(double* out) const
219  {
220  mainBlock->copyBlockFromMainDiagonal(out);
221  }
222 
223  inline void copyBlockToMainDiagonal(const double* in)
224  {
225  mainBlock->copyBlockToMainDiagonal(in);
226  }
227 
228  inline void copyFromMainDiagonal(double* out) const
229  {
230  mainBlock->copyFromMainDiagonal(out);
231  }
232 
233  inline void copyToMainDiagonal(const double* in)
234  {
235  mainBlock->copyToMainDiagonal(in);
236  }
237 
238  inline void setValues(double value)
239  {
240  mainBlock->setValues(value);
241  col_coupleBlock->setValues(value);
242  row_coupleBlock->setValues(value);
243  is_balanced = false;
244  }
245 
246  inline void saveMM(const char* filename) const
247  {
248  if (mpi_info->size > 1) {
249  Esys_setError(IO_ERROR, "SystemMatrix::saveMM: Only single rank supported.");
250  } else {
251  mainBlock->saveMM(filename);
252  }
253  }
254 
255  inline void saveHB(const char *filename) const
256  {
257  if (mpi_info->size > 1) {
258  Esys_setError(TYPE_ERROR, "SystemMatrix::saveHB: Only single rank supported.");
259  } else if (!(type & MATRIX_FORMAT_CSC)) {
260  Esys_setError(TYPE_ERROR, "SystemMatrix::saveHB: Only CSC format supported.");
261  } else {
262  mainBlock->saveHB_CSC(filename);
263  }
264  }
265 
266  inline void rowSum(double* row_sum) const
267  {
269  Esys_setError(TYPE_ERROR, "SystemMatrix::rowSum: No normalization "
270  "available for compressed sparse column or index offset 1.");
271  } else {
272  const dim_t nrow = mainBlock->numRows*row_block_size;
273 #pragma omp parallel for
274  for (index_t irow=0; irow<nrow; ++irow) {
275  row_sum[irow]=0.;
276  }
277  mainBlock->addRow_CSR_OFFSET0(row_sum);
278  col_coupleBlock->addRow_CSR_OFFSET0(row_sum);
279  }
280  }
281 
282  static SystemMatrix_ptr loadMM_toCSR(const char* filename);
283 
284  static SystemMatrix_ptr loadMM_toCSC(const char* filename);
285 
286  static index_t getSystemMatrixTypeId(index_t solver,
287  index_t preconditioner,
288  index_t package, bool symmetry,
289  const esysUtils::JMPI& mpi_info);
290 
291  SystemMatrixType type;
293 
296 
300 
304 
307 
316 
318 
324  double* balance_vector;
325 
327  mutable index_t* global_id;
328 
331 
333  void* solver_p;
334 
337 };
338 
339 
340 void SystemMatrix_MatrixVector(double alpha, SystemMatrix_ptr A, const double* in, double beta, double* out);
341 
342 void SystemMatrix_MatrixVector_CSR_OFFSET0(double alpha, SystemMatrix_ptr A, const double* in, double beta, double* out);
343 
344 void RHS_loadMM_toCSR(const char* filename, double* b, dim_t size);
345 
346 
347 } // namespace paso
348 
349 #endif // __PASO_SYSTEMMATRIX_H__
350 
void copyToMainDiagonal(const double *in)
Definition: SystemMatrix.h:233
#define PASO_DLL_API
Definition: Paso.h:41
Definition: error.h:48
void nullifyRowsAndCols(double *mask_row, double *mask_col, double main_diagonal_value)
Definition: SystemMatrix.cpp:243
Distribution_ptr col_distribution
Definition: SystemMatrix.h:302
SparseMatrix_ptr row_coupleBlock
coupling to neighbouring processors (col - row)
Definition: SystemMatrix.h:313
void mergeMainAndCouple_CSR_OFFSET0_Block(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:173
SparseMatrix_ptr mergeSystemMatrix() const
Definition: SystemMatrix.cpp:500
SparseMatrix_ptr mainBlock
main block
Definition: SystemMatrix.h:309
void copyBlockFromMainDiagonal(double *out) const
Definition: SystemMatrix.h:218
dim_t getTotalNumRows() const
Definition: SystemMatrix.h:161
Definition: Options.h:90
dim_t logical_col_block_size
Definition: SystemMatrix.h:295
void startColCollect(const double *in)
Definition: SystemMatrix.h:131
void * trilinos_data
this is only used for a trilinos matrix
Definition: SystemMatrix.h:336
dim_t logical_row_block_size
Definition: SystemMatrix.h:294
void * solver_p
pointer to data needed by a solver
Definition: SystemMatrix.h:333
double * finishCollect()
Definition: SystemMatrix.h:126
#define MATRIX_FORMAT_OFFSET1
Definition: Paso.h:56
double * finishRowCollect()
Definition: SystemMatrix.h:146
void applyBalance(double *x_out, const double *x, bool RHS) const
Definition: SystemMatrix.cpp:390
SystemMatrixPattern_ptr pattern
Definition: SystemMatrix.h:292
void mergeMainAndCouple_CSR_OFFSET0(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:59
Distribution_ptr row_distribution
Definition: SystemMatrix.h:301
dim_t getNumRows() const
Definition: SystemMatrix.h:151
index_t * global_id
stores the global ids for all cols in col_coupleBlock
Definition: SystemMatrix.h:327
static SystemMatrix_ptr loadMM_toCSR(const char *filename)
Definition: SystemMatrix_loadMM.cpp:103
void Esys_setError(Esys_ErrorCodeType err, __const char *msg)
Definition: error.cpp:38
dim_t getGlobalNumCols() const
Definition: SystemMatrix.h:189
bool is_balanced
Definition: SystemMatrix.h:317
void startCollect(const double *in)
Definition: SystemMatrix.h:121
dim_t getGlobalNumRows() const
Definition: SystemMatrix.h:181
void freePreconditioner()
Definition: SystemMatrix.cpp:161
void saveHB(const char *filename) const
Definition: SystemMatrix.h:255
dim_t getGlobalTotalNumCols() const
Definition: SystemMatrix.h:202
boost::shared_ptr< SparseMatrix > SparseMatrix_ptr
Definition: SparseMatrix.h:35
boost::shared_ptr< Distribution > Distribution_ptr
Definition: Distribution.h:36
void balance()
Definition: SystemMatrix.cpp:409
SystemMatrixType type
Definition: SystemMatrix.h:291
dim_t getNumOutput() const
Definition: SystemMatrix.h:213
void rowSum(double *row_sum) const
Definition: SystemMatrix.h:266
void extendedRowsForST(dim_t *degree_ST, index_t *offset_ST, index_t *ST)
Definition: SystemMatrix_extendedRows.cpp:41
void SystemMatrix_MatrixVector_CSR_OFFSET0(double alpha, SystemMatrix_ptr A, const double *in, double beta, double *out)
Definition: SystemMatrix_MatrixVector.cpp:72
dim_t getRowOverlap() const
Definition: SystemMatrix.h:171
boost::shared_ptr< SystemMatrixPattern > SystemMatrixPattern_ptr
Definition: SystemMatrixPattern.h:38
void RHS_loadMM_toCSR(const char *filename, double *b, dim_t size)
Definition: SystemMatrix_loadMM.cpp:314
boost::shared_ptr< SystemMatrix > SystemMatrix_ptr
Definition: SystemMatrix.h:38
static index_t getSystemMatrixTypeId(index_t solver, index_t preconditioner, index_t package, bool symmetry, const esysUtils::JMPI &mpi_info)
Definition: SystemMatrix.cpp:461
Definition: AMG.cpp:38
dim_t getTotalNumCols() const
Definition: SystemMatrix.h:166
SparseMatrix_ptr col_coupleBlock
coupling to neighbouring processors (row - col)
Definition: SystemMatrix.h:311
~SystemMatrix()
Definition: SystemMatrix.cpp:141
void fillWithGlobalCoordinates(double f1)
Definition: SystemMatrix_debug.cpp:36
void setValues(double value)
Definition: SystemMatrix.h:238
double getSparsity() const
Definition: SystemMatrix.h:207
Definition: SystemMatrix.h:46
void SystemMatrix_MatrixVector(double alpha, SystemMatrix_ptr A, const double *in, double beta, double *out)
Definition: SystemMatrix_MatrixVector.cpp:35
void startRowCollect(const double *in)
Definition: SystemMatrix.h:141
Coupler_ptr col_coupler
Definition: SystemMatrix.h:305
dim_t row_block_size
Definition: SystemMatrix.h:297
dim_t block_size
Definition: SystemMatrix.h:299
SystemMatrix(SystemMatrixType, SystemMatrixPattern_ptr, dim_t, dim_t, bool patternIsUnrolled)
Definition: SystemMatrix.cpp:42
void copyRemoteCoupleBlock(bool recreatePattern)
Definition: SystemMatrix_copyRemoteCoupleBlock.cpp:39
double getGlobalSize() const
Definition: SystemMatrix.cpp:168
boost::shared_ptr< const SystemMatrix > const_SystemMatrix_ptr
Definition: SystemMatrix.h:40
void print() const
Definition: SystemMatrix_debug.cpp:100
int SystemMatrixType
Definition: SystemMatrix.h:42
double * balance_vector
Definition: SystemMatrix.h:324
void copyFromMainDiagonal(double *out) const
Definition: SystemMatrix.h:228
void add(dim_t, index_t *, dim_t, dim_t, index_t *, dim_t, double *)
void mergeMainAndCouple_CSC_OFFSET1(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:283
void nullifyRows(double *mask_row, double main_diagonal_value)
Definition: SystemMatrix.cpp:220
int index_t
Definition: types.h:24
dim_t getNumCols() const
Definition: SystemMatrix.h:156
static SystemMatrix_ptr loadMM_toCSC(const char *filename)
Definition: SystemMatrix_loadMM.cpp:209
dim_t col_block_size
Definition: SystemMatrix.h:298
index_t * borrowMainDiagonalPointer() const
Definition: SystemMatrix.cpp:184
SparseMatrix_ptr remote_coupleBlock
coupling of rows-cols on neighbouring processors (may not be valid)
Definition: SystemMatrix.h:315
boost::shared_ptr< Coupler > Coupler_ptr
Definition: Coupler.h:39
void solvePreconditioner(double *x, double *b)
Definition: SystemMatrix.cpp:155
index_t dim_t
Definition: types.h:27
void makeZeroRowSums(double *left_over)
Definition: SystemMatrix.cpp:198
void applyBalanceInPlace(double *x, bool RHS) const
Definition: SystemMatrix.cpp:371
#define MATRIX_FORMAT_CSC
Definition: Paso.h:54
dim_t getColOverlap() const
Definition: SystemMatrix.h:176
void mergeMainAndCouple(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:41
void copyBlockToMainDiagonal(const double *in)
Definition: SystemMatrix.h:223
Coupler_ptr row_coupler
Definition: SystemMatrix.h:306
void setPreconditioner(Options *options)
Definition: SystemMatrix.cpp:148
dim_t getGlobalTotalNumRows() const
Definition: SystemMatrix.h:197
boost::shared_ptr< JMPI_ > JMPI
Definition: Esys_MPI.h:79
index_t solver_package
package code controlling the solver pointer
Definition: SystemMatrix.h:330
Definition: error.h:46
double * finishColCollect()
Definition: SystemMatrix.h:136
void copyColCoupleBlock()
Definition: SystemMatrix.cpp:291
void saveMM(const char *filename) const
Definition: SystemMatrix.h:246
void copyMain_CSC_OFFSET1(index_t **p_ptr, index_t **p_idx, double **p_val)
esysUtils::JMPI mpi_info
Definition: SystemMatrix.h:303