Go to the documentation of this file.
29 #ifndef __PASO_SPARSEMATRIX_H__
30 #define __PASO_SPARSEMATRIX_H__
36 template <
typename T>
struct SparseMatrix;
44 struct SparseMatrix : boost::enable_shared_from_this<SparseMatrix<T> >
48 bool patternIsUnrolled);
63 const double* b)
const;
71 const index_t* new_col_index)
const;
79 void saveMM(
const char* filename)
const;
83 return pattern->borrowMainDiagonalPointer();
88 return pattern->borrowColoringPointer();
135 const double* mask_col,
136 double main_diagonal_value);
139 const double* mask_col,
140 double main_diagonal_value);
143 double main_diagonal_value);
146 double main_diagonal_value);
149 double main_diagonal_value);
185 const double beta,
double* out);
190 const double beta,
double* out);
195 const double beta,
double* out);
197 template <
typename T>
201 const double beta, T* out);
206 const double beta,
double* out);
224 struct Preconditioner_LocalSmoother;
243 template <
typename T>
246 bool patternIsUnrolled) :
252 if (patternIsUnrolled) {
254 throw PasoException(
"SparseMatrix: requested offset and pattern offset do not match.");
260 = (rowBlockSize != colBlockSize)
261 #ifndef ESYS_HAVE_LAPACK
263 || (colBlockSize > 3)
276 if (patternIsUnrolled) {
279 pattern = npattern->unrollBlocks(pattern_format_out,
280 colBlockSize, rowBlockSize);
285 pattern = npattern->unrollBlocks(pattern_format_out, 1, 1);
294 if (patternIsUnrolled) {
297 pattern = npattern->unrollBlocks(pattern_format_out,
298 rowBlockSize, colBlockSize);
303 pattern = npattern->unrollBlocks(pattern_format_out, 1, 1);
321 template <
typename T>
324 switch (solver_package) {
344 template <
typename T>
348 if (!pattern->isEmpty()) {
349 const dim_t nOut = pattern->numOutput;
350 #pragma omp parallel for
351 for (
dim_t i=0; i < nOut; ++i) {
352 for (
index_t iptr=pattern->ptr[i]-index_offset; iptr < pattern->ptr[i+1]-index_offset; ++iptr) {
353 for (
dim_t j=0; j<block_size; ++j)
354 val[iptr*block_size+j] = value;
361 template <
typename T>
367 const int totalRowSize = A->numRows * A->row_block_size;
368 if (std::abs(beta) > 0) {
370 #pragma omp parallel for schedule(static)
371 for (
index_t irow=0; irow < totalRowSize; irow++) {
376 #pragma omp parallel for schedule(static)
377 for (
index_t irow=0; irow < totalRowSize; irow++) {
382 if (std::abs(alpha) > 0) {
383 const int nRows = A->pattern->numOutput;
384 if (A->col_block_size==1 && A->row_block_size==1) {
385 #pragma omp parallel for schedule(static)
386 for (
index_t irow=0; irow < nRows; irow++) {
389 for (
index_t iptr=A->pattern->ptr[irow]-1;
390 iptr < A->pattern->ptr[irow+1]-1; ++iptr) {
391 reg += A->val[iptr] * in[A->pattern->index[iptr]-1];
393 out[irow] += alpha * reg;
395 }
else if (A->col_block_size==2 && A->row_block_size==2) {
396 #pragma omp parallel for schedule(static)
397 for (
index_t ir=0; ir < nRows; ir++) {
401 for (
index_t iptr=A->pattern->ptr[ir]-1;
402 iptr < A->pattern->ptr[ir+1]-1; iptr++) {
403 const index_t ic=2*(A->pattern->index[iptr]-1);
404 reg1 += A->val[iptr*4 ]*in[ic] + A->val[iptr*4+2]*in[1+ic];
405 reg2 += A->val[iptr*4+1]*in[ic] + A->val[iptr*4+3]*in[1+ic];
407 out[ 2*ir] += alpha * reg1;
408 out[1+2*ir] += alpha * reg2;
410 }
else if (A->col_block_size==3 && A->row_block_size==3) {
411 #pragma omp parallel for schedule(static)
412 for (
index_t ir=0; ir < nRows; ir++) {
417 for (
index_t iptr=A->pattern->ptr[ir]-1;
418 iptr < A->pattern->ptr[ir+1]-1; iptr++) {
419 const index_t ic=3*(A->pattern->index[iptr]-1);
420 reg1 += A->val[iptr*9 ]*in[ic] + A->val[iptr*9+3]*in[1+ic] + A->val[iptr*9+6]*in[2+ic];
421 reg2 += A->val[iptr*9+1]*in[ic] + A->val[iptr*9+4]*in[1+ic] + A->val[iptr*9+7]*in[2+ic];
422 reg3 += A->val[iptr*9+2]*in[ic] + A->val[iptr*9+5]*in[1+ic] + A->val[iptr*9+8]*in[2+ic];
424 out[ 3*ir] += alpha * reg1;
425 out[1+3*ir] += alpha * reg2;
426 out[2+3*ir] += alpha * reg3;
429 #pragma omp parallel for schedule(static)
430 for (
index_t ir=0; ir < nRows; ir++) {
431 for (
index_t iptr=A->pattern->ptr[ir]-1;
432 iptr < A->pattern->ptr[ir+1]-1; iptr++) {
433 for (
index_t irb=0; irb < A->row_block_size; irb++) {
436 for (
index_t icb=0; icb < A->col_block_size; icb++) {
437 const index_t icol=icb+A->col_block_size*(A->pattern->index[iptr]-1);
438 reg += A->val[iptr*A->block_size+irb+A->row_block_size*icb] * in[icol];
440 const index_t irow=irb+A->row_block_size*ir;
441 out[irow] += alpha * reg;
449 template <
typename T>
451 const double* mask_col,
452 double main_diagonal_value)
455 const int nOut = pattern->numOutput;
456 #pragma omp parallel for
457 for (
index_t icol=0; icol < nOut; icol++) {
459 for (
index_t iptr=pattern->ptr[icol]-index_offset; iptr < pattern->ptr[icol+1]-index_offset; iptr++) {
460 const index_t irow = pattern->index[iptr]-index_offset;
461 if (mask_col[icol]>0. || mask_row[irow]>0.) {
462 val[iptr] = (irow==icol ? main_diagonal_value : 0);
468 template <
typename T>
470 const double* mask_col,
471 double main_diagonal_value)
474 const int nOut = pattern->numOutput;
475 #pragma omp parallel for
476 for (
index_t ic=0; ic < nOut; ic++) {
477 for (
index_t iptr=pattern->ptr[ic]-index_offset; iptr < pattern->ptr[ic+1]-index_offset; iptr++) {
478 for (
index_t irb=0; irb < row_block_size; irb++) {
479 const index_t irow=irb+row_block_size*(pattern->index[iptr]-index_offset);
481 for (
index_t icb=0; icb < col_block_size; icb++) {
482 const index_t icol=icb+col_block_size*ic;
483 if (mask_col[icol]>0. || mask_row[irow]>0.) {
484 const index_t l=iptr*block_size+irb+row_block_size*icb;
485 val[l] = (irow==icol ? main_diagonal_value : 0);
493 template <
typename T>
495 const double* mask_col,
496 double main_diagonal_value)
499 const int nOut = pattern->numOutput;
500 #pragma omp parallel for
501 for (
index_t irow=0; irow < nOut; irow++) {
503 for (
index_t iptr=pattern->ptr[irow]-index_offset; iptr < pattern->ptr[irow+1]-index_offset; iptr++) {
504 const index_t icol = pattern->index[iptr]-index_offset;
505 if (mask_col[icol]>0. || mask_row[irow]>0.) {
506 val[iptr] = (irow==icol ? main_diagonal_value : 0);
512 template <
typename T>
514 const double* mask_col,
515 double main_diagonal_value)
518 const int nOut = pattern->numOutput;
519 #pragma omp parallel for
520 for (
index_t ir=0; ir < nOut; ir++) {
521 for (
index_t iptr=pattern->ptr[ir]-index_offset; iptr < pattern->ptr[ir+1]-index_offset; iptr++) {
522 for (
index_t irb=0; irb < row_block_size; irb++) {
523 const index_t irow=irb+row_block_size*ir;
525 for (
index_t icb=0; icb < col_block_size; icb++) {
526 const index_t icol=icb+col_block_size*(pattern->index[iptr]-index_offset);
527 if (mask_col[icol]>0. || mask_row[irow]>0.) {
528 const index_t l=iptr*block_size+irb+row_block_size*icb;
529 val[l] = (irow==icol ? main_diagonal_value : 0);
539 #endif // __PASO_SPARSEMATRIX_H__
#define MATRIX_FORMAT_BLK1
Definition: Paso.h:63
int mm_read_mtx_crd_size(std::istream &f, int *M, int *N, int *nz)
Definition: mmio.cpp:181
static dim_t N
Definition: SparseMatrix_saveHB.cpp:37
#define PASO_DLL_API
Definition: paso/src/system_dep.h:29
int mm_write_mtx_crd_size(std::ostream &f, int M, int N, int nz)
Definition: mmio.cpp:173
boost::shared_ptr< SparseMatrix< T > > SparseMatrix_ptr
Definition: SparseMatrix.h:37
void SparseMatrix_MatrixVector_CSC_OFFSET1(const double alpha, const_SparseMatrix_ptr< double > A, const double *in, const double beta, double *out)
Definition: SparseMatrix_MatrixVector.cpp:119
void nullifyRows_CSR_BLK1(const double *mask_row, double main_diagonal_value)
static Pattern_ptr fromIndexListArray(dim_t n0, dim_t n, const escript::IndexList *index_list_array, index_t range_min, index_t range_max, index_t index_offset)
Definition: Pattern.cpp:105
dim_t getNumCols() const
Definition: SparseMatrix.h:116
SparseMatrixType type
Definition: SparseMatrix.h:161
dim_t getNumRows() const
Definition: SparseMatrix.h:111
#define MATRIX_FORMAT_CSC
Definition: Paso.h:62
void applyDiagonal_CSR_OFFSET0(const double *left, const double *right)
SparseMatrix_ptr< double > getSubmatrix(dim_t n_row_sub, dim_t n_col_sub, const index_t *row_list, const index_t *new_col_index) const
dim_t len
Definition: SparseMatrix.h:168
index_t * borrowColoringPointer() const
Definition: SparseMatrix.h:86
void applyBlockMatrix(double *block_diag, index_t *pivot, double *x, const double *b) const
#define mm_set_real(typecode)
Definition: mmio.h:79
void BlockOps_invM_2(double *invA, const double *A, int *failed)
Definition: BlockOps.h:99
void MUMPS_free(SparseMatrix< T > *A)
frees any MUMPS related data from the matrix
Definition: MUMPS.h:118
SparseMatrix_ptr< double > getBlock(int blockid) const
dim_t maxDeg() const
Definition: SparseMatrix.h:96
dim_t col_block_size
Definition: SparseMatrix.h:163
~SparseMatrix()
Definition: SparseMatrix.h:322
#define mm_is_general(typecode)
Definition: mmio.h:63
void maxAbsRow_CSR_OFFSET0(double *array) const
void saveMM(const char *filename) const
dim_t getTotalNumCols() const
Definition: SparseMatrix.h:106
void BlockOps_invM_N(dim_t N, double *mat, index_t *pivot, int *failed)
LU factorization of NxN matrix mat with partial pivoting.
Definition: BlockOps.h:148
void SparseMatrix_MatrixVector_CSC_OFFSET0(const double alpha, const_SparseMatrix_ptr< double > A, const double *in, const double beta, double *out)
Definition: SparseMatrix_MatrixVector.cpp:43
void addAbsRow_CSR_OFFSET0(double *array) const
#define PASO_MKL
Definition: Options.h:50
#define MATRIX_FORMAT_OFFSET1
Definition: Paso.h:64
void nullifyRowsAndCols_CSC_BLK1(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition: SparseMatrix.h:450
SparseMatrix_ptr< double > unroll(SparseMatrixType type) const
void BlockOps_solveAll(dim_t n_block, dim_t n, double *D, index_t *pivot, double *x)
Definition: BlockOps.h:219
void copy(dim_t N, double *out, const double *in)
out = in
Definition: PasoUtil.h:88
void Preconditioner_LocalSmoother_free(Preconditioner_LocalSmoother *in)
Definition: Smoother.cpp:42
static dim_t nz
Definition: SparseMatrix_saveHB.cpp:37
boost::shared_ptr< const SparseMatrix< T > > const_SparseMatrix_ptr
Definition: SparseMatrix.h:38
int comparIndex(const void *index1, const void *index2)
this int-comparison function is used by qsort/bsearch in various places
Definition: PasoUtil.cpp:25
index_t dim_t
Definition: DataTypes.h:66
void SparseMatrix_MatrixVector_CSR_OFFSET0_DIAG(const double alpha, const_SparseMatrix_ptr< double > A, const double *in, const double beta, double *out)
Definition: SparseMatrix_MatrixVector.cpp:338
dim_t block_size
Definition: SparseMatrix.h:164
double getSize() const
Definition: SparseMatrix.h:121
void BlockOps_invM_3(double *invA, const double *A, int *failed)
Definition: BlockOps.h:117
void copyBlockFromMainDiagonal(double *out) const
void saveHB_CSC(const char *filename) const
#define mm_set_matrix(typecode)
Definition: mmio.h:72
char MM_typecode[4]
Definition: mmio.h:35
double getSparsity() const
Definition: SparseMatrix.h:126
void swap(index_t *r, index_t *c, double *v, int left, int right)
Definition: SparseMatrix.cpp:57
dim_t numCols
Definition: SparseMatrix.h:166
index_t solver_package
package controlling the solver pointer
Definition: SparseMatrix.h:174
dim_t getTotalNumRows() const
Definition: SparseMatrix.h:101
void MKL_free(SparseMatrix< double > *A)
Definition: MKL.cpp:35
int SparseMatrixType
Definition: SparseMatrix.h:40
#define PASO_PASO
Definition: Options.h:56
#define PASO_UMFPACK
Definition: Options.h:51
Pattern_ptr pattern
Definition: SparseMatrix.h:167
void nullifyRowsAndCols_CSR_BLK1(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition: SparseMatrix.h:494
void copyFromMainDiagonal(double *out) const
dim_t getNumColors() const
Definition: SparseMatrix.h:91
void BlockOps_Cpy_N(dim_t N, double *R, const double *V)
Definition: BlockOps.h:40
dim_t numRows
Definition: SparseMatrix.h:165
PasoException exception class.
Definition: PasoException.h:34
dim_t row_block_size
Definition: SparseMatrix.h:162
SparseMatrix(SparseMatrixType type, Pattern_ptr pattern, dim_t rowBlockSize, dim_t colBlockSize, bool patternIsUnrolled)
Definition: SparseMatrix.h:244
#define mm_initialize_typecode(typecode)
Definition: mmio.h:92
int index_t
type for array/matrix indices used both globally and on each rank
Definition: DataTypes.h:61
#define mm_is_real(typecode)
Definition: mmio.h:58
void nullifyRowsAndCols_CSC(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition: SparseMatrix.h:469
void copyBlockToMainDiagonal(const double *in)
void UMFPACK_free(SparseMatrix< double > *A)
frees any UMFPACK related data from the matrix
Definition: UMFPACK.cpp:35
static SparseMatrix_ptr< double > loadMM_toCSR(const char *filename)
int mm_read_banner(std::istream &f, MM_typecode *matcode)
Definition: mmio.cpp:103
void SparseMatrix_MatrixVector_CSR_OFFSET0(const double alpha, const_SparseMatrix_ptr< double > A, const double *in, const double beta, double *out)
Definition: SparseMatrix_MatrixVector.cpp:191
SparseMatrix_ptr< double > getTranspose() const
#define PASO_MUMPS
Definition: Options.h:57
index_t * borrowMainDiagonalPointer() const
Definition: SparseMatrix.h:81
void * solver_p
pointer to data needed by a solver
Definition: SparseMatrix.h:177
T * val
this is used for classical CSR or CSC
Definition: SparseMatrix.h:171
void nullifyRows_CSR(const double *mask_row, double main_diagonal_value)
int mm_write_banner(std::ostream &f, MM_typecode matcode)
Definition: mmio.cpp:356
static dim_t M
Definition: SparseMatrix_saveHB.cpp:37
void invMain(double *inv_diag, index_t *pivot) const
boost::shared_ptr< Pattern > Pattern_ptr
Definition: Pattern.h:39
SparseMatrix_ptr< double > SparseMatrix_MatrixMatrix(const_SparseMatrix_ptr< double > A, const_SparseMatrix_ptr< double > B)
Definition: SparseMatrix_MatrixMatrix.cpp:44
#define mm_set_coordinate(typecode)
Definition: mmio.h:73
void nullifyRowsAndCols_CSR(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition: SparseMatrix.h:513
Definition: escriptcore/src/IndexList.h:29
SparseMatrix_ptr< double > SparseMatrix_MatrixMatrixTranspose(const_SparseMatrix_ptr< double > A, const_SparseMatrix_ptr< double > B, const_SparseMatrix_ptr< double > T)
Definition: SparseMatrix_MatrixMatrixTranspose.cpp:52
Definition: BiCGStab.cpp:25
#define PASO_SMOOTHER
Definition: Options.h:75
void SparseMatrix_MatrixVector_CSR_OFFSET1(const double alpha, const_SparseMatrix_ptr< T > A, const T *in, const double beta, T *out)
Definition: SparseMatrix.h:362
void setValues(T value)
Definition: SparseMatrix.h:345
Definition: SparseMatrix.h:45
void q_sort(index_t *row, index_t *col, double *val, int begin, int end, int N)
Definition: SparseMatrix.cpp:75
Definition: Preconditioner.h:57
#define MATRIX_FORMAT_DEFAULT
Definition: Paso.h:61
#define MATRIX_FORMAT_DIAGONAL_BLOCK
Definition: Paso.h:65
void addRow_CSR_OFFSET0(double *array) const
#define mm_is_sparse(typecode)
Definition: mmio.h:52
void copyToMainDiagonal(const double *in)