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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2551 - (show annotations)
Thu Jul 23 09:19:15 2009 UTC (11 years, 2 months ago) by gross
File MIME type: text/plain
File size: 12370 byte(s)
a problem with the sparse matrix unrolling fixed.
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: SparseMatrix */
18
19 /**************************************************************/
20
21 /* Author: gross@access.edu.au */
22
23 /**************************************************************/
24
25 #include "Paso.h"
26 #include "SparseMatrix.h"
27 #include "TRILINOS.h"
28 #include "mmio.h"
29
30 /**************************************************************/
31 static void swap( index_t*, index_t*, double*, int, int );
32 static void q_sort( index_t*, index_t*, double*, int, int );
33 /*static void print_entries( index_t*, index_t*, double* );*/
34
35 static int M, N, nz;
36
37
38 /* debug: print the entries */
39 /*
40 void print_entries( index_t *r, index_t *c, double *v )
41 {
42 int i;
43
44 for( i=0; i<nz; i++ )
45 {
46 printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] );
47 }
48 }
49 */
50
51 /* swap function */
52 void swap( index_t *r, index_t *c, double *v, int left, int right )
53 {
54 double v_temp;
55 index_t temp;
56
57 temp = r[left];
58 r[left] = r[right];
59 r[right] = temp;
60
61 temp = c[left];
62 c[left] = c[right];
63 c[right] = temp;
64
65 v_temp = v[left];
66 v[left] = v[right];
67 v[right] = v_temp;
68 }
69
70 void q_sort( index_t *row, index_t *col, double *val, int begin, int end )
71 {
72 int l, r;
73 index_t pivot, lval;
74
75 if( end > begin )
76 {
77 pivot = N * row[begin] + col[begin];
78 l = begin + 1;
79 r = end;
80
81 while( l < r )
82 {
83 lval = N * row[l] + col[l];
84 if( lval < pivot )
85 l++;
86 else
87 {
88 r--;
89 swap( row, col, val, l, r );
90 }
91 }
92 l--;
93 swap( row, col, val, begin, l );
94 q_sort( row, col, val, begin, l );
95 q_sort( row, col, val, r, end );
96 }
97 }
98
99
100 /* allocates a SparseMatrix of type type using the given matrix pattern
101 Values are initialized by zero.
102 if patternIsUnrolled and type & MATRIX_FORMAT_BLK1, it is assumed that the pattern is allready unrolled to match the requested block size
103 and offsets otherwise unrolling and offset adjustment will be performed.
104 */
105 Paso_SparseMatrix* Paso_SparseMatrix_alloc(Paso_SparseMatrixType type,Paso_Pattern *pattern, int row_block_size, int col_block_size, const bool_t patternIsUnrolled) {
106
107 Paso_SparseMatrix*out=NULL;
108 Paso_SparseMatrixType pattern_format_out;
109 bool_t unroll=FALSE;
110
111 if (type & MATRIX_FORMAT_SYM) {
112 Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_alloc: symmetric matrix pattern are not supported.");
113 return NULL;
114 }
115 if (patternIsUnrolled) {
116 if ( (type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & PATTERN_FORMAT_OFFSET1) ) {
117 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_alloc: requested offset and pattern offset does not match.");
118 return NULL;
119 }
120 }
121 /* do we need to apply unrolling ? */
122 unroll
123 /* we don't like non-square blocks */
124 = (row_block_size!=col_block_size)
125 /* or any block size bigger than 3 */
126 || (col_block_size>3)
127 /* or if lock size one requested and the block size is not 1 */
128 || ((type & MATRIX_FORMAT_BLK1) && (col_block_size>1) ) ;
129
130 pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT;
131
132 Paso_resetError();
133 out=MEMALLOC(1,Paso_SparseMatrix);
134 if (! Paso_checkPtr(out)) {
135 out->pattern=NULL;
136 out->val=NULL;
137 out->reference_counter=1;
138 out->type=type;
139 out->solver=NULL;
140
141 /* ====== compressed sparse columns === */
142 if (type & MATRIX_FORMAT_CSC) {
143 if (unroll) {
144 if (patternIsUnrolled) {
145 out->pattern=Paso_Pattern_getReference(pattern);
146 } else {
147 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
148 }
149 out->row_block_size=1;
150 out->col_block_size=1;
151 } else {
152 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
153 out->row_block_size=row_block_size;
154 out->col_block_size=col_block_size;
155 }
156 if (Paso_noError()) {
157 out->numRows = out->pattern->numInput;
158 out->numCols = out->pattern->numOutput;
159 }
160 } else {
161 /* ====== compressed sparse row === */
162 if (unroll) {
163 if (patternIsUnrolled) {
164 out->pattern=Paso_Pattern_getReference(pattern);
165 } else {
166 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
167 }
168 out->row_block_size=1;
169 out->col_block_size=1;
170 } else {
171 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
172 out->row_block_size=row_block_size;
173 out->col_block_size=col_block_size;
174 }
175 if (Paso_noError()) {
176 out->numRows = out->pattern->numOutput;
177 out->numCols = out->pattern->numInput;
178 }
179 }
180 if (Paso_noError()) {
181 out->block_size=out->row_block_size*out->col_block_size;
182 out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
183
184 out->val=MEMALLOC(out->len,double);
185 if (! Paso_checkPtr(out->val)) Paso_SparseMatrix_setValues(out,DBLE(0));
186 }
187 }
188 /* all done: */
189 if (Paso_noError()) {
190 return out;
191 } else {
192 Paso_SparseMatrix_free(out);
193 return NULL;
194 }
195 }
196
197 /* returns a reference to Paso_SparseMatrix in */
198
199 Paso_SparseMatrix* Paso_SparseMatrix_getReference(Paso_SparseMatrix* in) {
200 if (in!=NULL) ++(in->reference_counter);
201 return in;
202 }
203
204 /* deallocates a SparseMatrix: */
205
206 void Paso_SparseMatrix_free(Paso_SparseMatrix* in) {
207 if (in!=NULL) {
208 in->reference_counter--;
209 if (in->reference_counter<=0) {
210 Paso_Pattern_free(in->pattern);
211 MEMFREE(in->val);
212 MEMFREE(in);
213 #ifdef Paso_TRACE
214 printf("Paso_SparseMatrix_free: system matrix as been deallocated.\n");
215 #endif
216 }
217 }
218 }
219
220 Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p )
221 {
222 index_t *col_ind = NULL;
223 index_t *row_ind = NULL;
224 index_t *row_ptr = NULL;
225 double *val = NULL;
226 FILE *fileHandle_p = NULL;
227 Paso_Pattern* mainPattern=NULL;
228 Paso_SparseMatrix *out = NULL;
229 int i, curr_row, scan_ret;
230 MM_typecode matrixCode;
231 Paso_resetError();
232
233 /* open the file */
234 fileHandle_p = fopen( fileName_p, "r" );
235 if( fileHandle_p == NULL )
236 {
237 Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Cannot read file for reading.");
238 return NULL;
239 }
240
241 /* process banner */
242 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
243 {
244 Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner.");
245 fclose( fileHandle_p );
246 return NULL;
247 }
248 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
249 {
250
251 Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
252 fclose( fileHandle_p );
253 return NULL;
254 }
255
256 /* get matrix size */
257 if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
258 {
259 Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size");
260 fclose( fileHandle_p );
261 return NULL;
262 }
263
264 /* prepare storage */
265 col_ind = MEMALLOC( nz, index_t );
266 row_ind = MEMALLOC( nz, index_t );
267 val = MEMALLOC( nz, double );
268
269 row_ptr = MEMALLOC( (M+1), index_t );
270
271 if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
272 {
273 Paso_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory" );
274 fclose( fileHandle_p );
275 return NULL;
276 }
277
278 /* perform actual read of elements */
279 for( i=0; i<nz; i++ )
280 {
281 scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
282 if (scan_ret!=3)
283 {
284 MEMFREE( val );
285 MEMFREE( row_ind );
286 MEMFREE( col_ind );
287 MEMFREE( row_ptr );
288 fclose(fileHandle_p);
289 return NULL;
290 }
291 row_ind[i]--;
292 col_ind[i]--;
293 }
294 fclose( fileHandle_p );
295
296 /* sort the entries */
297 q_sort( row_ind, col_ind, val, 0, nz );
298
299 /* setup row_ptr */
300 curr_row = 0;
301 for( i=0; (i<nz && curr_row<M); curr_row++ )
302 {
303 while( row_ind[i] != curr_row )
304 i++;
305 row_ptr[curr_row] = i;
306 }
307 row_ptr[M] = nz;
308
309 mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,M,N,row_ptr,col_ind);
310 out = Paso_SparseMatrix_alloc(MATRIX_FORMAT_DEFAULT, mainPattern, 1, 1, TRUE);
311 /* copy values and cleanup temps */
312 for( i=0; i<nz; i++ ) out->val[i] = val[i];
313
314 Paso_Pattern_free(mainPattern);
315 MEMFREE( val );
316 MEMFREE( row_ind );
317 return out;
318 }
319
320
321 void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) {
322 FILE * fileHandle_p = NULL;
323 dim_t N,M,i, iptr_ij;
324 MM_typecode matcode;
325
326 if (A_p->col_block_size !=A_p->row_block_size) {
327 Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported.");
328 return;
329 }
330 if (A_p->row_block_size>3) {
331 Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM: currently only block size 3 is supported.\n");
332 return;
333 }
334
335 if (A_p->type & MATRIX_FORMAT_SYM) {
336 Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support symmetric storage scheme");
337 return;
338 }
339
340 /* open the file */
341 fileHandle_p = fopen(fileName_p, "w");
342 if (fileHandle_p==NULL) {
343 Paso_setError(IO_ERROR,"file could not be opened for writing");
344 return;
345 }
346
347 if (A_p->type & MATRIX_FORMAT_CSC) {
348 Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet.");
349 } else {
350 mm_initialize_typecode(&matcode);
351 mm_set_matrix(&matcode);
352 mm_set_coordinate(&matcode);
353 mm_set_real(&matcode);
354
355 N= A_p->numRows;
356 M= A_p->numCols;
357 mm_write_banner(fileHandle_p, matcode);
358 mm_write_mtx_crd_size(fileHandle_p, N, M, A_p->pattern->ptr[N]);
359 if(A_p->row_block_size==1)
360 for (i=0; i<N; i++)
361 for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
362 fprintf(fileHandle_p, "%d %d %25.15e\n", i+1, A_p->pattern->index[iptr_ij]+1, A_p->val[iptr_ij]);
363 }
364
365 if(A_p->row_block_size==2)
366 for (i=0; i<N; i++) {
367 for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
368 fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1, 2*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*4]);
369 fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1, 2*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*4+2]);
370 fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1+1, 2*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*4+1]);
371 fprintf(fileHandle_p, "%d %d %25.15e\n", i*2+1+1, 2*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*4+3]);
372 }
373 }
374
375 if(A_p->row_block_size==3)
376 for (i=0; i<N; i++) {
377 for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
378 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9]);
379 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+3]);
380 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+6]);
381 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9+1]);
382 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+4]);
383 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+1+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+7]);
384 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+1, A_p->val[iptr_ij*9+2]);
385 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+1+1, A_p->val[iptr_ij*9+5]);
386 fprintf(fileHandle_p, "%d %d %25.15e\n", i*3+2+1, 3*(A_p->pattern->index[iptr_ij])+2+1, A_p->val[iptr_ij*9+8]);
387 }
388 }
389
390 }
391
392 /* close the file */
393 fclose(fileHandle_p);
394
395 return;
396 }
397

  ViewVC Help
Powered by ViewVC 1.1.26