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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2554 - (hide annotations)
Fri Jul 24 05:38:54 2009 UTC (10 years, 2 months ago) by gross
File MIME type: text/plain
File size: 9441 byte(s)
bug with MKL call 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 gross 2554 if ( ! XNOR(type & MATRIX_FORMAT_OFFSET1, pattern->type & PATTERN_FORMAT_OFFSET1) ) {
51 gross 2551 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 gross 2554 || ((type & MATRIX_FORMAT_BLK1) && (col_block_size>1) )
63     /* or the offsets are wrong */
64     || ((type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & PATTERN_FORMAT_OFFSET1));
65 gross 2551
66 jgs 150 out=MEMALLOC(1,Paso_SystemMatrix);
67     if (! Paso_checkPtr(out)) {
68 ksteube 1312 out->type=type;
69 jgs 150 out->pattern=NULL;
70 ksteube 1312 out->row_distribution=NULL;
71     out->col_distribution=NULL;
72     out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
73 gross 1552 out->row_coupler=NULL;
74     out->col_coupler=NULL;
75 ksteube 1312 out->mainBlock=NULL;
76 gross 1552 out->row_coupleBlock=NULL;
77     out->col_coupleBlock=NULL;
78 ksteube 1312 out->normalizer_is_valid=FALSE;
79     out->normalizer=NULL;
80 gross 425 out->solver_package=PASO_PASO;
81     out->solver=NULL;
82 ksteube 1312 out->trilinos_data=NULL;
83 jgs 150 out->reference_counter=1;
84 gross 2551 out->logical_row_block_size=row_block_size;
85     out->logical_col_block_size=col_block_size;
86 ksteube 1312
87 gross 1552
88 gross 415 if (type & MATRIX_FORMAT_CSC) {
89 gross 2551 if (unroll) {
90     if (patternIsUnrolled) {
91     out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
92     } else {
93     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
94     }
95     out->row_block_size=1;
96     out->col_block_size=1;
97     } else {
98 ksteube 1312 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
99 jgs 150 out->row_block_size=row_block_size;
100     out->col_block_size=col_block_size;
101 gross 2551 }
102     if (Paso_noError()) {
103     out->row_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
104     out->col_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
105     }
106 gross 415 } else {
107 gross 2551 if (unroll) {
108     if (patternIsUnrolled) {
109     out->pattern=Paso_SystemMatrixPattern_getReference(pattern);
110     } else {
111     out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
112     }
113 gross 415 out->row_block_size=1;
114     out->col_block_size=1;
115 gross 2551 } else {
116 ksteube 1312 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
117 gross 415 out->row_block_size=row_block_size;
118     out->col_block_size=col_block_size;
119 gross 2551 }
120     if (Paso_noError()) {
121     out->row_distribution=Paso_Distribution_getReference(out->pattern->output_distribution);
122     out->col_distribution=Paso_Distribution_getReference(out->pattern->input_distribution);
123     }
124     }
125     if (Paso_noError()) {
126     out->block_size=out->row_block_size*out->col_block_size;
127     out->col_coupler=Paso_Coupler_alloc(pattern->col_connector,out->col_block_size);
128     out->row_coupler=Paso_Coupler_alloc(pattern->row_connector,out->row_block_size);
129     /* this should be bypassed if trilinos is used */
130     if (type & MATRIX_FORMAT_TRILINOS_CRS) {
131     #ifdef TRILINOS
132     out->trilinos_data=Paso_TRILINOS_alloc();
133     #endif
134     } else {
135     out->solver_package=PASO_PASO;
136     out->mainBlock=Paso_SparseMatrix_alloc(type,out->pattern->mainPattern,row_block_size,col_block_size,TRUE);
137     out->col_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->col_couplePattern,row_block_size,col_block_size,TRUE);
138     out->row_coupleBlock=Paso_SparseMatrix_alloc(type,out->pattern->row_couplePattern,row_block_size,col_block_size,TRUE);
139     if (Paso_noError()) {
140     /* allocate memory for matrix entries */
141     if (type & MATRIX_FORMAT_CSC) {
142     n_norm = out->mainBlock->numCols * out->col_block_size;
143     } else {
144     n_norm = out->mainBlock->numRows * out->row_block_size;
145     }
146     out->normalizer=MEMALLOC(n_norm,double);
147     out->normalizer_is_valid=FALSE;
148     if (! Paso_checkPtr(out->normalizer)) {
149     #pragma omp parallel for private(i) schedule(static)
150     for (i=0;i<n_norm;++i) out->normalizer[i]=0.;
151     }
152 gross 415 }
153     }
154 jgs 150 }
155     }
156     /* all done: */
157     if (! Paso_noError()) {
158 ksteube 1312 Paso_SystemMatrix_free(out);
159 jgs 150 return NULL;
160     } else {
161     return out;
162     }
163     }
164    
165     /* returns a reference to Paso_SystemMatrix in */
166    
167 gross 2551 Paso_SystemMatrix* Paso_SystemMatrix_getReference(Paso_SystemMatrix* in) {
168 jgs 150 if (in!=NULL) ++(in->reference_counter);
169 gross 1792 return in;
170 jgs 150 }
171    
172     /* deallocates a SystemMatrix: */
173    
174 ksteube 1312 void Paso_SystemMatrix_free(Paso_SystemMatrix* in) {
175 jgs 150 if (in!=NULL) {
176     in->reference_counter--;
177     if (in->reference_counter<=0) {
178 ksteube 1312 Paso_SystemMatrixPattern_free(in->pattern);
179     Paso_Distribution_free(in->row_distribution);
180     Paso_Distribution_free(in->col_distribution);
181     Paso_MPIInfo_free(in->mpi_info);
182 gross 1552 Paso_Coupler_free(in->row_coupler);
183     Paso_Coupler_free(in->col_coupler);
184 ksteube 1312 Paso_SparseMatrix_free(in->mainBlock);
185 gross 1552 Paso_SparseMatrix_free(in->col_coupleBlock);
186     Paso_SparseMatrix_free(in->row_coupleBlock);
187 jgs 150 MEMFREE(in->normalizer);
188     Paso_solve_free(in);
189 ksteube 1312 #ifdef TRILINOS
190     Paso_TRILINOS_free(in->trilinos_data);
191     #endif
192 jgs 150 MEMFREE(in);
193     #ifdef Paso_TRACE
194 ksteube 1312 printf("Paso_SystemMatrix_free: system matrix as been deallocated.\n");
195 jgs 150 #endif
196     }
197     }
198     }
199 gross 1639 void Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A,const double* in)
200 ksteube 1312 {
201 gross 1552 Paso_SystemMatrix_startColCollect(A,in);
202 ksteube 1312 }
203     double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)
204     {
205 gross 1552 return Paso_SystemMatrix_finishColCollect(A);
206 ksteube 1312 }
207    
208 gross 1639 void Paso_SystemMatrix_startColCollect(Paso_SystemMatrix* A,const double* in)
209 gross 1552 {
210     Paso_Coupler_startCollect(A->col_coupler, in);
211     }
212     double* Paso_SystemMatrix_finishColCollect(Paso_SystemMatrix* A)
213     {
214     Paso_Coupler_finishCollect(A->col_coupler);
215     return A->col_coupler->recv_buffer;
216     }
217 gross 1639 void Paso_SystemMatrix_startRowCollect(Paso_SystemMatrix* A,const double* in)
218 gross 1552 {
219     Paso_Coupler_startCollect(A->row_coupler, in);
220     }
221     double* Paso_SystemMatrix_finishRowCollect(Paso_SystemMatrix* A)
222     {
223     Paso_Coupler_finishCollect(A->row_coupler);
224     return A->row_coupler->recv_buffer;
225     }
226    
227 gross 1407 dim_t Paso_SystemMatrix_getTotalNumRows(const Paso_SystemMatrix* A){
228 ksteube 1312 return A->mainBlock->numRows * A->row_block_size;
229     }
230    
231 gross 1407 dim_t Paso_SystemMatrix_getTotalNumCols(const Paso_SystemMatrix* A){
232 ksteube 1312 return A->mainBlock->numCols * A->col_block_size;
233     }
234     dim_t Paso_SystemMatrix_getGlobalNumRows(Paso_SystemMatrix* A) {
235     if (A->type & MATRIX_FORMAT_CSC) {
236     return Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
237     } else {
238     return Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
239     }
240     }
241     dim_t Paso_SystemMatrix_getGlobalNumCols(Paso_SystemMatrix* A) {
242     if (A->type & MATRIX_FORMAT_CSC) {
243     return Paso_Distribution_getGlobalNumComponents(A->pattern->output_distribution);
244     } else {
245     return Paso_Distribution_getGlobalNumComponents(A->pattern->input_distribution);
246     }
247    
248     }
249 gross 1361 dim_t Paso_SystemMatrix_getNumOutput(Paso_SystemMatrix* A) {
250     return Paso_SystemMatrixPattern_getNumOutput(A->pattern);
251     }
252 ksteube 1312

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26