SystemMatrix.h
Go to the documentation of this file.
/*****************************************************************************
*
* Copyright (c) 2003-2017 by The University of Queensland
* http://www.uq.edu.au
*
* Primary Business: Queensland, Australia
* Licensed under the Apache License, version 2.0
* http://www.apache.org/licenses/LICENSE-2.0
*
* Development until 2012 by Earth Systems Science Computational Center (ESSCC)
* Development 2012-2013 by School of Earth Sciences
* Development from 2014 by Centre for Geoscience Computing (GeoComp)
*
*****************************************************************************/
/****************************************************************************/
/* Paso: SystemMatrix */
/****************************************************************************/
/* Copyrights by ACcESS Australia 2003,2004,2005,2006 */
/* Author: Lutz Gross, l.gross@uq.edu.au */
/****************************************************************************/
#ifndef __PASO_SYSTEMMATRIX_H__
#define __PASO_SYSTEMMATRIX_H__
#include "SparseMatrix.h"
#include "SystemMatrixPattern.h"
#include <escript/AbstractSystemMatrix.h>
namespace paso {
class Options;
class SystemMatrix;
typedef boost::shared_ptr<SystemMatrix> SystemMatrix_ptr;
typedef boost::shared_ptr<const SystemMatrix> const_SystemMatrix_ptr;
typedef int SystemMatrixType;
class SystemMatrix : public escript::AbstractSystemMatrix
{
public:
SystemMatrix();
SystemMatrix(SystemMatrixType type, SystemMatrixPattern_ptr pattern,
dim_t rowBlockSize, dim_t columnBlockSize,
bool patternIsUnrolled, const escript::FunctionSpace& rowFS,
const escript::FunctionSpace& colFS);
~SystemMatrix();
virtual void nullifyRowsAndCols(escript::Data& mask_row,
escript::Data& mask_col,
double main_diagonal_value);
virtual inline void saveMM(const std::string& filename) const
{
if (mpi_info->size > 1) {
//throw PasoException("SystemMatrix::saveMM: Only single rank supported.");
SparseMatrix_ptr merged(mergeSystemMatrix());
if (mpi_info->rank == 0)
merged->saveMM(filename.c_str());
} else {
mainBlock->saveMM(filename.c_str());
}
}
virtual inline void saveHB(const std::string& filename) const
{
if (mpi_info->size > 1) {
throw PasoException("SystemMatrix::saveHB: Only single rank supported.");
} else if (!(type & MATRIX_FORMAT_CSC)) {
throw PasoException("SystemMatrix::saveHB: Only CSC format supported.");
} else {
mainBlock->saveHB_CSC(filename.c_str());
}
}
virtual void resetValues(bool preserveSolverData = false);
void nullifyRows(double* mask_row, double main_diagonal_value);
void add(dim_t, index_t*, dim_t, dim_t, index_t*, dim_t, double*);
void makeZeroRowSums(double* left_over);
void copyColCoupleBlock();
void copyRemoteCoupleBlock(bool recreatePattern);
void fillWithGlobalCoordinates(double f1);
void print() const;
SparseMatrix_ptr mergeSystemMatrix() const;
void mergeMainAndCouple(index_t** p_ptr, index_t** p_idx, double** p_val) const;
void mergeMainAndCouple_CSR_OFFSET0(index_t** p_ptr, index_t** p_idx, double** p_val) const;
void mergeMainAndCouple_CSR_OFFSET0_Block(index_t** p_ptr, index_t** p_idx, double** p_val) const;
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);
void extendedRowsForST(dim_t* degree_ST, index_t* offset_ST, index_t* ST);
void applyBalanceInPlace(double* x, bool RHS) const;
void applyBalance(double* x_out, const double* x, bool RHS) const;
void balance();
double getGlobalSize() const;
void setPreconditioner(Options* options);
void solvePreconditioner(double* x, double* b);
void freePreconditioner();
index_t* borrowMainDiagonalPointer() const;
inline void startCollect(const double* in) const
{
startColCollect(in);
}
inline double* finishCollect() const
{
return finishColCollect();
}
inline void startColCollect(const double* in) const
{
col_coupler->startCollect(in);
}
inline double* finishColCollect() const
{
return col_coupler->finishCollect();
}
inline void startRowCollect(const double* in)
{
row_coupler->startCollect(in);
}
inline double* finishRowCollect()
{
return row_coupler->finishCollect();
}
inline dim_t getNumRows() const
{
return mainBlock->numRows;
}
inline dim_t getNumCols() const
{
return mainBlock->numCols;
}
inline dim_t getTotalNumRows() const
{
return getNumRows() * row_block_size;
}
inline dim_t getTotalNumCols() const
{
return getNumCols() * col_block_size;
}
inline dim_t getRowOverlap() const
{
return row_coupler->getNumOverlapComponents();
}
inline dim_t getColOverlap() const
{
return col_coupler->getNumOverlapComponents();
}
inline dim_t getGlobalNumRows() const
{
if (type & MATRIX_FORMAT_CSC) {
return pattern->input_distribution->getGlobalNumComponents();
}
return pattern->output_distribution->getGlobalNumComponents();
}
inline dim_t getGlobalNumCols() const
{
if (type & MATRIX_FORMAT_CSC) {
return pattern->output_distribution->getGlobalNumComponents();
}
return pattern->input_distribution->getGlobalNumComponents();
}
inline dim_t getGlobalTotalNumRows() const
{
return getGlobalNumRows() * row_block_size;
}
inline dim_t getGlobalTotalNumCols() const
{
return getGlobalNumCols() * col_block_size;
}
inline double getSparsity() const
{
return getGlobalSize() /
((double)getGlobalTotalNumRows()*getGlobalTotalNumCols());
}
inline dim_t getNumOutput() const
{
return pattern->getNumOutput();
}
inline void copyBlockFromMainDiagonal(double* out) const
{
mainBlock->copyBlockFromMainDiagonal(out);
}
inline void copyBlockToMainDiagonal(const double* in)
{
mainBlock->copyBlockToMainDiagonal(in);
}
inline void copyFromMainDiagonal(double* out) const
{
mainBlock->copyFromMainDiagonal(out);
}
inline void copyToMainDiagonal(const double* in)
{
mainBlock->copyToMainDiagonal(in);
}
inline void setValues(double value)
{
mainBlock->setValues(value);
col_coupleBlock->setValues(value);
row_coupleBlock->setValues(value);
is_balanced = false;
}
inline void rowSum(double* row_sum) const
{
if ((type & MATRIX_FORMAT_CSC) || (type & MATRIX_FORMAT_OFFSET1)) {
throw PasoException("SystemMatrix::rowSum: No normalization "
"available for compressed sparse column or index offset 1.");
} else {
const dim_t nrow = mainBlock->numRows*row_block_size;
#pragma omp parallel for
for (index_t irow=0; irow<nrow; ++irow) {
row_sum[irow]=0.;
}
mainBlock->addRow_CSR_OFFSET0(row_sum);
col_coupleBlock->addRow_CSR_OFFSET0(row_sum);
}
}
void MatrixVector(double alpha, const double* in, double beta,
double* out) const;
void MatrixVector_CSR_OFFSET0(double alpha, const double* in, double beta,
double* out) const;
static SystemMatrix_ptr loadMM_toCSR(const char* filename);
static SystemMatrix_ptr loadMM_toCSC(const char* filename);
static int getSystemMatrixTypeId(int solver, int preconditioner,
int package, bool symmetry,
const escript::JMPI& mpi_info);
SystemMatrixType type;
SystemMatrixPattern_ptr pattern;
dim_t logical_row_block_size;
dim_t logical_col_block_size;
dim_t row_block_size;
dim_t col_block_size;
dim_t block_size;
escript::Distribution_ptr row_distribution;
escript::Distribution_ptr col_distribution;
escript::JMPI mpi_info;
Coupler_ptr<real_t> col_coupler;
Coupler_ptr<real_t> row_coupler;
SparseMatrix_ptr mainBlock;
SparseMatrix_ptr col_coupleBlock;
SparseMatrix_ptr row_coupleBlock;
SparseMatrix_ptr remote_coupleBlock;
bool is_balanced;
double* balance_vector;
mutable index_t* global_id;
mutable index_t solver_package;
void* solver_p;
private:
virtual void setToSolution(escript::Data& out, escript::Data& in,
boost::python::object& options) const;
virtual void ypAx(escript::Data& y, escript::Data& x) const;
void solve(double* out, double* in, Options* options) const;
};
void RHS_loadMM_toCSR(const char* filename, double* b, dim_t size);
} // namespace paso
#endif // __PASO_SYSTEMMATRIX_H__
paso::SystemMatrix::startColCollect
void startColCollect(const double *in) const
Definition: SystemMatrix.h:164
paso::SystemMatrix::row_coupler
Coupler_ptr< real_t > row_coupler
Definition: SystemMatrix.h:324
paso::SystemMatrix::col_coupleBlock
SparseMatrix_ptr col_coupleBlock
coupling to neighbouring processors (row - col)
Definition: SystemMatrix.h:329
escript::FunctionSpace
Definition: FunctionSpace.h:34
paso::SystemMatrix::getSystemMatrixTypeId
static int getSystemMatrixTypeId(int solver, int preconditioner, int package, bool symmetry, const escript::JMPI &mpi_info)
Definition: SystemMatrix.cpp:532
paso::SystemMatrix::getNumRows
dim_t getNumRows() const
Definition: SystemMatrix.h:184
paso::SystemMatrix::extendedRowsForST
void extendedRowsForST(dim_t *degree_ST, index_t *offset_ST, index_t *ST)
Definition: SystemMatrix_extendedRows.cpp:41
paso::SystemMatrix::getTotalNumCols
dim_t getTotalNumCols() const
Definition: SystemMatrix.h:199
paso::SystemMatrix::balance_vector
double * balance_vector
Definition: SystemMatrix.h:342
paso::SystemMatrix::copyBlockToMainDiagonal
void copyBlockToMainDiagonal(const double *in)
Definition: SystemMatrix.h:256
paso::SystemMatrix::row_coupleBlock
SparseMatrix_ptr row_coupleBlock
coupling to neighbouring processors (col - row)
Definition: SystemMatrix.h:331
paso::SystemMatrix::getColOverlap
dim_t getColOverlap() const
Definition: SystemMatrix.h:209
paso::SystemMatrix::~SystemMatrix
~SystemMatrix()
Definition: SystemMatrix.cpp:142
paso::SystemMatrix::nullifyRows
void nullifyRows(double *mask_row, double main_diagonal_value)
Definition: SystemMatrix.cpp:223
paso::SystemMatrix::print
void print() const
Definition: SystemMatrix_debug.cpp:100
paso::SystemMatrix::logical_col_block_size
dim_t logical_col_block_size
Definition: SystemMatrix.h:313
paso::Options
Definition: Options.h:89
paso::SystemMatrix::solve
void solve(double *out, double *in, Options *options) const
Definition: solve.cpp:39
MATRIX_FORMAT_OFFSET1
#define MATRIX_FORMAT_OFFSET1
Definition: Paso.h:64
paso::SystemMatrix::mergeMainAndCouple
void mergeMainAndCouple(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:40
escript::Distribution_ptr
boost::shared_ptr< Distribution > Distribution_ptr
Definition: Distribution.h:24
paso::SystemMatrix::pattern
SystemMatrixPattern_ptr pattern
Definition: SystemMatrix.h:310
paso::SystemMatrix::copyMain_CSC_OFFSET1
void copyMain_CSC_OFFSET1(index_t **p_ptr, index_t **p_idx, double **p_val)
paso::SystemMatrix::getTotalNumRows
dim_t getTotalNumRows() const
Definition: SystemMatrix.h:194
paso::SystemMatrix::mergeMainAndCouple_CSC_OFFSET1
void mergeMainAndCouple_CSC_OFFSET1(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:276
paso::SystemMatrix::copyBlockFromMainDiagonal
void copyBlockFromMainDiagonal(double *out) const
Definition: SystemMatrix.h:251
paso::SystemMatrix::mergeMainAndCouple_CSR_OFFSET0_Block
void mergeMainAndCouple_CSR_OFFSET0_Block(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:168
paso::SystemMatrix::block_size
dim_t block_size
Definition: SystemMatrix.h:317
escript::JMPI
boost::shared_ptr< JMPI_ > JMPI
Definition: EsysMPI.h:71
paso::SystemMatrix::col_distribution
escript::Distribution_ptr col_distribution
Definition: SystemMatrix.h:320
paso::SystemMatrix::getGlobalNumRows
dim_t getGlobalNumRows() const
Definition: SystemMatrix.h:214
paso::SystemMatrix::startCollect
void startCollect(const double *in) const
Definition: SystemMatrix.h:154
paso::SystemMatrix::borrowMainDiagonalPointer
index_t * borrowMainDiagonalPointer() const
Definition: SystemMatrix.cpp:187
paso::SystemMatrix::applyBalanceInPlace
void applyBalanceInPlace(double *x, bool RHS) const
Definition: SystemMatrix.cpp:443
paso::SparseMatrix_ptr
boost::shared_ptr< SparseMatrix > SparseMatrix_ptr
Definition: SparseMatrix.h:35
paso::SystemMatrix::finishCollect
double * finishCollect() const
Definition: SystemMatrix.h:159
paso::SystemMatrix::add
void add(dim_t, index_t *, dim_t, dim_t, index_t *, dim_t, double *)
paso::SystemMatrix::copyToMainDiagonal
void copyToMainDiagonal(const double *in)
Definition: SystemMatrix.h:266
paso::SystemMatrix::getNumOutput
dim_t getNumOutput() const
Definition: SystemMatrix.h:246
paso::SystemMatrix::saveMM
virtual void saveMM(const std::string &filename) const
writes the matrix to a file using the Matrix Market file format
Definition: SystemMatrix.h:68
paso::SystemMatrix::setValues
void setValues(double value)
Definition: SystemMatrix.h:271
paso::SystemMatrix::getSparsity
double getSparsity() const
Definition: SystemMatrix.h:240
paso::SystemMatrix::saveHB
virtual void saveHB(const std::string &filename) const
writes the matrix to a file using the Harwell-Boeing file format
Definition: SystemMatrix.h:80
SparseMatrix.h
paso::SystemMatrix::row_distribution
escript::Distribution_ptr row_distribution
Definition: SystemMatrix.h:319
paso::SystemMatrixPattern_ptr
boost::shared_ptr< SystemMatrixPattern > SystemMatrixPattern_ptr
Definition: SystemMatrixPattern.h:39
paso::RHS_loadMM_toCSR
void RHS_loadMM_toCSR(const char *filename, double *b, dim_t size)
Definition: SystemMatrix_loadMM.cpp:299
paso::SystemMatrix::is_balanced
bool is_balanced
Definition: SystemMatrix.h:335
paso::SystemMatrix_ptr
boost::shared_ptr< SystemMatrix > SystemMatrix_ptr
Definition: SystemMatrix.h:40
paso::SystemMatrix::copyFromMainDiagonal
void copyFromMainDiagonal(double *out) const
Definition: SystemMatrix.h:261
paso::SystemMatrix::mainBlock
SparseMatrix_ptr mainBlock
main block
Definition: SystemMatrix.h:327
paso::SystemMatrix::row_block_size
dim_t row_block_size
Definition: SystemMatrix.h:315
paso::SystemMatrix::solver_p
void * solver_p
pointer to data needed by a solver
Definition: SystemMatrix.h:351
paso::SystemMatrix::freePreconditioner
void freePreconditioner()
Definition: SystemMatrix.cpp:164
paso
Definition: AMG.cpp:45
paso::SystemMatrix::loadMM_toCSC
static SystemMatrix_ptr loadMM_toCSC(const char *filename)
Definition: SystemMatrix_loadMM.cpp:200
escript::DataTypes::index_t
int index_t
type for array/matrix indices used both globally and on each rank
Definition: DataTypes.h:59
paso::SystemMatrix::copyRemoteCoupleBlock
void copyRemoteCoupleBlock(bool recreatePattern)
Definition: SystemMatrix_copyRemoteCoupleBlock.cpp:39
paso::SystemMatrix::setToSolution
virtual void setToSolution(escript::Data &out, escript::Data &in, boost::python::object &options) const
solves the linear system this*out=in
Definition: SystemMatrix.cpp:313
paso::SystemMatrix::col_coupler
Coupler_ptr< real_t > col_coupler
Definition: SystemMatrix.h:323
paso::SystemMatrix::mergeMainAndCouple_CSR_OFFSET0
void mergeMainAndCouple_CSR_OFFSET0(index_t **p_ptr, index_t **p_idx, double **p_val) const
Definition: SystemMatrix_mergeMainAndCouple.cpp:56
paso::SystemMatrix::startRowCollect
void startRowCollect(const double *in)
Definition: SystemMatrix.h:174
paso::SystemMatrix::setPreconditioner
void setPreconditioner(Options *options)
Definition: SystemMatrix.cpp:149
paso::SystemMatrix
this class holds a (distributed) stiffness matrix
Definition: SystemMatrix.h:47
escript::Data
Data represents a collection of datapoints.
Definition: Data.h:63
paso::SystemMatrix::solvePreconditioner
void solvePreconditioner(double *x, double *b)
Definition: SystemMatrix.cpp:157
paso::SystemMatrix::ypAx
virtual void ypAx(escript::Data &y, escript::Data &x) const
performs y+=this*x
Definition: SystemMatrix.cpp:341
paso::SystemMatrix::balance
void balance()
Definition: SystemMatrix.cpp:481
paso::SystemMatrix::getGlobalSize
double getGlobalSize() const
Definition: SystemMatrix.cpp:171
paso::const_SystemMatrix_ptr
boost::shared_ptr< const SystemMatrix > const_SystemMatrix_ptr
Definition: SystemMatrix.h:42
paso::SystemMatrix::global_id
index_t * global_id
stores the global ids for all cols in col_coupleBlock
Definition: SystemMatrix.h:345
paso::SystemMatrixType
int SystemMatrixType
Definition: SystemMatrix.h:44
paso::SystemMatrix::remote_coupleBlock
SparseMatrix_ptr remote_coupleBlock
coupling of rows-cols on neighbouring processors (may not be valid)
Definition: SystemMatrix.h:333
paso::SystemMatrix::logical_row_block_size
dim_t logical_row_block_size
Definition: SystemMatrix.h:312
SystemMatrixPattern.h
paso::SystemMatrix::fillWithGlobalCoordinates
void fillWithGlobalCoordinates(double f1)
Definition: SystemMatrix_debug.cpp:37
paso::SystemMatrix::nullifyRowsAndCols
virtual void nullifyRowsAndCols(escript::Data &mask_row, escript::Data &mask_col, double main_diagonal_value)
Definition: SystemMatrix.cpp:244
paso::SystemMatrix::applyBalance
void applyBalance(double *x_out, const double *x, bool RHS) const
Definition: SystemMatrix.cpp:462
paso::SystemMatrix::getNumCols
dim_t getNumCols() const
Definition: SystemMatrix.h:189
paso::SystemMatrix::type
SystemMatrixType type
Definition: SystemMatrix.h:309
paso::SystemMatrix::getGlobalNumCols
dim_t getGlobalNumCols() const
Definition: SystemMatrix.h:222
paso::SystemMatrix::getGlobalTotalNumRows
dim_t getGlobalTotalNumRows() const
Definition: SystemMatrix.h:230
escript::AbstractSystemMatrix
Base class for escript system matrices.
Definition: AbstractSystemMatrix.h:42
paso::SystemMatrix::solver_package
index_t solver_package
package code controlling the solver pointer
Definition: SystemMatrix.h:348
MATRIX_FORMAT_CSC
#define MATRIX_FORMAT_CSC
Definition: Paso.h:62
paso::SystemMatrix::mpi_info
escript::JMPI mpi_info
Definition: SystemMatrix.h:321
paso::PasoException
PasoException exception class.
Definition: PasoException.h:32
paso::SystemMatrix::getRowOverlap
dim_t getRowOverlap() const
Definition: SystemMatrix.h:204
paso::SystemMatrix::MatrixVector
void MatrixVector(double alpha, const double *in, double beta, double *out) const
Definition: SystemMatrix_MatrixVector.cpp:34
paso::SystemMatrix::makeZeroRowSums
void makeZeroRowSums(double *left_over)
Definition: SystemMatrix.cpp:201
paso::SystemMatrix::finishColCollect
double * finishColCollect() const
Definition: SystemMatrix.h:169
paso::SystemMatrix::rowSum
void rowSum(double *row_sum) const
Definition: SystemMatrix.h:279
paso::Coupler_ptr
boost::shared_ptr< Coupler< T > > Coupler_ptr
Definition: Coupler.h:42
paso::SystemMatrix::col_block_size
dim_t col_block_size
Definition: SystemMatrix.h:316
paso::SystemMatrix::finishRowCollect
double * finishRowCollect()
Definition: SystemMatrix.h:179
paso::SystemMatrix::SystemMatrix
SystemMatrix()
default constructor - throws exception.
Definition: SystemMatrix.cpp:41
paso::SystemMatrix::MatrixVector_CSR_OFFSET0
void MatrixVector_CSR_OFFSET0(double alpha, const double *in, double beta, double *out) const
Definition: SystemMatrix_MatrixVector.cpp:63
paso::SystemMatrix::resetValues
virtual void resetValues(bool preserveSolverData=false)
resets the matrix entries
Definition: SystemMatrix.cpp:306
escript::DataTypes::dim_t
index_t dim_t
Definition: DataTypes.h:64
paso::SystemMatrix::copyColCoupleBlock
void copyColCoupleBlock()
Definition: SystemMatrix.cpp:365
paso::SystemMatrix::loadMM_toCSR
static SystemMatrix_ptr loadMM_toCSR(const char *filename)
Definition: SystemMatrix_loadMM.cpp:100
paso::SystemMatrix::mergeSystemMatrix
SparseMatrix_ptr mergeSystemMatrix() const
Definition: SystemMatrix.cpp:566
paso::SystemMatrix::getGlobalTotalNumCols
dim_t getGlobalTotalNumCols() const
Definition: SystemMatrix.h:235