30#ifndef __PASO_SYSTEMMATRIX_H__
31#define __PASO_SYSTEMMATRIX_H__
36#include <escript/AbstractSystemMatrix.h>
41template <
class T>
class SystemMatrix;
68 double main_diagonal_value);
70 virtual inline void saveMM(
const std::string& filename)
const
76 merged->saveMM(filename.c_str());
82 virtual inline void saveHB(
const std::string& filename)
const
85 throw PasoException(
"SystemMatrix::saveHB: Only single rank supported.");
87 throw PasoException(
"SystemMatrix::saveHB: Only CSC format supported.");
93 virtual void resetValues(
bool preserveSolverData =
false);
99 void nullifyRows(
double* mask_row,
double main_diagonal_value);
219 return pattern->input_distribution->getGlobalNumComponents();
221 return pattern->output_distribution->getGlobalNumComponents();
227 return pattern->output_distribution->getGlobalNumComponents();
229 return pattern->input_distribution->getGlobalNumComponents();
250 return pattern->getNumOutput();
255 mainBlock->copyBlockFromMainDiagonal(out);
281 inline void rowSum(
double* row_sum)
const
284 throw PasoException(
"SystemMatrix::rowSum: No normalization "
285 "available for compressed sparse column or index offset 1.");
288#pragma omp parallel for
289 for (
index_t irow=0; irow<nrow; ++irow) {
308 int package,
bool is_complex,
bool symmetry,
357 boost::python::object& options)
const;
373#include <escript/Data.h>
395 throw PasoException(
"SystemMatrix: Illegal to generate default SystemMatrix.");
406 dim_t colBlockSize,
bool patternIsUnrolled,
409 escript::AbstractSystemMatrix(rowBlockSize, rowFS, colBlockSize, colFS),
411 logical_row_block_size(rowBlockSize),
412 logical_col_block_size(colBlockSize),
414 balance_vector(NULL),
419 if (patternIsUnrolled) {
421 throw PasoException(
"SystemMatrix: requested offset and pattern offset do not match.");
427 = (rowBlockSize != colBlockSize)
428#ifndef ESYS_HAVE_LAPACK
430 || (colBlockSize > 3)
444 if (patternIsUnrolled) {
447 pattern = npattern->unrollBlocks(pattern_format_out,
448 colBlockSize, rowBlockSize);
453 pattern = npattern->unrollBlocks(pattern_format_out, 1, 1);
461 if (patternIsUnrolled) {
464 pattern = npattern->unrollBlocks(pattern_format_out,
465 rowBlockSize, colBlockSize);
470 pattern = npattern->unrollBlocks(pattern_format_out, 1, 1);
489#pragma omp parallel for
490 for (
dim_t i=0; i<n_norm; ++i)
499 delete[] balance_vector;
505 int package,
bool is_complex,
bool symmetry,
513 switch(true_package) {
523 if (mpi_info->size > 1) {
525 "requires CSC format which is not supported with "
526 "more than one rank.");
539 if (out > 0 && is_complex)
547 double main_diagonal_value)
551 throw PasoException(
"SystemMatrix::nullifyRowsAndCols: complex arguments not supported");
554 throw PasoException(
"nullifyRowsAndCols: column block size does not match the number of components of column mask.");
556 throw PasoException(
"nullifyRowsAndCols: row block size does not match the number of components of row mask.");
558 throw PasoException(
"nullifyRowsAndCols: column function space and function space of column mask don't match.");
560 throw PasoException(
"nullifyRowsAndCols: row function space and function space of row mask don't match.");
569 if (mpi_info->size > 1) {
572 "CSC is not supported with MPI.");
575 startColCollect(mask_col);
576 startRowCollect(mask_row);
577 if (col_block_size==1 && row_block_size==1) {
578 mainBlock->nullifyRowsAndCols_CSR_BLK1(mask_row, mask_col, main_diagonal_value);
579 double* remote_values = finishColCollect();
580 col_coupleBlock->nullifyRowsAndCols_CSR_BLK1(mask_row, remote_values, 0.);
581 remote_values = finishRowCollect();
582 row_coupleBlock->nullifyRowsAndCols_CSR_BLK1(remote_values, mask_col, 0.);
584 mainBlock->nullifyRowsAndCols_CSR(mask_row, mask_col, main_diagonal_value);
585 double* remote_values = finishColCollect();
586 col_coupleBlock->nullifyRowsAndCols_CSR(mask_row, remote_values, 0.);
587 remote_values = finishRowCollect();
588 row_coupleBlock->nullifyRowsAndCols_CSR(remote_values, mask_col, 0.);
591 if (col_block_size==1 && row_block_size==1) {
593 mainBlock->nullifyRowsAndCols_CSC_BLK1(mask_row, mask_col, main_diagonal_value);
595 mainBlock->nullifyRowsAndCols_CSR_BLK1(mask_row, mask_col, main_diagonal_value);
599 mainBlock->nullifyRowsAndCols_CSC(mask_row, mask_col, main_diagonal_value);
601 mainBlock->nullifyRowsAndCols_CSR(mask_row, mask_col, main_diagonal_value);
611 if (!preserveSolverData)
617 boost::python::object& options)
const
619#if !defined(ESYS_HAVE_MUMPS)
622 throw PasoException(
"SystemMatrix::setToSolution: complex arguments not supported.");
625 options.attr(
"resetDiagnostics")();
628 throw PasoException(
"solve: column block size does not match the number of components of solution.");
630 throw PasoException(
"solve: row block size does not match the number of components of right hand side.");
632 throw PasoException(
"solve: column function space and function space of solution don't match.");
634 throw PasoException(
"solve: row function space and function space of right hand side don't match.");
642 solve(out_dp, in_dp, &paso_options);
649#if !defined(ESYS_HAVE_MUMPS)
652 throw PasoException(
"SystemMatrix::ypAx: complex arguments not supported.");
656 throw PasoException(
"matrix vector product: column block size does not match the number of components in input.");
658 throw PasoException(
"matrix vector product: row block size does not match the number of components in output.");
660 throw PasoException(
"matrix vector product: column function space and function space of input don't match.");
662 throw PasoException(
"matrix vector product: row function space and function space of output don't match.");
670 MatrixVector(1., x_dp, 1., y_dp);
#define PASO_MUMPS
Definition: Options.h:57
#define PASO_UMFPACK
Definition: Options.h:51
#define PASO_PASO
Definition: Options.h:56
#define PASO_MKL
Definition: Options.h:50
#define MATRIX_FORMAT_BLK1
Definition: Paso.h:63
#define MATRIX_FORMAT_DEFAULT
Definition: Paso.h:61
#define MATRIX_FORMAT_OFFSET1
Definition: Paso.h:64
#define MATRIX_FORMAT_CSC
Definition: Paso.h:62
#define MATRIX_FORMAT_DIAGONAL_BLOCK
Definition: Paso.h:65
#define MATRIX_FORMAT_COMPLEX
Definition: Paso.h:66
Base class for escript system matrices.
Definition: AbstractSystemMatrix.h:44
ElementType * data()
Definition: DataVectorAlt.h:200
Data represents a collection of datapoints.
Definition: Data.h:64
const FunctionSpace & getFunctionSpace() const
Return the function space.
Definition: Data.h:463
int getDataPointSize() const
Return the size of the data point. It is the product of the data point shape dimensions.
Definition: Data.cpp:1360
DataTypes::RealVectorType & getExpandedVectorReference(DataTypes::real_t dummy=0)
Ensures that the Data is expanded and returns its underlying vector Does not check for exclusive writ...
Definition: Data.cpp:5837
void expand()
Whatever the current Data type make this into a DataExpanded.
Definition: Data.cpp:1180
bool isComplex() const
True if components of this data are stored as complex.
Definition: Data.cpp:1160
void requireWrite()
Ensures data is ready for write access. This means that the data will be resolved if lazy and will be...
Definition: Data.cpp:1239
Definition: FunctionSpace.h:36
PasoException exception class.
Definition: PasoException.h:34
this class holds a (distributed) stiffness matrix
Definition: SystemMatrix.h:50
SparseMatrix_ptr< T > col_coupleBlock
coupling to neighbouring processors (row - col)
Definition: SystemMatrix.h:331
void * solver_p
pointer to data needed by a solver
Definition: SystemMatrix.h:353
dim_t getGlobalTotalNumRows() const
Definition: SystemMatrix.h:232
SystemMatrixType type
Definition: SystemMatrix.h:311
SparseMatrix_ptr< T > row_coupleBlock
coupling to neighbouring processors (col - row)
Definition: SystemMatrix.h:333
dim_t getGlobalNumCols() const
Definition: SystemMatrix.h:224
escript::Distribution_ptr col_distribution
Definition: SystemMatrix.h:322
void copyBlockToMainDiagonal(const double *in)
Definition: SystemMatrix.h:258
static SystemMatrix_ptr< double > loadMM_toCSR(const char *filename)
dim_t getColOverlap() const
Definition: SystemMatrix.h:211
virtual void saveHB(const std::string &filename) const
writes the matrix to a file using the Harwell-Boeing file format
Definition: SystemMatrix.h:82
void applyBalance(double *x_out, const double *x, bool RHS) const
void mergeMainAndCouple_CSR_OFFSET0(index_t **p_ptr, index_t **p_idx, double **p_val) const
static int getSystemMatrixTypeId(int solver, int preconditioner, int package, bool is_complex, bool symmetry, const escript::JMPI &mpi_info)
Definition: SystemMatrix.h:504
void copyBlockFromMainDiagonal(double *out) const
Definition: SystemMatrix.h:253
void add(dim_t, index_t *, dim_t, dim_t, index_t *, dim_t, double *)
SparseMatrix_ptr< T > mergeSystemMatrix() const
void mergeMainAndCouple_CSR_OFFSET0_Block(index_t **p_ptr, index_t **p_idx, double **p_val) const
void solve(T *out, T *in, Options *options) const
void MatrixVector_CSR_OFFSET0(double alpha, const double *in, double beta, double *out) const
SystemMatrixPattern_ptr pattern
Definition: SystemMatrix.h:312
escript::Distribution_ptr row_distribution
Definition: SystemMatrix.h:321
static SystemMatrix_ptr< double > loadMM_toCSC(const char *filename)
void extendedRowsForST(dim_t *degree_ST, index_t *offset_ST, index_t *ST)
dim_t getNumCols() const
Definition: SystemMatrix.h:191
escript::JMPI mpi_info
Definition: SystemMatrix.h:323
SparseMatrix_ptr< T > remote_coupleBlock
coupling of rows-cols on neighbouring processors (may not be valid)
Definition: SystemMatrix.h:335
void copyColCoupleBlock()
void applyBalanceInPlace(double *x, bool RHS) const
void startColCollect(const double *in) const
Definition: SystemMatrix.h:166
void MatrixVector(double alpha, const T *in, double beta, T *out) const
void fillWithGlobalCoordinates(double f1)
dim_t col_block_size
Definition: SystemMatrix.h:318
index_t * borrowMainDiagonalPointer() const
index_t solver_package
package code controlling the solver pointer
Definition: SystemMatrix.h:350
void nullifyRows(double *mask_row, double main_diagonal_value)
dim_t getGlobalNumRows() const
Definition: SystemMatrix.h:216
dim_t getRowOverlap() const
Definition: SystemMatrix.h:206
void startCollect(const double *in) const
Definition: SystemMatrix.h:156
dim_t logical_row_block_size
Definition: SystemMatrix.h:314
void copyFromMainDiagonal(double *out) const
Definition: SystemMatrix.h:263
dim_t getGlobalTotalNumCols() const
Definition: SystemMatrix.h:237
SystemMatrix()
default constructor - throws exception.
Definition: SystemMatrix.h:393
void startRowCollect(const double *in)
Definition: SystemMatrix.h:176
dim_t getTotalNumCols() const
Definition: SystemMatrix.h:201
dim_t getNumRows() const
Definition: SystemMatrix.h:186
SparseMatrix_ptr< T > mainBlock
main block
Definition: SystemMatrix.h:329
void copyToMainDiagonal(const double *in)
Definition: SystemMatrix.h:268
double getGlobalSize() const
void makeZeroRowSums(double *left_over)
dim_t getNumOutput() const
Definition: SystemMatrix.h:248
void mergeMainAndCouple_CSC_OFFSET1(index_t **p_ptr, index_t **p_idx, double **p_val) const
void copyMain_CSC_OFFSET1(index_t **p_ptr, index_t **p_idx, double **p_val)
virtual void ypAx(escript::Data &y, escript::Data &x) const
performs y+=this*x
Definition: SystemMatrix.h:647
index_t * global_id
stores the global ids for all cols in col_coupleBlock
Definition: SystemMatrix.h:347
double * balance_vector
Definition: SystemMatrix.h:344
virtual void resetValues(bool preserveSolverData=false)
resets the matrix entries
Definition: SystemMatrix.h:608
~SystemMatrix()
Definition: SystemMatrix.h:496
void copyRemoteCoupleBlock(bool recreatePattern)
dim_t row_block_size
Definition: SystemMatrix.h:317
void setPreconditioner(Options *options)
dim_t block_size
Definition: SystemMatrix.h:319
dim_t getTotalNumRows() const
Definition: SystemMatrix.h:196
double * finishColCollect() const
Definition: SystemMatrix.h:171
Coupler_ptr< real_t > col_coupler
Definition: SystemMatrix.h:325
virtual void setToSolution(escript::Data &out, escript::Data &in, boost::python::object &options) const
solves the linear system this*out=in
Definition: SystemMatrix.h:616
virtual void nullifyRowsAndCols(escript::Data &mask_row, escript::Data &mask_col, double main_diagonal_value)
Definition: SystemMatrix.h:545
double * finishRowCollect()
Definition: SystemMatrix.h:181
dim_t logical_col_block_size
Definition: SystemMatrix.h:315
bool is_balanced
Definition: SystemMatrix.h:337
double * finishCollect() const
Definition: SystemMatrix.h:161
Coupler_ptr< real_t > row_coupler
Definition: SystemMatrix.h:326
void rowSum(double *row_sum) const
Definition: SystemMatrix.h:281
virtual void saveMM(const std::string &filename) const
writes the matrix to a file using the Matrix Market file format
Definition: SystemMatrix.h:70
void solvePreconditioner(double *x, double *b)
void freePreconditioner()
void mergeMainAndCouple(index_t **p_ptr, index_t **p_idx, double **p_val) const
void setValues(double value)
Definition: SystemMatrix.h:273
double getSparsity() const
Definition: SystemMatrix.h:242
std::complex< real_t > cplx_t
complex data type
Definition: DataTypes.h:55
double real_t
type of all real-valued scalars in escript
Definition: DataTypes.h:52
index_t dim_t
Definition: DataTypes.h:66
int index_t
type for array/matrix indices used both globally and on each rank
Definition: DataTypes.h:61
Definition: AbstractContinuousDomain.cpp:23
boost::shared_ptr< Distribution > Distribution_ptr
Definition: Distribution.h:26
boost::shared_ptr< JMPI_ > JMPI
Definition: EsysMPI.h:76
Definition: BiCGStab.cpp:25
boost::shared_ptr< SystemMatrix< T > > SystemMatrix_ptr
Definition: SystemMatrix.h:42
boost::shared_ptr< SparseMatrix< T > > SparseMatrix_ptr
Definition: SparseMatrix.h:37
boost::shared_ptr< SystemMatrixPattern > SystemMatrixPattern_ptr
Definition: SystemMatrixPattern.h:41
boost::shared_ptr< const SystemMatrix< T > > const_SystemMatrix_ptr
Definition: SystemMatrix.h:43
void RHS_loadMM_toCSR(const char *filename, double *b, dim_t size)
Definition: SystemMatrix_loadMM.cpp:301
void solve_free(SystemMatrix< T > *A)
Definition: Solver.h:79
int SystemMatrixType
Definition: SystemMatrix.h:45
boost::shared_ptr< Coupler< T > > Coupler_ptr
Definition: Coupler.h:43
#define PASO_DLL_API
Definition: paso/src/system_dep.h:29
Definition: Coupler.h:100
static int getPackage(int solver, int package, bool symmetry, const escript::JMPI &mpi_info)
Definition: Options.cpp:315
void updateEscriptDiagnostics(boost::python::object &options) const
updates SolverBuddy diagnostics from this
Definition: Options.cpp:441
static int mapEscriptOption(int escriptOption)
returns the corresponding paso option code for an escript option code
Definition: Options.cpp:360
Definition: SparseMatrix.h:45