/[escript]/branches/symbolic_from_3470/paso/src/SystemMatrix.c
ViewVC logotype

Contents of /branches/symbolic_from_3470/paso/src/SystemMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3868 - (show annotations)
Thu Mar 15 06:07:08 2012 UTC (7 years, 1 month ago) by caltinay
File MIME type: text/plain
File size: 13824 byte(s)
Update to latest trunk

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2010 by University of Queensland
5 * Earth Systems Science Computational Center (ESSCC)
6 * http://www.uq.edu.au/esscc
7 *
8 * Primary Business: Queensland, Australia
9 * Licensed under the Open Software License version 3.0
10 * http://www.opensource.org/licenses/osl-3.0.php
11 *
12 *******************************************************/
13
14
15 /**************************************************************/
16
17 /* Paso: SystemMatrix */
18
19 /**************************************************************/
20
21 /* Author: Lutz Gross, l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25 #include "SystemMatrix.h"
26 #include "Preconditioner.h"
27
28 /**************************************************************/
29
30 /* Allocates a SystemMatrix of given type using the given matrix pattern.
31 Values are initialized with zero.
32 If patternIsUnrolled and type & MATRIX_FORMAT_BLK1, it is assumed
33 that the pattern is already unrolled to match the requested block size
34 and offsets. Otherwise unrolling and offset adjustment will be performed. */
35
36 Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType type,Paso_SystemMatrixPattern *pattern, int row_block_size, int col_block_size,
37 const bool_t patternIsUnrolled) {
38
39 Paso_SystemMatrix*out=NULL;
40 dim_t n_norm,i;
41 Paso_SystemMatrixType pattern_format_out;
42 bool_t unroll=FALSE;
43
44 pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? MATRIX_FORMAT_OFFSET1: MATRIX_FORMAT_DEFAULT;
45 Esys_resetError();
46 if (patternIsUnrolled) {
47 if ( ! XNOR(type & MATRIX_FORMAT_OFFSET1, pattern->type & MATRIX_FORMAT_OFFSET1) ) {
48 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_alloc: requested offset and pattern offset do not match.");
49 return NULL;
50 }
51 }
52 /* do we need to apply unrolling? */
53 unroll
54 /* we don't like non-square blocks */
55 = (row_block_size!=col_block_size)
56 #ifndef USE_LAPACK
57 /* or any block size bigger than 3 */
58 || (col_block_size>3)
59 # endif
60 /* or if lock size one requested and the block size is not 1 */
61 || ((type & MATRIX_FORMAT_BLK1) && (col_block_size>1) )
62 /* or the offsets are wrong */
63 || ((type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & MATRIX_FORMAT_OFFSET1));
64
65 out=MEMALLOC(1,Paso_SystemMatrix);
66 if (! Esys_checkPtr(out)) {
67 out->type=type;
68 out->pattern=NULL;
69 out->row_distribution=NULL;
70 out->col_distribution=NULL;
71 out->mpi_info=Esys_MPIInfo_getReference(pattern->mpi_info);
72 out->row_coupler=NULL;
73 out->col_coupler=NULL;
74 out->mainBlock=NULL;
75 out->row_coupleBlock=NULL;
76 out->col_coupleBlock=NULL;
77 out->remote_coupleBlock=NULL;
78 out->is_balanced=FALSE;
79 out->balance_vector=NULL;
80 out->global_id=NULL;
81 out->solver_package=PASO_PASO;
82 out->solver_p=NULL;
83 out->trilinos_data=NULL;
84 out->reference_counter=1;
85 out->logical_row_block_size=row_block_size;
86 out->logical_col_block_size=col_block_size;
87
88
89 if (type & MATRIX_FORMAT_CSC) {
90 if (unroll) {
91 if (patternIsUnrolled) {
92 out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
93 } else {
94 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
95 }
96 out->row_block_size=1;
97 out->col_block_size=1;
98 } else {
99 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
100 out->row_block_size=row_block_size;
101 out->col_block_size=col_block_size;
102 }
103 if (Esys_noError()) {
104 out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
105 out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
106 }
107 } else {
108 if (unroll) {
109 if (patternIsUnrolled) {
110 out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
111 } else {
112 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
113 }
114 out->row_block_size=1;
115 out->col_block_size=1;
116 } else {
117 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
118 out->row_block_size=row_block_size;
119 out->col_block_size=col_block_size;
120 }
121 if (Esys_noError()) {
122 out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
123 out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
124 }
125 }
126 if (Esys_noError()) {
127 if (type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
128 out->block_size=MIN(out->row_block_size,out->col_block_size);
129 } else {
130 out->block_size=out->row_block_size*out->col_block_size;
131 }
132 out->col_coupler=Paso_Coupler_alloc(out->pattern->col_connector,out->col_block_size);
133 out->row_coupler=Paso_Coupler_alloc(out->pattern->row_connector,out->row_block_size);
134 /* this should be bypassed if trilinos is used */
135 if (type & MATRIX_FORMAT_TRILINOS_CRS) {
136 #ifdef TRILINOS
137 out->trilinos_data=Paso_TRILINOS_alloc();
138 #endif
139 } else {
140 out->solver_package=PASO_PASO;
141 out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size,TRUE);
142 out->col_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->col_couplePattern,row_block_size,col_block_size,TRUE);
143 out->row_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->row_couplePattern,row_block_size,col_block_size,TRUE);
144 if (Esys_noError()) {
145 /* allocate memory for matrix entries */
146 n_norm = MAX(out->mainBlock->numCols * out->col_block_size, out->mainBlock->numRows * out->row_block_size);
147 out->balance_vector=MEMALLOC(n_norm,double);
148 out->is_balanced=FALSE;
149 if (! Esys_checkPtr(out->balance_vector)) {
150 #pragma omp parallel for private(i) schedule(static)
151 for (i=0;i<n_norm;++i) out->balance_vector[i]=1.;
152 }
153 }
154 }
155 }
156 }
157 /* all done: */
158 if (! Esys_noError()) {
159 Paso_SystemMatrix_free(out);
160 return NULL;
161 } else {
162 return out;
163 }
164 }
165
166 /* returns a reference to Paso_SystemMatrix in */
167
168 Paso_SystemMatrix* Paso_SystemMatrix_getReference(Paso_SystemMatrix* in) {
169 if (in!=NULL) ++(in->reference_counter);
170 return in;
171 }
172
173 /* deallocates a SystemMatrix: */
174
175 void Paso_SystemMatrix_free(Paso_SystemMatrix* in) {
176 if (in!=NULL) {
177 in->reference_counter--;
178 if (in->reference_counter<=0) {
179
180 Paso_solve_free(in);
181 Paso_SystemMatrixPattern_free(in->pattern);
182 Paso_Distribution_free(in->row_distribution);
183 Paso_Distribution_free(in->col_distribution);
184 Esys_MPIInfo_free(in->mpi_info);
185 Paso_Coupler_free(in->row_coupler);
186 Paso_Coupler_free(in->col_coupler);
187 Paso_SparseMatrix_free(in->mainBlock);
188 Paso_SparseMatrix_free(in->col_coupleBlock);
189 Paso_SparseMatrix_free(in->row_coupleBlock);
190 Paso_SparseMatrix_free(in->remote_coupleBlock);
191
192 MEMFREE(in->balance_vector);
193 if (in->global_id) MEMFREE(in->global_id);
194 #ifdef TRILINOS
195 Paso_TRILINOS_free(in->trilinos_data);
196 #endif
197 MEMFREE(in);
198 #ifdef Paso_TRACE
199 printf("Paso_SystemMatrix_free: system matrix has been deallocated.\n");
200 #endif
201 }
202 }
203 }
204 void Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A,const double* in)
205 {
206 Paso_SystemMatrix_startColCollect(A,in);
207 }
208 double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)
209 {
210 return Paso_SystemMatrix_finishColCollect(A);
211 }
212
213 void Paso_SystemMatrix_startColCollect(Paso_SystemMatrix* A,const double* in)
214 {
215 Paso_Coupler_startCollect(A->col_coupler, in);
216 }
217 double* Paso_SystemMatrix_finishColCollect(Paso_SystemMatrix* A)
218 {
219 Paso_Coupler_finishCollect(A->col_coupler);
220 return A->col_coupler->recv_buffer;
221 }
222 void Paso_SystemMatrix_startRowCollect(Paso_SystemMatrix* A,const double* in)
223 {
224 Paso_Coupler_startCollect(A->row_coupler, in);
225 }
226 double* Paso_SystemMatrix_finishRowCollect(Paso_SystemMatrix* A)
227 {
228 Paso_Coupler_finishCollect(A->row_coupler);
229 return A->row_coupler->recv_buffer;
230 }
231
232
233 dim_t Paso_SystemMatrix_getNumRows(const Paso_SystemMatrix* A){
234 return A->mainBlock->numRows;
235 }
236
237 dim_t Paso_SystemMatrix_getNumCols(const Paso_SystemMatrix* A){
238 return A->mainBlock->numCols;
239 }
240
241 dim_t Paso_SystemMatrix_getTotalNumRows(const Paso_SystemMatrix* A){
242 return Paso_SystemMatrix_getNumRows(A) * A->row_block_size;
243 }
244 dim_t Paso_SystemMatrix_getTotalNumCols(const Paso_SystemMatrix* A){
245 return Paso_SystemMatrix_getNumCols(A) * A->col_block_size;
246 }
247
248 dim_t Paso_SystemMatrix_getRowOverlap(const Paso_SystemMatrix* A)
249 {
250 return Paso_Coupler_getNumOverlapComponents(A->row_coupler);
251 }
252 dim_t Paso_SystemMatrix_getColOverlap(const Paso_SystemMatrix* A)
253 {
254 return Paso_Coupler_getNumOverlapComponents(A->col_coupler);
255 }
256
257
258
259 dim_t Paso_SystemMatrix_getGlobalNumRows(const Paso_SystemMatrix* A) {
260 if (A->type & MATRIX_FORMAT_CSC) {
261 return Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
262 } else {
263 return Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
264 }
265 }
266 dim_t Paso_SystemMatrix_getGlobalNumCols(const Paso_SystemMatrix* A) {
267 if (A->type & MATRIX_FORMAT_CSC) {
268 return Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
269 } else {
270 return Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
271 }
272
273 }
274 dim_t Paso_SystemMatrix_getGlobalTotalNumRows(const Paso_SystemMatrix* A){
275 return Paso_SystemMatrix_getGlobalNumRows(A) * A->row_block_size;
276 }
277
278 dim_t Paso_SystemMatrix_getGlobalTotalNumCols(const Paso_SystemMatrix* A){
279 return Paso_SystemMatrix_getGlobalNumCols(A) * A->col_block_size;
280 }
281 double Paso_SystemMatrix_getGlobalSize(const Paso_SystemMatrix*A) {
282 double global_size=0;
283 if (A!=NULL) {
284 double my_size=Paso_SparseMatrix_getSize(A->mainBlock)+Paso_SparseMatrix_getSize(A->col_coupleBlock);
285 if (A->mpi_info->size > 1) {
286
287 #ifdef ESYS_MPI
288 MPI_Allreduce(&my_size,&global_size, 1, MPI_DOUBLE, MPI_SUM, A->mpi_info->comm);
289 #else
290 global_size=my_size;
291 #endif
292
293 } else {
294 global_size=my_size;
295 }
296 }
297 return global_size;
298 }
299 double Paso_SystemMatrix_getSparsity(const Paso_SystemMatrix*A) {
300 return Paso_SystemMatrix_getGlobalSize(A)/(DBLE(Paso_SystemMatrix_getGlobalTotalNumRows(A))*DBLE(Paso_SystemMatrix_getGlobalTotalNumCols(A)));
301 }
302
303 dim_t Paso_SystemMatrix_getNumOutput(Paso_SystemMatrix* A) {
304 return Paso_SystemMatrixPattern_getNumOutput(A->pattern);
305 }
306
307 index_t* Paso_SystemMatrix_borrowMainDiagonalPointer(Paso_SystemMatrix * A_p)
308 {
309 index_t* out=NULL;
310 int fail=0;
311 out=Paso_SparseMatrix_borrowMainDiagonalPointer(A_p->mainBlock);
312 if (out==NULL) fail=1;
313 #ifdef ESYS_MPI
314 {
315 int fail_loc = fail;
316 MPI_Allreduce(&fail_loc, &fail, 1, MPI_INT, MPI_MAX, A_p->mpi_info->comm);
317 }
318 #endif
319 if (fail>0) Esys_setError(VALUE_ERROR, "Paso_SystemMatrix_borrowMainDiagonalPointer: no main diagonal");
320 return out;
321 }
322
323 void Paso_SystemMatrix_makeZeroRowSums(Paso_SystemMatrix * A_p, double* left_over)
324 {
325 index_t ir, ib, irow;
326 register double rtmp1, rtmp2;
327 const dim_t n = Paso_SystemMatrixPattern_getNumOutput(A_p->pattern);
328 const dim_t nblk = A_p->block_size;
329 const dim_t blk = A_p->row_block_size;
330 const index_t* main_ptr=Paso_SystemMatrix_borrowMainDiagonalPointer(A_p);
331
332
333 Paso_SystemMatrix_rowSum(A_p, left_over); /* left_over now holds the row sum */
334
335 #pragma omp parallel for private(ir,ib, rtmp1, rtmp2) schedule(static)
336 for (ir=0;ir< n;ir++) {
337 for (ib=0;ib<blk; ib++) {
338 irow=ib+blk*ir;
339 rtmp1=left_over[irow];
340 rtmp2=A_p->mainBlock->val[main_ptr[ir]*nblk+ib+blk*ib];
341 A_p->mainBlock->val[main_ptr[ir]*nblk+ib+blk*ib] = -rtmp1;
342 left_over[irow]=rtmp2+rtmp1;
343 }
344 }
345 }
346 void Paso_SystemMatrix_copyBlockFromMainDiagonal(Paso_SystemMatrix * A_p, double* out)
347 {
348 Paso_SparseMatrix_copyBlockFromMainDiagonal(A_p->mainBlock, out);
349 return;
350 }
351 void Paso_SystemMatrix_copyBlockToMainDiagonal(Paso_SystemMatrix * A_p, const double* in)
352 {
353 Paso_SparseMatrix_copyBlockToMainDiagonal(A_p->mainBlock, in);
354 return;
355 }
356 void Paso_SystemMatrix_copyFromMainDiagonal(Paso_SystemMatrix * A_p, double* out)
357 {
358 Paso_SparseMatrix_copyFromMainDiagonal(A_p->mainBlock, out);
359 return;
360 }
361 void Paso_SystemMatrix_copyToMainDiagonal(Paso_SystemMatrix * A_p, const double* in)
362 {
363 Paso_SparseMatrix_copyToMainDiagonal(A_p->mainBlock, in);
364 return;
365 }
366
367 void Paso_SystemMatrix_setPreconditioner(Paso_SystemMatrix* A,Paso_Options* options) {
368 if (A->solver_p==NULL) {
369 A->solver_p=Paso_Preconditioner_alloc(A,options);
370 }
371 }
372
373 /* Applies the preconditioner. */
374 /* Has to be called within a parallel region. */
375 /* Barrier synchronization is performed before the evaluation to make sure
376 that the input vector is available */
377 void Paso_SystemMatrix_solvePreconditioner(Paso_SystemMatrix* A,double* x,double* b){
378 Paso_Preconditioner* prec=(Paso_Preconditioner*) A->solver_p;
379 Paso_Preconditioner_solve(prec, A,x,b);
380 }
381 void Paso_SystemMatrix_freePreconditioner(Paso_SystemMatrix* A) {
382 Paso_Preconditioner* prec=NULL;
383 if (A!=NULL) {
384 prec=(Paso_Preconditioner*) A->solver_p;
385 Paso_Preconditioner_free(prec);
386 A->solver_p=NULL;
387 }
388 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26