escript
Revision_
|
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