escript  Revision_
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 
35 #include <escript/AbstractSystemMatrix.h>
36 
37 namespace paso {
38 
39 class Options;
41 typedef boost::shared_ptr<SystemMatrix> SystemMatrix_ptr;
42 typedef boost::shared_ptr<const SystemMatrix> const_SystemMatrix_ptr;
43 
44 typedef int SystemMatrixType;
45 
48 {
49 public:
51  SystemMatrix();
52 
54  dim_t rowBlockSize, dim_t columnBlockSize,
55  bool patternIsUnrolled, const escript::FunctionSpace& rowFS,
56  const escript::FunctionSpace& colFS);
57 
58  ~SystemMatrix();
59 
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.");
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 
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 
111  void copyColCoupleBlock();
112 
113  void copyRemoteCoupleBlock(bool recreatePattern);
114 
115  void fillWithGlobalCoordinates(double f1);
116 
117  void print() const;
118 
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 
148  void solvePreconditioner(double* x, double* b);
149 
150  void freePreconditioner();
151 
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 
231  {
232  return getGlobalNumRows() * row_block_size;
233  }
234 
236  {
237  return getGlobalNumCols() * col_block_size;
238  }
239 
240  inline double getSparsity() const
241  {
242  return getGlobalSize() /
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;
311 
314 
318 
322 
325 
334 
336 
342  double* balance_vector;
343 
345  mutable index_t* global_id;
346 
349 
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 
void startColCollect(const double *in) const
Definition: SystemMatrix.h:164
SparseMatrix_ptr col_coupleBlock
coupling to neighbouring processors (row - col)
Definition: SystemMatrix.h:329
Definition: FunctionSpace.h:34
static int getSystemMatrixTypeId(int solver, int preconditioner, int package, bool symmetry, const escript::JMPI &mpi_info)
Definition: SystemMatrix.cpp:520
dim_t getNumRows() const
Definition: SystemMatrix.h:184
void extendedRowsForST(dim_t *degree_ST, index_t *offset_ST, index_t *ST)
Definition: SystemMatrix_extendedRows.cpp:41
dim_t getTotalNumCols() const
Definition: SystemMatrix.h:199
double * balance_vector
Definition: SystemMatrix.h:342
void copyBlockToMainDiagonal(const double *in)
Definition: SystemMatrix.h:256
SparseMatrix_ptr row_coupleBlock
coupling to neighbouring processors (col - row)
Definition: SystemMatrix.h:331
dim_t getColOverlap() const
Definition: SystemMatrix.h:209
~SystemMatrix()
Definition: SystemMatrix.cpp:142
void nullifyRows(double *mask_row, double main_diagonal_value)
Definition: SystemMatrix.cpp:223
void print() const
Definition: SystemMatrix_debug.cpp:100
dim_t logical_col_block_size
Definition: SystemMatrix.h:313
Definition: Options.h:89
void solve(double *out, double *in, Options *options) const
Definition: solve.cpp:39
#define MATRIX_FORMAT_OFFSET1
Definition: Paso.h:63
void mergeMainAndCouple(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:40
boost::shared_ptr< Distribution > Distribution_ptr
Definition: Distribution.h:24
SystemMatrixPattern_ptr pattern
Definition: SystemMatrix.h:310
void copyMain_CSC_OFFSET1(index_t **p_ptr, index_t **p_idx, double **p_val)
dim_t getTotalNumRows() const
Definition: SystemMatrix.h:194
void mergeMainAndCouple_CSC_OFFSET1(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:276
void copyBlockFromMainDiagonal(double *out) const
Definition: SystemMatrix.h:251
void mergeMainAndCouple_CSR_OFFSET0_Block(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:168
dim_t block_size
Definition: SystemMatrix.h:317
boost::shared_ptr< JMPI_ > JMPI
Definition: EsysMPI.h:71
escript::Distribution_ptr col_distribution
Definition: SystemMatrix.h:320
dim_t getGlobalNumRows() const
Definition: SystemMatrix.h:214
void startCollect(const double *in) const
Definition: SystemMatrix.h:154
index_t * borrowMainDiagonalPointer() const
Definition: SystemMatrix.cpp:187
void applyBalanceInPlace(double *x, bool RHS) const
Definition: SystemMatrix.cpp:431
boost::shared_ptr< SparseMatrix > SparseMatrix_ptr
Definition: SparseMatrix.h:35
double * finishCollect() const
Definition: SystemMatrix.h:159
void add(dim_t, index_t *, dim_t, dim_t, index_t *, dim_t, double *)
void copyToMainDiagonal(const double *in)
Definition: SystemMatrix.h:266
dim_t getNumOutput() const
Definition: SystemMatrix.h:246
virtual void saveMM(const std::string &filename) const
writes the matrix to a file using the Matrix Market file format
Definition: SystemMatrix.h:68
void setValues(double value)
Definition: SystemMatrix.h:271
double getSparsity() const
Definition: SystemMatrix.h:240
virtual void saveHB(const std::string &filename) const
writes the matrix to a file using the Harwell-Boeing file format
Definition: SystemMatrix.h:80
escript::Distribution_ptr row_distribution
Definition: SystemMatrix.h:319
boost::shared_ptr< SystemMatrixPattern > SystemMatrixPattern_ptr
Definition: SystemMatrixPattern.h:39
void RHS_loadMM_toCSR(const char *filename, double *b, dim_t size)
Definition: SystemMatrix_loadMM.cpp:299
bool is_balanced
Definition: SystemMatrix.h:335
boost::shared_ptr< SystemMatrix > SystemMatrix_ptr
Definition: SystemMatrix.h:40
void copyFromMainDiagonal(double *out) const
Definition: SystemMatrix.h:261
SparseMatrix_ptr mainBlock
main block
Definition: SystemMatrix.h:327
dim_t row_block_size
Definition: SystemMatrix.h:315
void * solver_p
pointer to data needed by a solver
Definition: SystemMatrix.h:351
void freePreconditioner()
Definition: SystemMatrix.cpp:164
Definition: AMG.cpp:45
static SystemMatrix_ptr loadMM_toCSC(const char *filename)
Definition: SystemMatrix_loadMM.cpp:200
int index_t
type for array/matrix indices used both globally and on each rank
Definition: DataTypes.h:59
void copyRemoteCoupleBlock(bool recreatePattern)
Definition: SystemMatrix_copyRemoteCoupleBlock.cpp:39
virtual void setToSolution(escript::Data &out, escript::Data &in, boost::python::object &options) const
solves the linear system this*out=in
Definition: SystemMatrix.cpp:309
void mergeMainAndCouple_CSR_OFFSET0(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:56
void startRowCollect(const double *in)
Definition: SystemMatrix.h:174
void setPreconditioner(Options *options)
Definition: SystemMatrix.cpp:149
this class holds a (distributed) stiffness matrix
Definition: SystemMatrix.h:47
Data represents a collection of datapoints.
Definition: Data.h:63
void solvePreconditioner(double *x, double *b)
Definition: SystemMatrix.cpp:157
virtual void ypAx(escript::Data &y, escript::Data &x) const
performs y+=this*x
Definition: SystemMatrix.cpp:333
void balance()
Definition: SystemMatrix.cpp:469
double getGlobalSize() const
Definition: SystemMatrix.cpp:171
boost::shared_ptr< const SystemMatrix > const_SystemMatrix_ptr
Definition: SystemMatrix.h:42
index_t * global_id
stores the global ids for all cols in col_coupleBlock
Definition: SystemMatrix.h:345
int SystemMatrixType
Definition: SystemMatrix.h:44
SparseMatrix_ptr remote_coupleBlock
coupling of rows-cols on neighbouring processors (may not be valid)
Definition: SystemMatrix.h:333
dim_t logical_row_block_size
Definition: SystemMatrix.h:312
Coupler_ptr row_coupler
Definition: SystemMatrix.h:324
void fillWithGlobalCoordinates(double f1)
Definition: SystemMatrix_debug.cpp:37
virtual void nullifyRowsAndCols(escript::Data &mask_row, escript::Data &mask_col, double main_diagonal_value)
Definition: SystemMatrix.cpp:244
void applyBalance(double *x_out, const double *x, bool RHS) const
Definition: SystemMatrix.cpp:450
dim_t getNumCols() const
Definition: SystemMatrix.h:189
SystemMatrixType type
Definition: SystemMatrix.h:309
boost::shared_ptr< Coupler > Coupler_ptr
Definition: Coupler.h:41
dim_t getGlobalNumCols() const
Definition: SystemMatrix.h:222
dim_t getGlobalTotalNumRows() const
Definition: SystemMatrix.h:230
Base class for escript system matrices.
Definition: AbstractSystemMatrix.h:42
index_t solver_package
package code controlling the solver pointer
Definition: SystemMatrix.h:348
#define MATRIX_FORMAT_CSC
Definition: Paso.h:61
escript::JMPI mpi_info
Definition: SystemMatrix.h:321
PasoException exception class.
Definition: PasoException.h:32
dim_t getRowOverlap() const
Definition: SystemMatrix.h:204
void MatrixVector(double alpha, const double *in, double beta, double *out) const
Definition: SystemMatrix_MatrixVector.cpp:34
void makeZeroRowSums(double *left_over)
Definition: SystemMatrix.cpp:201
double * finishColCollect() const
Definition: SystemMatrix.h:169
void rowSum(double *row_sum) const
Definition: SystemMatrix.h:279
dim_t col_block_size
Definition: SystemMatrix.h:316
Coupler_ptr col_coupler
Definition: SystemMatrix.h:323
double * finishRowCollect()
Definition: SystemMatrix.h:179
SystemMatrix()
default constructor - throws exception.
Definition: SystemMatrix.cpp:41
void MatrixVector_CSR_OFFSET0(double alpha, const double *in, double beta, double *out) const
Definition: SystemMatrix_MatrixVector.cpp:63
virtual void resetValues(bool preserveSolverData=false)
resets the matrix entries
Definition: SystemMatrix.cpp:302
index_t dim_t
Definition: DataTypes.h:64
void copyColCoupleBlock()
Definition: SystemMatrix.cpp:353
static SystemMatrix_ptr loadMM_toCSR(const char *filename)
Definition: SystemMatrix_loadMM.cpp:100
SparseMatrix_ptr mergeSystemMatrix() const
Definition: SystemMatrix.cpp:554
dim_t getGlobalTotalNumCols() const
Definition: SystemMatrix.h:235