/[escript]/trunk/paso/src/SystemMatrix.h
ViewVC logotype

Diff of /trunk/paso/src/SystemMatrix.h

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4835 by caltinay, Thu Apr 3 04:02:53 2014 UTC revision 4836 by caltinay, Mon Apr 7 05:51:55 2014 UTC
# Line 17  Line 17 
17    
18  /****************************************************************************/  /****************************************************************************/
19    
20  /*   Paso: SystemMatrix and SystemVector */  /*   Paso: SystemMatrix */
21    
22  /****************************************************************************/  /****************************************************************************/
23    
# Line 26  Line 26 
26    
27  /****************************************************************************/  /****************************************************************************/
28    
29  #ifndef INC_PASO_SYSTEMMATRIX  #ifndef __PASO_SYSTEMMATRIX_H__
30  #define INC_PASO_SYSTEMMATRIX  #define __PASO_SYSTEMMATRIX_H__
31    
32  #include "Common.h"  #include "Common.h"
33    #include "Coupler.h"
34  #include "SparseMatrix.h"  #include "SparseMatrix.h"
35  #include "SystemMatrixPattern.h"  #include "SystemMatrixPattern.h"
36  #include "Options.h"  #include "Options.h"
37  #include "esysUtils/Esys_MPI.h"  #include "esysUtils/Esys_MPI.h"
 #include "Paso.h"  
 #include "Coupler.h"  
   
   
 typedef int Paso_SystemMatrixType;  
   
 //  this struct holds a stiffness matrix  
 struct Paso_SystemMatrix  
 {  
   Paso_SystemMatrixType type;  
   paso::SystemMatrixPattern_ptr pattern;  
   
   dim_t reference_counter;  
   
   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;  
   
   paso::Distribution_ptr row_distribution;  
   paso::Distribution_ptr col_distribution;  
   Esys_MPIInfo *mpi_info;  
   
   paso::Coupler_ptr col_coupler;  
   paso::Coupler_ptr row_coupler;  
   
   /* this comes into play when PASO is used */  
   /// main block  
   paso::SparseMatrix_ptr mainBlock;  
   /// coupling to neighbouring processors (row - col)  
   paso::SparseMatrix_ptr col_coupleBlock;  
   /// coupling to neighbouring processors (col - row)  
   paso::SparseMatrix_ptr row_coupleBlock;  
   /// coupling of rows-cols on neighbouring processors,  
   /// don't assume that this is set  
   paso::SparseMatrix_ptr remote_coupleBlock;  
   
   bool is_balanced;  
   double *balance_vector; /* matrix may be balanced by a diagonal matrix D=diagonal(balance_vector)  
                              if is_balanced is set, the matrix stored is D*A*D where A is the original matrix.  
                              When the system of linear equations is solved we solve D*A*D*y=c.  
                              So to solve A*x=b one needs to set c=D*b and x=D*y. */  
   
   index_t *global_id; /* store the global ids for all cols in col_couplerBlock */  
   
     
     
   index_t solver_package;  /* package controlling the solver pointer */  
   void* solver_p;  /* pointer to data needed by a solver */  
   
   /* this is only used for a trilinos matrix */  
   void *trilinos_data;  
 };  
   
 /*  interfaces: */  
   
 PASO_DLL_API  
 Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType, paso::SystemMatrixPattern_ptr, dim_t, dim_t, bool patternIsUnrolled);  
   
 PASO_DLL_API  
 Paso_SystemMatrix* Paso_SystemMatrix_getReference(Paso_SystemMatrix*);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_free(Paso_SystemMatrix*);  
   
   
 PASO_DLL_API  
 void Paso_SystemMatrix_MatrixVector(const double alpha, Paso_SystemMatrix* A, const double* in, const double beta, double* out);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_MatrixVector_CSR_OFFSET0(double alpha, Paso_SystemMatrix* A, const double* in, const double beta, double* out);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_nullifyRowsAndCols(Paso_SystemMatrix* A, double* mask_row, double* mask_col, double main_diagonal_value);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_applyBalanceInPlace(const Paso_SystemMatrix* A, double* x, const bool RHS);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_applyBalance(const Paso_SystemMatrix* A, double* x_out, const double* x, const bool RHS);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_balance(Paso_SystemMatrix* A);  
   
   
 PASO_DLL_API  
 void Paso_solve(Paso_SystemMatrix* A, double* out, double* in, Paso_Options* options);  
   
 PASO_DLL_API  
 void Paso_solve_free(Paso_SystemMatrix* in);  
   
 PASO_DLL_API  
 void  Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A,const double* in);  
   
 PASO_DLL_API  
 double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A);  
   
 PASO_DLL_API  
 void  Paso_SystemMatrix_startColCollect(Paso_SystemMatrix* A,const double* in);  
   
 PASO_DLL_API  
 double* Paso_SystemMatrix_finishColCollect(Paso_SystemMatrix* A);  
38    
39  PASO_DLL_API  namespace paso {
 void  Paso_SystemMatrix_startRowCollect(Paso_SystemMatrix* A,const double* in);  
   
 PASO_DLL_API  
 double* Paso_SystemMatrix_finishRowCollect(Paso_SystemMatrix* A);  
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getTotalNumRows(const Paso_SystemMatrix* A);  
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getTotalNumCols(const Paso_SystemMatrix*);  
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getGlobalNumRows(const Paso_SystemMatrix*);  
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getGlobalNumCols(const Paso_SystemMatrix*);  
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getGlobalTotalNumRows(const Paso_SystemMatrix* A);  
40    
41  PASO_DLL_API  struct SystemMatrix;
42  dim_t Paso_SystemMatrix_getGlobalTotalNumCols(const Paso_SystemMatrix* A);  typedef boost::shared_ptr<SystemMatrix> SystemMatrix_ptr;
43    typedef boost::shared_ptr<const SystemMatrix> const_SystemMatrix_ptr;
44    
45  PASO_DLL_API  typedef int SystemMatrixType;
 double Paso_SystemMatrix_getGlobalSize(const Paso_SystemMatrix*A);  
46    
47    //  this struct holds a (distributed) stiffness matrix
48  PASO_DLL_API  PASO_DLL_API
49  double Paso_SystemMatrix_getSparsity(const Paso_SystemMatrix*A);  struct SystemMatrix : boost::enable_shared_from_this<SystemMatrix>
50    {
51  PASO_DLL_API      SystemMatrix(SystemMatrixType, SystemMatrixPattern_ptr, dim_t, dim_t,
52  dim_t Paso_SystemMatrix_getNumRows(const Paso_SystemMatrix* A);                   bool patternIsUnrolled);
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getNumCols(const Paso_SystemMatrix* A);  
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getRowOverlap(const Paso_SystemMatrix* A);  
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getColOverlap(const Paso_SystemMatrix* A);  
   
   
   
   
 PASO_DLL_API  
 void Paso_SystemMatrix_saveMM(Paso_SystemMatrix *, char *);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_saveHB(Paso_SystemMatrix *, char *);  
   
 PASO_DLL_API  
 Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSR(char *);  
   
 PASO_DLL_API  
 Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSC(char *);  
   
 PASO_DLL_API  
 void Paso_RHS_loadMM_toCSR( char *fileName_p, double *b, dim_t size);  
   
   
 PASO_DLL_API  
 int Paso_SystemMatrix_getSystemMatrixTypeId(const index_t solver,const index_t preconditioner, const  index_t package,const  bool symmetry, Esys_MPIInfo *mpi_info);  
   
 PASO_DLL_API  
 dim_t Paso_SystemMatrix_getNumOutput(Paso_SystemMatrix* A);  
   
   
 PASO_DLL_API  
 void Paso_SystemMatrix_setValues(Paso_SystemMatrix*,double);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_add(Paso_SystemMatrix*,dim_t,index_t*, dim_t,dim_t,index_t*,dim_t, double*);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_rowSum(Paso_SystemMatrix* A, double* row_sum);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_nullifyRows(Paso_SystemMatrix* A, double* mask_row, double main_diagonal_value);  
   
   
 PASO_DLL_API  
 void Paso_SystemMatrix_makeZeroRowSums(Paso_SystemMatrix * A_p, double* left_over);  
   
   
 PASO_DLL_API  
 void Paso_SystemMatrix_copyBlockFromMainDiagonal(Paso_SystemMatrix * A_p, double* out);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_copyBlockToMainDiagonal(Paso_SystemMatrix * A_p, const double* in);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_copyFromMainDiagonal(Paso_SystemMatrix * A_p, double* out);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_copyToMainDiagonal(Paso_SystemMatrix * A_p, const double* in);  
   
   
 PASO_DLL_API  
 void Paso_SystemMatrix_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options);  
53    
54  PASO_DLL_API      ~SystemMatrix();
 void Paso_SystemMatrix_freePreconditioner(Paso_SystemMatrix* A);  
55    
56  PASO_DLL_API      /// Nullifies rows and columns in the matrix.
57  void Paso_SystemMatrix_copyColCoupleBlock(Paso_SystemMatrix *A);      /// The rows and columns are marked by positive values in mask_row and
58        /// mask_col. Values on the main diagonal which are marked to set to
59        /// zero by both mask_row and mask_col are set to main_diagonal_value.
60        void nullifyRowsAndCols(double* mask_row, double* mask_col,
61                                double main_diagonal_value);
62    
63        /// Nullifies rows in the matrix.
64        /// The rows are marked by positive values in mask_row. Values on the
65        /// main diagonal which are marked to set to zero by mask_row are set
66        /// to main_diagonal_value.
67        void nullifyRows(double* mask_row, double main_diagonal_value);
68    
69        void add(dim_t, index_t*, dim_t, dim_t, index_t*, dim_t, double*);
70    
71        void makeZeroRowSums(double* left_over);
72    
73        /// copies the col_coupleBlock into row_coupleBlock.
74        /// WARNING: this method uses mpi_requests of the coupler attached to the
75        /// matrix. No reordering on the received columns is performed.
76        /// In practice this means that components in
77        /// row_coupleBlock->pattern->index  and
78        /// row_coupler->connector->recv->shared
79        /// are ordered by increasing value.
80        /// Note that send and receive row_coupler->connectors are swapping roles.
81        void copyColCoupleBlock();
82    
83        void copyRemoteCoupleBlock(bool recreatePattern);
84    
85        void fillWithGlobalCoordinates(double f1);
86    
87        void print() const;
88    
89        void mergeMainAndCouple(index_t** p_ptr, index_t** p_idx, double** p_val);
90    
91        void mergeMainAndCouple_CSR_OFFSET0(index_t** p_ptr, index_t** p_idx, double** p_val);
92        void mergeMainAndCouple_CSR_OFFSET0_Block(index_t** p_ptr, index_t** p_idx, double** p_val);
93    
94        void mergeMainAndCouple_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val);
95    
96        void copyMain_CSC_OFFSET1(index_t** p_ptr, index_t** p_idx, double** p_val);
97    
98        void extendedRowsForST(dim_t* degree_ST, index_t* offset_ST, index_t* ST);
99    
100        void applyBalanceInPlace(double* x, bool RHS) const;
101    
102        void applyBalance(double* x_out, const double* x, bool RHS) const;
103    
104        void balance();
105    
106        double getGlobalSize() const;
107    
108        void setPreconditioner(Paso_Options* options);
109    
110        /// Applies the preconditioner.
111        /// This method needs to be called within a parallel region.
112        /// Barrier synchronization is performed before the evaluation to make
113        /// sure that the input vector is available
114        void solvePreconditioner(double* x, double* b);
115    
116        void freePreconditioner();
117    
118        index_t* borrowMainDiagonalPointer() const;
119    
120        inline void startCollect(const double* in)
121        {
122            startColCollect(in);
123        }
124    
125        inline double* finishCollect()
126        {
127            return finishColCollect();
128        }
129    
130        inline void startColCollect(const double* in)
131        {
132            col_coupler->startCollect(in);
133        }
134    
135        inline double* finishColCollect()
136        {
137            return col_coupler->finishCollect();
138        }
139    
140        inline void startRowCollect(const double* in)
141        {
142            row_coupler->startCollect(in);
143        }
144    
145        inline double* finishRowCollect()
146        {
147            return row_coupler->finishCollect();
148        }
149    
150        inline dim_t getNumRows() const
151        {
152            return mainBlock->numRows;
153        }
154    
155        inline dim_t getNumCols() const
156        {
157            return mainBlock->numCols;
158        }
159    
160        inline dim_t getTotalNumRows() const
161        {
162            return getNumRows() * row_block_size;
163        }
164    
165        inline dim_t getTotalNumCols() const
166        {
167            return getNumCols() * col_block_size;
168        }
169    
170        inline dim_t getRowOverlap()  const
171        {
172            return row_coupler->getNumOverlapComponents();
173        }
174    
175        inline dim_t getColOverlap() const
176        {
177            return col_coupler->getNumOverlapComponents();
178        }
179    
180        inline dim_t getGlobalNumRows() const
181        {
182            if (type & MATRIX_FORMAT_CSC) {
183                return pattern->input_distribution->getGlobalNumComponents();
184            }
185            return pattern->output_distribution->getGlobalNumComponents();
186        }
187    
188        inline dim_t getGlobalNumCols() const
189        {
190            if (type & MATRIX_FORMAT_CSC) {
191                return pattern->output_distribution->getGlobalNumComponents();
192            }
193            return pattern->input_distribution->getGlobalNumComponents();
194        }
195    
196        inline dim_t getGlobalTotalNumRows() const
197        {
198            return getGlobalNumRows() * row_block_size;
199        }
200    
201        inline dim_t getGlobalTotalNumCols() const
202        {
203            return getGlobalNumCols() * col_block_size;
204        }
205    
206        inline double getSparsity() const
207        {
208            return getGlobalSize() /
209                     (DBLE(getGlobalTotalNumRows())*getGlobalTotalNumCols());
210        }
211    
212        inline dim_t getNumOutput() const
213        {
214           return pattern->getNumOutput();
215        }
216    
217        inline void copyBlockFromMainDiagonal(double* out) const
218        {
219            mainBlock->copyBlockFromMainDiagonal(out);
220        }
221    
222        inline void copyBlockToMainDiagonal(const double* in)
223        {
224            mainBlock->copyBlockToMainDiagonal(in);
225        }
226    
227        inline void copyFromMainDiagonal(double* out) const
228        {
229            mainBlock->copyFromMainDiagonal(out);
230        }
231    
232        inline void copyToMainDiagonal(const double* in)
233        {
234            mainBlock->copyToMainDiagonal(in);
235        }
236    
237        inline void setValues(double value)
238        {
239            mainBlock->setValues(value);
240            col_coupleBlock->setValues(value);
241            row_coupleBlock->setValues(value);
242            is_balanced = false;
243        }
244    
245        inline void saveMM(const char* filename) const
246        {
247            if (mpi_info->size > 1) {
248                Esys_setError(IO_ERROR, "SystemMatrix::saveMM: Only single rank supported.");
249            } else {
250                mainBlock->saveMM(filename);
251            }
252        }
253    
254        inline void saveHB(const char *filename) const
255        {
256            if (mpi_info->size > 1) {
257                Esys_setError(TYPE_ERROR, "SystemMatrix::saveHB: Only single rank supported.");
258            } else if (!(type & MATRIX_FORMAT_CSC)) {
259                Esys_setError(TYPE_ERROR, "SystemMatrix::saveHB: Only CSC format supported.");
260            } else {
261                mainBlock->saveHB_CSC(filename);
262            }
263        }
264    
265        inline void rowSum(double* row_sum) const
266        {
267            if ((type & MATRIX_FORMAT_CSC) || (type & MATRIX_FORMAT_OFFSET1)) {
268                Esys_setError(TYPE_ERROR, "SystemMatrix::rowSum: No normalization "
269                      "available for compressed sparse column or index offset 1.");
270            } else {
271                const dim_t nrow = mainBlock->numRows*row_block_size;
272    #pragma omp parallel for
273                for (index_t irow=0; irow<nrow; ++irow) {
274                    row_sum[irow]=0.;
275                }
276                mainBlock->addRow_CSR_OFFSET0(row_sum);
277                col_coupleBlock->addRow_CSR_OFFSET0(row_sum);
278            }
279        }
280    
281        static SystemMatrix_ptr loadMM_toCSR(const char* filename);
282    
283        static SystemMatrix_ptr loadMM_toCSC(const char* filename);
284    
285        static index_t getSystemMatrixTypeId(index_t solver,
286                                             index_t preconditioner,
287                                             index_t package, bool symmetry,
288                                             Esys_MPIInfo* mpi_info);
289    
290        SystemMatrixType type;
291        SystemMatrixPattern_ptr pattern;
292    
293        dim_t logical_row_block_size;
294        dim_t logical_col_block_size;
295    
296        dim_t row_block_size;
297        dim_t col_block_size;
298        dim_t block_size;
299    
300        Distribution_ptr row_distribution;
301        Distribution_ptr col_distribution;
302        Esys_MPIInfo *mpi_info;
303    
304        Coupler_ptr col_coupler;
305        Coupler_ptr row_coupler;
306    
307        /// main block
308        SparseMatrix_ptr mainBlock;
309        /// coupling to neighbouring processors (row - col)
310        SparseMatrix_ptr col_coupleBlock;
311        /// coupling to neighbouring processors (col - row)
312        SparseMatrix_ptr row_coupleBlock;
313        /// coupling of rows-cols on neighbouring processors (may not be valid)
314        SparseMatrix_ptr remote_coupleBlock;
315    
316        bool is_balanced;
317    
318        /// matrix may be balanced by a diagonal matrix D=diagonal(balance_vector)
319        /// if is_balanced is true, the matrix stored is D*A*D where A is the
320        /// original matrix.
321        /// When the system of linear equations is solved we solve D*A*D*y=c.
322        /// So to solve A*x=b one needs to set c=D*b and x=D*y.
323        double* balance_vector;
324    
325        /// stores the global ids for all cols in col_coupleBlock
326        index_t* global_id;
327    
328  PASO_DLL_API      /// package code controlling the solver pointer
329  void Paso_SystemMatrix_copyRemoteCoupleBlock(Paso_SystemMatrix *A, bool recreatePattern);      index_t solver_package;
330    
331  PASO_DLL_API      /// pointer to data needed by a solver
332  void Paso_SystemMatrix_fillWithGlobalCoordinates(Paso_SystemMatrix *A, const double f1);      void* solver_p;
333    
334  PASO_DLL_API      /// this is only used for a trilinos matrix
335  void Paso_SystemMatrix_print(Paso_SystemMatrix *A);      void* trilinos_data;
336    };
337    
338    
339  PASO_DLL_API  void SystemMatrix_MatrixVector(double alpha, SystemMatrix_ptr A, const double* in, double beta, double* out);
 void Paso_SystemMatrix_mergeMainAndCouple(Paso_SystemMatrix *A, index_t **p_ptr, index_t **p_idx, double **p_val);  
340    
341  PASO_DLL_API  void SystemMatrix_MatrixVector_CSR_OFFSET0(double alpha, SystemMatrix_ptr A, const double* in, double beta, double* out);
 void Paso_SystemMatrix_mergeMainAndCouple_CSR_OFFSET0(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val);  
 void Paso_SystemMatrix_mergeMainAndCouple_CSR_OFFSET0_Block(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val);  
342    
343  PASO_DLL_API  void Paso_RHS_loadMM_toCSR(const char* filename, double* b, dim_t size);
 void Paso_SystemMatrix_mergeMainAndCouple_CSC_OFFSET1(Paso_SystemMatrix *A, index_t **p_ptr, index_t **p_idx, double **p_val);  
   
 PASO_DLL_API  
 void Paso_SystemMatrix_copyMain_CSC_OFFSET1(Paso_SystemMatrix* A, index_t** p_ptr, index_t** p_idx, double** p_val);  
344    
 void Paso_SystemMatrix_extendedRowsForST(Paso_SystemMatrix* A, dim_t* degree_ST, index_t* offset_ST, index_t* ST);  
345    
346    } // namespace paso
347        
348  #endif /* #ifndef INC_PASO_SYSTEMMATRIX */  #endif // __PASO_SYSTEMMATRIX_H__
349    

Legend:
Removed from v.4835  
changed lines
  Added in v.4836

  ViewVC Help
Powered by ViewVC 1.1.26