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

Annotation of /trunk/paso/src/SystemMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2551 - (hide annotations)
Thu Jul 23 09:19:15 2009 UTC (10 years, 1 month ago) by gross
File MIME type: text/plain
File size: 9316 byte(s)
a problem with the sparse matrix unrolling fixed.
1 ksteube 1312
2     /*******************************************************
3 ksteube 1811 *
4 jfenwick 2548 * Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1811 * 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 dhawcroft 631
14 ksteube 1811
15 jgs 150 /**************************************************************/
16    
17     /* Paso: SystemMatrix */
18    
19     /**************************************************************/
20    
21     /* Author: gross@access.edu.au */
22    
23     /**************************************************************/
24    
25     #include "SystemMatrix.h"
26    
27     /**************************************************************/
28    
29     /* allocates a SystemMatrix of type type using the given matrix pattern
30 gross 2551 Values are initialized by zero.
31     if patternIsUnrolled and type & MATRIX_FORMAT_BLK1, it is assumed that the pattern is allready unrolled to match the requested block size
32     and offsets otherwise unrolling and offset adjustment will be performed.
33     */
34 jgs 150
35 gross 2551 Paso_SystemMatrix* Paso_SystemMatrix_alloc(Paso_SystemMatrixType type,Paso_SystemMatrixPattern *pattern, int row_block_size, int col_block_size,
36     const bool_t patternIsUnrolled) {
37 gross 432
38 jgs 150 Paso_SystemMatrix*out=NULL;
39 gross 1028 dim_t n_norm,i;
40 ksteube 1312 Paso_SystemMatrixType pattern_format_out;
41 gross 2551 bool_t unroll=FALSE;
42 gross 1028
43 gross 2551 pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT;
44 jgs 150 Paso_resetError();
45 gross 2551 if (type & MATRIX_FORMAT_SYM) {
46     Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_alloc: Symmetric matrix patterns are not supported.");
47     return NULL;
48     }
49     if (patternIsUnrolled) {
50     if ( (type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & PATTERN_FORMAT_OFFSET1) ) {
51     Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_alloc: requested offset and pattern offset does not match.");
52     return NULL;
53     }
54     }
55     /* do we need to apply unrolling ? */
56     unroll
57     /* we don't like non-square blocks */
58     = (row_block_size!=col_block_size)
59     /* or any block size bigger than 3 */
60     || (col_block_size>3)
61     /* or if lock size one requested and the block size is not 1 */
62     || ((type & MATRIX_FORMAT_BLK1) && (col_block_size>1) ) ;
63    
64 jgs 150 out=MEMALLOC(1,Paso_SystemMatrix);
65     if (! Paso_checkPtr(out)) {
66 ksteube 1312 out->type=type;
67 jgs 150 out->pattern=NULL;
68 ksteube 1312 out->row_distribution=NULL;
69     out->col_distribution=NULL;
70     out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
71 gross 1552 out->row_coupler=NULL;
72     out->col_coupler=NULL;
73 ksteube 1312 out->mainBlock=NULL;
74 gross 1552 out->row_coupleBlock=NULL;
75     out->col_coupleBlock=NULL;
76 ksteube 1312 out->normalizer_is_valid=FALSE;
77     out->normalizer=NULL;
78 gross 425 out->solver_package=PASO_PASO;
79     out->solver=NULL;
80 ksteube 1312 out->trilinos_data=NULL;
81 jgs 150 out->reference_counter=1;
82 gross 2551 out->logical_row_block_size=row_block_size;
83     out->logical_col_block_size=col_block_size;
84 ksteube 1312
85 gross 1552
86 gross 415 if (type & MATRIX_FORMAT_CSC) {
87 gross 2551 if (unroll) {
88     if (patternIsUnrolled) {
89     out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
90     } else {
91     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
92     }
93     out->row_block_size=1;
94     out->col_block_size=1;
95     } else {
96 ksteube 1312 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
97 jgs 150 out->row_block_size=row_block_size;
98     out->col_block_size=col_block_size;
99 gross 2551 }
100     if (Paso_noError()) {
101     out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
102     out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
103     }
104 gross 415 } else {
105 gross 2551 if (unroll) {
106     if (patternIsUnrolled) {
107     out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
108     } else {
109     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
110     }
111 gross 415 out->row_block_size=1;
112     out->col_block_size=1;
113 gross 2551 } else {
114 ksteube 1312 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
115 gross 415 out->row_block_size=row_block_size;
116     out->col_block_size=col_block_size;
117 gross 2551 }
118     if (Paso_noError()) {
119     out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
120     out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
121     }
122     }
123     if (Paso_noError()) {
124     out->block_size=out->row_block_size*out->col_block_size;
125     out->col_coupler=Paso_Coupler_alloc(pattern->col_connector,out->col_block_size);
126     out->row_coupler=Paso_Coupler_alloc(pattern->row_connector,out->row_block_size);
127     /* this should be bypassed if trilinos is used */
128     if (type & MATRIX_FORMAT_TRILINOS_CRS) {
129     #ifdef TRILINOS
130     out->trilinos_data=Paso_TRILINOS_alloc();
131     #endif
132     } else {
133     out->solver_package=PASO_PASO;
134     out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size,TRUE);
135     out->col_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->col_couplePattern,row_block_size,col_block_size,TRUE);
136     out->row_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->row_couplePattern,row_block_size,col_block_size,TRUE);
137     if (Paso_noError()) {
138     /* allocate memory for matrix entries */
139     if (type & MATRIX_FORMAT_CSC) {
140     n_norm = out->mainBlock->numCols * out->col_block_size;
141     } else {
142     n_norm = out->mainBlock->numRows * out->row_block_size;
143     }
144     out->normalizer=MEMALLOC(n_norm,double);
145     out->normalizer_is_valid=FALSE;
146     if (! Paso_checkPtr(out->normalizer)) {
147     #pragma omp parallel for private(i) schedule(static)
148     for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
149     }
150 gross 415 }
151     }
152 jgs 150 }
153     }
154     /* all done: */
155     if (! Paso_noError()) {
156 ksteube 1312 Paso_SystemMatrix_free(out);
157 jgs 150 return NULL;
158     } else {
159     return out;
160     }
161     }
162    
163     /* returns a reference to Paso_SystemMatrix in */
164    
165 gross 2551 Paso_SystemMatrix* Paso_SystemMatrix_getReference(Paso_SystemMatrix* in) {
166 jgs 150 if (in!=NULL) ++(in->reference_counter);
167 gross 1792 return in;
168 jgs 150 }
169    
170     /* deallocates a SystemMatrix: */
171    
172 ksteube 1312 void Paso_SystemMatrix_free(Paso_SystemMatrix* in) {
173 jgs 150 if (in!=NULL) {
174     in->reference_counter--;
175     if (in->reference_counter<=0) {
176 ksteube 1312 Paso_SystemMatrixPattern_free(in->pattern);
177     Paso_Distribution_free(in->row_distribution);
178     Paso_Distribution_free(in->col_distribution);
179     Paso_MPIInfo_free(in->mpi_info);
180 gross 1552 Paso_Coupler_free(in->row_coupler);
181     Paso_Coupler_free(in->col_coupler);
182 ksteube 1312 Paso_SparseMatrix_free(in->mainBlock);
183 gross 1552 Paso_SparseMatrix_free(in->col_coupleBlock);
184     Paso_SparseMatrix_free(in->row_coupleBlock);
185 jgs 150 MEMFREE(in->normalizer);
186     Paso_solve_free(in);
187 ksteube 1312 #ifdef TRILINOS
188     Paso_TRILINOS_free(in->trilinos_data);
189     #endif
190 jgs 150 MEMFREE(in);
191     #ifdef Paso_TRACE
192 ksteube 1312 printf("Paso_SystemMatrix_free: system matrix as been deallocated.\n");
193 jgs 150 #endif
194     }
195     }
196     }
197 gross 1639 void Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A,const double* in)
198 ksteube 1312 {
199 gross 1552 Paso_SystemMatrix_startColCollect(A,in);
200 ksteube 1312 }
201     double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)
202     {
203 gross 1552 return Paso_SystemMatrix_finishColCollect(A);
204 ksteube 1312 }
205    
206 gross 1639 void Paso_SystemMatrix_startColCollect(Paso_SystemMatrix* A,const double* in)
207 gross 1552 {
208     Paso_Coupler_startCollect(A->col_coupler, in);
209     }
210     double* Paso_SystemMatrix_finishColCollect(Paso_SystemMatrix* A)
211     {
212     Paso_Coupler_finishCollect(A->col_coupler);
213     return A->col_coupler->recv_buffer;
214     }
215 gross 1639 void Paso_SystemMatrix_startRowCollect(Paso_SystemMatrix* A,const double* in)
216 gross 1552 {
217     Paso_Coupler_startCollect(A->row_coupler, in);
218     }
219     double* Paso_SystemMatrix_finishRowCollect(Paso_SystemMatrix* A)
220     {
221     Paso_Coupler_finishCollect(A->row_coupler);
222     return A->row_coupler->recv_buffer;
223     }
224    
225 gross 1407 dim_t Paso_SystemMatrix_getTotalNumRows(const Paso_SystemMatrix* A){
226 ksteube 1312 return A->mainBlock->numRows * A->row_block_size;
227     }
228    
229 gross 1407 dim_t Paso_SystemMatrix_getTotalNumCols(const Paso_SystemMatrix* A){
230 ksteube 1312 return A->mainBlock->numCols * A->col_block_size;
231     }
232     dim_t Paso_SystemMatrix_getGlobalNumRows(Paso_SystemMatrix* A) {
233     if (A->type & MATRIX_FORMAT_CSC) {
234     return Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
235     } else {
236     return Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
237     }
238     }
239     dim_t Paso_SystemMatrix_getGlobalNumCols(Paso_SystemMatrix* A) {
240     if (A->type & MATRIX_FORMAT_CSC) {
241     return Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
242     } else {
243     return Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
244     }
245    
246     }
247 gross 1361 dim_t Paso_SystemMatrix_getNumOutput(Paso_SystemMatrix* A) {
248     return Paso_SystemMatrixPattern_getNumOutput(A->pattern);
249     }
250 ksteube 1312

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26