29#ifndef __PASO_SPARSEMATRIX_H__
30#define __PASO_SPARSEMATRIX_H__
36template <
typename T>
struct SparseMatrix;
44struct SparseMatrix : boost::enable_shared_from_this<SparseMatrix<T> >
47 dim_t rowBlockSize, dim_t colBlockSize,
48 bool patternIsUnrolled);
63 const double* b)
const;
65 void invMain(
double* inv_diag, index_t* pivot)
const;
70 const index_t* row_list,
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);
201 const double beta, T* out);
206 const double beta,
double* out);
224struct Preconditioner_LocalSmoother;
245 dim_t rowBlockSize, dim_t colBlockSize,
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);
324 switch (solver_package) {
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;
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;
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);
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);
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);
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);
#define PASO_SMOOTHER
Definition Options.h:75
#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
PasoException exception class.
Definition PasoException.h:34
Definition BiCGStab.cpp:25
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
boost::shared_ptr< SparseMatrix< T > > SparseMatrix_ptr
Definition SparseMatrix.h:37
boost::shared_ptr< Pattern > Pattern_ptr
Definition Pattern.h:40
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
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 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 UMFPACK_free(SparseMatrix< double > *A)
frees any UMFPACK related data from the matrix
Definition UMFPACK.cpp:35
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
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
boost::shared_ptr< const SparseMatrix< T > > const_SparseMatrix_ptr
Definition SparseMatrix.h:38
void Preconditioner_LocalSmoother_free(Preconditioner_LocalSmoother *in)
Definition Smoother.cpp:42
void MUMPS_free(SparseMatrix< T > *A)
frees any MUMPS related data from the matrix
Definition MUMPS.h:118
void MKL_free(SparseMatrix< double > *A)
Definition MKL.cpp:35
SparseMatrix_ptr< double > SparseMatrix_MatrixMatrix(const_SparseMatrix_ptr< double > A, const_SparseMatrix_ptr< double > B)
Definition SparseMatrix_MatrixMatrix.cpp:44
int SparseMatrixType
Definition SparseMatrix.h:40
#define PASO_DLL_API
Definition paso/src/system_dep.h:29
Definition Preconditioner.h:57
Definition SparseMatrix.h:45
index_t solver_package
package controlling the solver pointer
Definition SparseMatrix.h:174
void invMain(double *inv_diag, index_t *pivot) const
SparseMatrixType type
Definition SparseMatrix.h:161
void copyBlockToMainDiagonal(const double *in)
void copyFromMainDiagonal(double *out) const
dim_t len
Definition SparseMatrix.h:168
double getSparsity() const
Definition SparseMatrix.h:126
dim_t block_size
Definition SparseMatrix.h:164
dim_t getTotalNumCols() const
Definition SparseMatrix.h:106
SparseMatrix_ptr< double > unroll(SparseMatrixType type) const
void nullifyRowsAndCols_CSC_BLK1(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition SparseMatrix.h:450
Pattern_ptr pattern
Definition SparseMatrix.h:167
dim_t numRows
Definition SparseMatrix.h:165
T * val
this is used for classical CSR or CSC
Definition SparseMatrix.h:171
dim_t getTotalNumRows() const
Definition SparseMatrix.h:101
void applyBlockMatrix(double *block_diag, index_t *pivot, double *x, const double *b) const
void saveMM(const char *filename) const
void applyDiagonal_CSR_OFFSET0(const double *left, const double *right)
dim_t numCols
Definition SparseMatrix.h:166
void addRow_CSR_OFFSET0(double *array) const
void saveHB_CSC(const char *filename) const
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
static SparseMatrix_ptr< double > loadMM_toCSR(const char *filename)
void copyBlockFromMainDiagonal(double *out) const
dim_t maxDeg() const
Definition SparseMatrix.h:96
void nullifyRows_CSR(const double *mask_row, double main_diagonal_value)
void * solver_p
pointer to data needed by a solver
Definition SparseMatrix.h:177
SparseMatrix_ptr< double > getTranspose() const
void nullifyRowsAndCols_CSR_BLK1(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition SparseMatrix.h:494
SparseMatrix(SparseMatrixType type, Pattern_ptr pattern, dim_t rowBlockSize, dim_t colBlockSize, bool patternIsUnrolled)
Definition SparseMatrix.h:244
dim_t col_block_size
Definition SparseMatrix.h:163
index_t * borrowMainDiagonalPointer() const
Definition SparseMatrix.h:81
void maxAbsRow_CSR_OFFSET0(double *array) const
void nullifyRows_CSR_BLK1(const double *mask_row, double main_diagonal_value)
~SparseMatrix()
Definition SparseMatrix.h:322
dim_t getNumRows() const
Definition SparseMatrix.h:111
dim_t row_block_size
Definition SparseMatrix.h:162
index_t * borrowColoringPointer() const
Definition SparseMatrix.h:86
void nullifyRowsAndCols_CSC(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition SparseMatrix.h:469
void nullifyRowsAndCols_CSR(const double *mask_row, const double *mask_col, double main_diagonal_value)
Definition SparseMatrix.h:513
void setValues(T value)
Definition SparseMatrix.h:345
double getSize() const
Definition SparseMatrix.h:121
void copyToMainDiagonal(const double *in)
void addAbsRow_CSR_OFFSET0(double *array) const
dim_t getNumColors() const
Definition SparseMatrix.h:91
dim_t getNumCols() const
Definition SparseMatrix.h:116
SparseMatrix_ptr< double > getBlock(int blockid) const