escript  Revision_
SystemMatrix.h
Go to the documentation of this file.
00001 
00002 /*****************************************************************************
00003 *
00004 * Copyright (c) 2003-2014 by University of Queensland
00005 * http://www.uq.edu.au
00006 *
00007 * Primary Business: Queensland, Australia
00008 * Licensed under the Open Software License version 3.0
00009 * http://www.opensource.org/licenses/osl-3.0.php
00010 *
00011 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
00012 * Development 2012-2013 by School of Earth Sciences
00013 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
00014 *
00015 *****************************************************************************/
00016 
00017 
00018 /****************************************************************************/
00019 
00020 /*   Paso: SystemMatrix */
00021 
00022 /****************************************************************************/
00023 
00024 /*   Copyrights by ACcESS Australia 2003,2004,2005,2006 */
00025 /*   Author: Lutz Gross, l.gross@uq.edu.au */
00026 
00027 /****************************************************************************/
00028 
00029 #ifndef __PASO_SYSTEMMATRIX_H__
00030 #define __PASO_SYSTEMMATRIX_H__
00031 
00032 #include "SparseMatrix.h"
00033 #include "SystemMatrixPattern.h"
00034 #include "Options.h"
00035 
00036 namespace paso {
00037 
00038 struct SystemMatrix;
00039 typedef boost::shared_ptr<SystemMatrix> SystemMatrix_ptr;
00040 typedef boost::shared_ptr<const SystemMatrix> const_SystemMatrix_ptr;
00041 
00042 typedef int SystemMatrixType;
00043 
00044 //  this struct holds a (distributed) stiffness matrix
00045 PASO_DLL_API
00046 struct SystemMatrix : boost::enable_shared_from_this<SystemMatrix>
00047 {
00048     SystemMatrix(SystemMatrixType, SystemMatrixPattern_ptr, dim_t, dim_t,
00049                  bool patternIsUnrolled);
00050 
00051     ~SystemMatrix();
00052 
00057     void nullifyRowsAndCols(double* mask_row, double* mask_col,
00058                             double main_diagonal_value);
00059 
00064     void nullifyRows(double* mask_row, double main_diagonal_value);
00065 
00066     void add(dim_t, index_t*, dim_t, dim_t, index_t*, dim_t, double*);
00067 
00068     void makeZeroRowSums(double* left_over);
00069 
00078     void copyColCoupleBlock();
00079 
00080     void copyRemoteCoupleBlock(bool recreatePattern);
00081 
00082     void fillWithGlobalCoordinates(double f1);
00083 
00084     void print() const;
00085 
00088     SparseMatrix_ptr mergeSystemMatrix() const;
00089 
00090     void mergeMainAndCouple(index_t** p_ptr, index_t** p_idx, double** p_val) const;
00091 
00092     void mergeMainAndCouple_CSR_OFFSET0(index_t** p_ptr, index_t** p_idx, double** p_val) const;
00093     void mergeMainAndCouple_CSR_OFFSET0_Block(index_t** p_ptr, index_t** p_idx, double** p_val) const;
00094 
00095     void mergeMainAndCouple_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val) const;
00096 
00097     void copyMain_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val);
00098 
00099     void extendedRowsForST(dim_t* degree_ST, index_t* offset_ST, index_t* ST);
00100 
00101     void applyBalanceInPlace(double* x, bool RHS) const;
00102 
00103     void applyBalance(double* x_out, const double* x, bool RHS) const;
00104 
00105     void balance();
00106 
00107     double getGlobalSize() const;
00108 
00109     void setPreconditioner(Options* options);
00110 
00115     void solvePreconditioner(double* x, double* b);
00116 
00117     void freePreconditioner();
00118 
00119     index_t* borrowMainDiagonalPointer() const;
00120 
00121     inline void startCollect(const double* in)
00122     {
00123         startColCollect(in);
00124     }
00125 
00126     inline double* finishCollect()
00127     {
00128         return finishColCollect();
00129     }
00130 
00131     inline void startColCollect(const double* in)
00132     {
00133         col_coupler->startCollect(in);
00134     }
00135 
00136     inline double* finishColCollect()
00137     {
00138         return col_coupler->finishCollect();
00139     }
00140 
00141     inline void startRowCollect(const double* in)
00142     {
00143         row_coupler->startCollect(in);
00144     }
00145 
00146     inline double* finishRowCollect()
00147     {
00148         return row_coupler->finishCollect();
00149     }
00150 
00151     inline dim_t getNumRows() const
00152     {
00153         return mainBlock->numRows;
00154     }
00155 
00156     inline dim_t getNumCols() const
00157     {
00158         return mainBlock->numCols;
00159     }
00160 
00161     inline dim_t getTotalNumRows() const
00162     {
00163         return getNumRows() * row_block_size;
00164     }
00165 
00166     inline dim_t getTotalNumCols() const
00167     {
00168         return getNumCols() * col_block_size;
00169     }
00170 
00171     inline dim_t getRowOverlap()  const
00172     {
00173         return row_coupler->getNumOverlapComponents();
00174     }
00175 
00176     inline dim_t getColOverlap() const
00177     {
00178         return col_coupler->getNumOverlapComponents();
00179     }
00180 
00181     inline dim_t getGlobalNumRows() const
00182     {
00183         if (type & MATRIX_FORMAT_CSC) {
00184             return pattern->input_distribution->getGlobalNumComponents();
00185         }
00186         return pattern->output_distribution->getGlobalNumComponents();
00187     }
00188 
00189     inline dim_t getGlobalNumCols() const
00190     {
00191         if (type & MATRIX_FORMAT_CSC) {
00192             return pattern->output_distribution->getGlobalNumComponents();
00193         }
00194         return pattern->input_distribution->getGlobalNumComponents();
00195     }
00196 
00197     inline dim_t getGlobalTotalNumRows() const
00198     {
00199         return getGlobalNumRows() * row_block_size;
00200     }
00201 
00202     inline dim_t getGlobalTotalNumCols() const
00203     {
00204         return getGlobalNumCols() * col_block_size;
00205     }
00206 
00207     inline double getSparsity() const
00208     {
00209         return getGlobalSize() /
00210                  ((double)getGlobalTotalNumRows()*getGlobalTotalNumCols());
00211     }
00212 
00213     inline dim_t getNumOutput() const
00214     {
00215        return pattern->getNumOutput();
00216     }
00217 
00218     inline void copyBlockFromMainDiagonal(double* out) const
00219     {
00220         mainBlock->copyBlockFromMainDiagonal(out);
00221     }
00222 
00223     inline void copyBlockToMainDiagonal(const double* in)
00224     {
00225         mainBlock->copyBlockToMainDiagonal(in);
00226     }
00227 
00228     inline void copyFromMainDiagonal(double* out) const
00229     {
00230         mainBlock->copyFromMainDiagonal(out);
00231     }
00232 
00233     inline void copyToMainDiagonal(const double* in)
00234     {
00235         mainBlock->copyToMainDiagonal(in);
00236     }
00237 
00238     inline void setValues(double value)
00239     {
00240         mainBlock->setValues(value);
00241         col_coupleBlock->setValues(value);
00242         row_coupleBlock->setValues(value);
00243         is_balanced = false;
00244     }
00245 
00246     inline void saveMM(const char* filename) const
00247     {
00248         if (mpi_info->size > 1) {
00249             Esys_setError(IO_ERROR, "SystemMatrix::saveMM: Only single rank supported.");
00250         } else {
00251             mainBlock->saveMM(filename);
00252         }
00253     }
00254 
00255     inline void saveHB(const char *filename) const
00256     {
00257         if (mpi_info->size > 1) {
00258             Esys_setError(TYPE_ERROR, "SystemMatrix::saveHB: Only single rank supported.");
00259         } else if (!(type & MATRIX_FORMAT_CSC)) {
00260             Esys_setError(TYPE_ERROR, "SystemMatrix::saveHB: Only CSC format supported.");
00261         } else {
00262             mainBlock->saveHB_CSC(filename);
00263         }
00264     }
00265 
00266     inline void rowSum(double* row_sum) const
00267     {
00268         if ((type & MATRIX_FORMAT_CSC) || (type & MATRIX_FORMAT_OFFSET1)) {
00269             Esys_setError(TYPE_ERROR, "SystemMatrix::rowSum: No normalization "
00270                   "available for compressed sparse column or index offset 1.");
00271         } else {
00272             const dim_t nrow = mainBlock->numRows*row_block_size;
00273 #pragma omp parallel for
00274             for (index_t irow=0; irow<nrow; ++irow) {
00275                 row_sum[irow]=0.;
00276             }
00277             mainBlock->addRow_CSR_OFFSET0(row_sum);
00278             col_coupleBlock->addRow_CSR_OFFSET0(row_sum);
00279         }
00280     }
00281 
00282     static SystemMatrix_ptr loadMM_toCSR(const char* filename);
00283 
00284     static SystemMatrix_ptr loadMM_toCSC(const char* filename);
00285 
00286     static index_t getSystemMatrixTypeId(index_t solver,
00287                                          index_t preconditioner,
00288                                          index_t package, bool symmetry,
00289                                          Esys_MPIInfo* mpi_info);
00290 
00291     SystemMatrixType type;
00292     SystemMatrixPattern_ptr pattern;
00293 
00294     dim_t logical_row_block_size;
00295     dim_t logical_col_block_size;
00296 
00297     dim_t row_block_size;
00298     dim_t col_block_size;
00299     dim_t block_size;
00300 
00301     Distribution_ptr row_distribution;
00302     Distribution_ptr col_distribution;
00303     Esys_MPIInfo *mpi_info;
00304 
00305     Coupler_ptr col_coupler;
00306     Coupler_ptr row_coupler;
00307 
00309     SparseMatrix_ptr mainBlock;
00311     SparseMatrix_ptr col_coupleBlock;
00313     SparseMatrix_ptr row_coupleBlock;
00315     SparseMatrix_ptr remote_coupleBlock;
00316 
00317     bool is_balanced;
00318 
00324     double* balance_vector;
00325 
00327     mutable index_t* global_id;
00328 
00330     index_t solver_package;
00331 
00333     void* solver_p;
00334 
00336     void* trilinos_data;
00337 };
00338 
00339 
00340 void SystemMatrix_MatrixVector(double alpha, SystemMatrix_ptr A, const double* in, double beta, double* out);
00341 
00342 void SystemMatrix_MatrixVector_CSR_OFFSET0(double alpha, SystemMatrix_ptr A, const double* in, double beta, double* out);
00343 
00344 void RHS_loadMM_toCSR(const char* filename, double* b, dim_t size);
00345 
00346 
00347 } // namespace paso
00348 
00349 #endif // __PASO_SYSTEMMATRIX_H__
00350