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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4836 - (show annotations)
Mon Apr 7 05:51:55 2014 UTC (6 years, 1 month ago) by caltinay
File MIME type: text/plain
File size: 10451 byte(s)
"Some" SystemMatrix clean up.....

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2014 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************/
19
20 /* Paso: SystemMatrix */
21
22 /****************************************************************************/
23
24 /* Copyrights by ACcESS Australia 2003,2004,2005,2006 */
25 /* Author: Lutz Gross, l.gross@uq.edu.au */
26
27 /****************************************************************************/
28
29 #ifndef __PASO_SYSTEMMATRIX_H__
30 #define __PASO_SYSTEMMATRIX_H__
31
32 #include "Common.h"
33 #include "Coupler.h"
34 #include "SparseMatrix.h"
35 #include "SystemMatrixPattern.h"
36 #include "Options.h"
37 #include "esysUtils/Esys_MPI.h"
38
39 namespace paso {
40
41 struct SystemMatrix;
42 typedef boost::shared_ptr<SystemMatrix> SystemMatrix_ptr;
43 typedef boost::shared_ptr<const SystemMatrix> const_SystemMatrix_ptr;
44
45 typedef int SystemMatrixType;
46
47 // this struct holds a (distributed) stiffness matrix
48 PASO_DLL_API
49 struct SystemMatrix : boost::enable_shared_from_this<SystemMatrix>
50 {
51 SystemMatrix(SystemMatrixType, SystemMatrixPattern_ptr, dim_t, dim_t,
52 bool patternIsUnrolled);
53
54 ~SystemMatrix();
55
56 /// Nullifies rows and columns in the matrix.
57 /// 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 /// package code controlling the solver pointer
329 index_t solver_package;
330
331 /// pointer to data needed by a solver
332 void* solver_p;
333
334 /// this is only used for a trilinos matrix
335 void* trilinos_data;
336 };
337
338
339 void SystemMatrix_MatrixVector(double alpha, SystemMatrix_ptr A, const double* in, double beta, double* out);
340
341 void SystemMatrix_MatrixVector_CSR_OFFSET0(double alpha, SystemMatrix_ptr A, const double* in, double beta, double* out);
342
343 void Paso_RHS_loadMM_toCSR(const char* filename, double* b, dim_t size);
344
345
346 } // namespace paso
347
348 #endif // __PASO_SYSTEMMATRIX_H__
349

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26