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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2608 - (show annotations)
Tue Aug 18 01:25:18 2009 UTC (10 years, 1 month ago) by jfenwick
File MIME type: text/plain
File size: 9451 byte(s)
Updating Lutz' email

1
2 /*******************************************************
3 *
4 * Copyright (c) 2003-2009 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
27 /**************************************************************/
28
29 /* allocates a SystemMatrix of type type using the given matrix pattern
30 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
35 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
38 Paso_SystemMatrix*out=NULL;
39 dim_t n_norm,i;
40 Paso_SystemMatrixType pattern_format_out;
41 bool_t unroll=FALSE;
42
43 pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT;
44 Paso_resetError();
45 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 ( ! XNOR(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 /* or the offsets are wrong */
64 || ((type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & PATTERN_FORMAT_OFFSET1));
65
66 out=MEMALLOC(1,Paso_SystemMatrix);
67 if (! Paso_checkPtr(out)) {
68 out->type=type;
69 out->pattern=NULL;
70 out->row_distribution=NULL;
71 out->col_distribution=NULL;
72 out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info);
73 out->row_coupler=NULL;
74 out->col_coupler=NULL;
75 out->mainBlock=NULL;
76 out->row_coupleBlock=NULL;
77 out->col_coupleBlock=NULL;
78 out->normalizer_is_valid=FALSE;
79 out->normalizer=NULL;
80 out->solver_package=PASO_PASO;
81 out->solver=NULL;
82 out->trilinos_data=NULL;
83 out->reference_counter=1;
84 out->logical_row_block_size=row_block_size;
85 out->logical_col_block_size=col_block_size;
86
87
88 if (type & MATRIX_FORMAT_CSC) {
89 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 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
99 out->row_block_size=row_block_size;
100 out->col_block_size=col_block_size;
101 }
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 } else {
107 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 out->row_block_size=1;
114 out->col_block_size=1;
115 } else {
116 out->pattern=Paso_SystemMatrixPattern_unrollBlocks(pattern,pattern_format_out,1,1);
117 out->row_block_size=row_block_size;
118 out->col_block_size=col_block_size;
119 }
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 }
153 }
154 }
155 }
156 /* all done: */
157 if (! Paso_noError()) {
158 Paso_SystemMatrix_free(out);
159 return NULL;
160 } else {
161 return out;
162 }
163 }
164
165 /* returns a reference to Paso_SystemMatrix in */
166
167 Paso_SystemMatrix* Paso_SystemMatrix_getReference(Paso_SystemMatrix* in) {
168 if (in!=NULL) ++(in->reference_counter);
169 return in;
170 }
171
172 /* deallocates a SystemMatrix: */
173
174 void Paso_SystemMatrix_free(Paso_SystemMatrix* in) {
175 if (in!=NULL) {
176 in->reference_counter--;
177 if (in->reference_counter<=0) {
178 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 Paso_Coupler_free(in->row_coupler);
183 Paso_Coupler_free(in->col_coupler);
184 Paso_SparseMatrix_free(in->mainBlock);
185 Paso_SparseMatrix_free(in->col_coupleBlock);
186 Paso_SparseMatrix_free(in->row_coupleBlock);
187 MEMFREE(in->normalizer);
188 Paso_solve_free(in);
189 #ifdef TRILINOS
190 Paso_TRILINOS_free(in->trilinos_data);
191 #endif
192 MEMFREE(in);
193 #ifdef Paso_TRACE
194 printf("Paso_SystemMatrix_free: system matrix as been deallocated.\n");
195 #endif
196 }
197 }
198 }
199 void Paso_SystemMatrix_startCollect(Paso_SystemMatrix* A,const double* in)
200 {
201 Paso_SystemMatrix_startColCollect(A,in);
202 }
203 double* Paso_SystemMatrix_finishCollect(Paso_SystemMatrix* A)
204 {
205 return Paso_SystemMatrix_finishColCollect(A);
206 }
207
208 void Paso_SystemMatrix_startColCollect(Paso_SystemMatrix* A,const double* in)
209 {
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 void Paso_SystemMatrix_startRowCollect(Paso_SystemMatrix* A,const double* in)
218 {
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 dim_t Paso_SystemMatrix_getTotalNumRows(const Paso_SystemMatrix* A){
228 return A->mainBlock->numRows * A->row_block_size;
229 }
230
231 dim_t Paso_SystemMatrix_getTotalNumCols(const Paso_SystemMatrix* A){
232 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 dim_t Paso_SystemMatrix_getNumOutput(Paso_SystemMatrix* A) {
250 return Paso_SystemMatrixPattern_getNumOutput(A->pattern);
251 }
252

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26