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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3283 - (show annotations)
Mon Oct 18 22:39:28 2010 UTC (10 years, 4 months ago) by gross
File MIME type: text/plain
File size: 15318 byte(s)
AMG reengineered, BUG is Smoother fixed.


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: SparseMatrix */
18
19 /**************************************************************/
20
21 /* Author: Lutz Gross, l.gross@uq.edu.au */
22
23 /**************************************************************/
24
25 #include "Paso.h"
26 #include "SparseMatrix.h"
27 #include "MKL.h"
28 #include "Preconditioner.h"
29 #include "UMFPACK.h"
30 #include "TRILINOS.h"
31 #include "mmio.h"
32
33 /**************************************************************/
34 static void swap( index_t*, index_t*, double*, int, int );
35 static void q_sort( index_t*, index_t*, double*, int, int );
36 /*static void print_entries( index_t*, index_t*, double* );*/
37
38 static int M, N, nz;
39
40
41 /* debug: print the entries */
42 /*
43 void print_entries( index_t *r, index_t *c, double *v )
44 {
45 int i;
46
47 for( i=0; i<nz; i++ )
48 {
49 printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] );
50 }
51 }
52 */
53
54 /* swap function */
55 void swap( index_t *r, index_t *c, double *v, int left, int right )
56 {
57 double v_temp;
58 index_t temp;
59
60 temp = r[left];
61 r[left] = r[right];
62 r[right] = temp;
63
64 temp = c[left];
65 c[left] = c[right];
66 c[right] = temp;
67
68 v_temp = v[left];
69 v[left] = v[right];
70 v[right] = v_temp;
71 }
72
73 void q_sort( index_t *row, index_t *col, double *val, int begin, int end )
74 {
75 int l, r;
76 index_t pivot, lval;
77
78 if( end > begin )
79 {
80 pivot = N * row[begin] + col[begin];
81 l = begin + 1;
82 r = end;
83
84 while( l < r )
85 {
86 lval = N * row[l] + col[l];
87 if( lval < pivot )
88 l++;
89 else
90 {
91 r--;
92 swap( row, col, val, l, r );
93 }
94 }
95 l--;
96 swap( row, col, val, begin, l );
97 q_sort( row, col, val, begin, l );
98 q_sort( row, col, val, r, end );
99 }
100 }
101
102
103 /* allocates a SparseMatrix of type type using the given matrix pattern
104 Values are initialized by zero.
105 if patternIsUnrolled and type & MATRIX_FORMAT_BLK1, it is assumed that the pattern is allready unrolled to match the requested block size
106 and offsets otherwise unrolling and offset adjustment will be performed.
107 */
108 Paso_SparseMatrix* Paso_SparseMatrix_alloc(Paso_SparseMatrixType type,Paso_Pattern *pattern, int row_block_size, int col_block_size, const bool_t patternIsUnrolled) {
109
110 Paso_SparseMatrix*out=NULL;
111 Paso_SparseMatrixType pattern_format_out;
112 bool_t unroll=FALSE;
113
114 if (patternIsUnrolled) {
115 if (! XNOR(type & MATRIX_FORMAT_OFFSET1, pattern->type & PATTERN_FORMAT_OFFSET1) ) {
116 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_alloc: requested offset and pattern offset does not match.");
117 return NULL;
118 }
119 }
120 /* do we need to apply unrolling ? */
121 unroll
122 /* we don't like non-square blocks */
123 = (row_block_size!=col_block_size)
124 /* or any block size bigger than 3 */
125 || (col_block_size>3)
126 /* or if lock size one requested and the block size is not 1 */
127 || ((type & MATRIX_FORMAT_BLK1) && (col_block_size>1) )
128 /* offsets don't match */
129 || ( (type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & PATTERN_FORMAT_OFFSET1) ) ;
130
131 pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT;
132
133 Esys_resetError();
134 out=MEMALLOC(1,Paso_SparseMatrix);
135 if (! Esys_checkPtr(out)) {
136 out->pattern=NULL;
137 out->val=NULL;
138 out->reference_counter=1;
139 out->type=type;
140 out->solver_package=PASO_PASO;
141 out->solver_p=NULL;
142
143 /* ====== compressed sparse columns === */
144 if (type & MATRIX_FORMAT_CSC) {
145 if (unroll) {
146 if (patternIsUnrolled) {
147 out->pattern=Paso_Pattern_getReference(pattern);
148 } else {
149 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
150 }
151 out->row_block_size=1;
152 out->col_block_size=1;
153 } else {
154 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
155 out->row_block_size=row_block_size;
156 out->col_block_size=col_block_size;
157 }
158 if (Esys_noError()) {
159 out->numRows = out->pattern->numInput;
160 out->numCols = out->pattern->numOutput;
161 }
162 } else {
163 /* ====== compressed sparse row === */
164 if (unroll) {
165 if (patternIsUnrolled) {
166 out->pattern=Paso_Pattern_getReference(pattern);
167 } else {
168 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
169 }
170 out->row_block_size=1;
171 out->col_block_size=1;
172 } else {
173 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
174 out->row_block_size=row_block_size;
175 out->col_block_size=col_block_size;
176 }
177 if (Esys_noError()) {
178 out->numRows = out->pattern->numOutput;
179 out->numCols = out->pattern->numInput;
180 }
181 }
182 if (Esys_noError()) {
183 out->block_size=out->row_block_size*out->col_block_size;
184 out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
185
186 out->val=MEMALLOC(out->len,double);
187 if (! Esys_checkPtr(out->val)) Paso_SparseMatrix_setValues(out,DBLE(0));
188 }
189 }
190 /* all done: */
191 if (Esys_noError()) {
192 return out;
193 } else {
194 Paso_SparseMatrix_free(out);
195 return NULL;
196 }
197 }
198
199 /* returns a reference to Paso_SparseMatrix in */
200
201 Paso_SparseMatrix* Paso_SparseMatrix_getReference(Paso_SparseMatrix* in) {
202 if (in!=NULL) ++(in->reference_counter);
203 return in;
204 }
205
206 /* deallocates a SparseMatrix: */
207
208 void Paso_SparseMatrix_free(Paso_SparseMatrix* in) {
209 if (in!=NULL) {
210 in->reference_counter--;
211 if (in->reference_counter<=0) {
212 switch(in->solver_package) {
213
214 case PASO_SMOOTHER:
215 Paso_Preconditioner_LocalSmoother_free((Paso_Preconditioner_LocalSmoother*) in->solver_p);
216 break;
217
218 case PASO_MKL:
219 Paso_MKL_free(in); /* releases solver_p */
220 break;
221
222 case PASO_UMFPACK:
223 Paso_UMFPACK_free(in); /* releases solver_p */
224 break;
225 }
226 MEMFREE(in->val);
227 Paso_Pattern_free(in->pattern);
228 MEMFREE(in);
229 }
230 }
231 }
232
233 Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p )
234 {
235 index_t *col_ind = NULL;
236 index_t *row_ind = NULL;
237 index_t *row_ptr = NULL;
238 double *val = NULL;
239 FILE *fileHandle_p = NULL;
240 Paso_Pattern* mainPattern=NULL;
241 Paso_SparseMatrix *out = NULL;
242 int i, curr_row, scan_ret;
243 MM_typecode matrixCode;
244 Esys_resetError();
245
246 /* open the file */
247 fileHandle_p = fopen( fileName_p, "r" );
248 if( fileHandle_p == NULL )
249 {
250 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Cannot read file for reading.");
251 return NULL;
252 }
253
254 /* process banner */
255 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
256 {
257 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner.");
258 fclose( fileHandle_p );
259 return NULL;
260 }
261 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
262 {
263
264 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
265 fclose( fileHandle_p );
266 return NULL;
267 }
268
269 /* get matrix size */
270 if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
271 {
272 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size");
273 fclose( fileHandle_p );
274 return NULL;
275 }
276
277 /* prepare storage */
278 col_ind = MEMALLOC( nz, index_t );
279 row_ind = MEMALLOC( nz, index_t );
280 val = MEMALLOC( nz, double );
281
282 row_ptr = MEMALLOC( (M+1), index_t );
283
284 if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
285 {
286 Esys_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory" );
287 fclose( fileHandle_p );
288 return NULL;
289 }
290
291 /* perform actual read of elements */
292 for( i=0; i<nz; i++ )
293 {
294 scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
295 if (scan_ret!=3)
296 {
297 MEMFREE( val );
298 MEMFREE( row_ind );
299 MEMFREE( col_ind );
300 MEMFREE( row_ptr );
301 fclose(fileHandle_p);
302 return NULL;
303 }
304 row_ind[i]--;
305 col_ind[i]--;
306 }
307 fclose( fileHandle_p );
308
309 /* sort the entries */
310 q_sort( row_ind, col_ind, val, 0, nz );
311
312 /* setup row_ptr */
313 curr_row = 0;
314 for( i=0; (i<nz && curr_row<M); curr_row++ )
315 {
316 while( row_ind[i] != curr_row )
317 i++;
318 row_ptr[curr_row] = i;
319 }
320 row_ptr[M] = nz;
321
322 mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,M,N,row_ptr,col_ind);
323 out = Paso_SparseMatrix_alloc(MATRIX_FORMAT_DEFAULT, mainPattern, 1, 1, TRUE);
324 /* copy values and cleanup temps */
325 for( i=0; i<nz; i++ ) out->val[i] = val[i];
326
327 Paso_Pattern_free(mainPattern);
328 MEMFREE( val );
329 MEMFREE( row_ind );
330 return out;
331 }
332
333
334 void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) {
335 FILE * fileHandle_p = NULL;
336 dim_t N,M,i, iptr_ij;
337 MM_typecode matcode;
338
339 if (A_p->col_block_size !=A_p->row_block_size) {
340 Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported.");
341 return;
342 }
343 if (A_p->row_block_size>3) {
344 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM: currently only block size 3 is supported.\n");
345 return;
346 }
347
348 /* open the file */
349 fileHandle_p = fopen(fileName_p, "w");
350 if (fileHandle_p==NULL) {
351 Esys_setError(IO_ERROR,"file could not be opened for writing");
352 return;
353 }
354
355 if (A_p->type & MATRIX_FORMAT_CSC) {
356 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet.");
357 } else {
358 mm_initialize_typecode(&matcode);
359 mm_set_matrix(&matcode);
360 mm_set_coordinate(&matcode);
361 mm_set_real(&matcode);
362
363 N= A_p->numRows;
364 M= A_p->numCols;
365 mm_write_banner(fileHandle_p, matcode);
366 mm_write_mtx_crd_size(fileHandle_p, N, M, A_p->pattern->ptr[N]);
367 if(A_p->row_block_size==1)
368 for (i=0; i<N; i++)
369 for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
370 fprintf(fileHandle_p, "%d %d %25.15e\n", i+1, A_p->pattern->index[iptr_ij]+1, A_p->val[iptr_ij]);
371 }
372
373 if(A_p->row_block_size==2)
374 for (i=0; i<N; i++) {
375 for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
376 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]);
377 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]);
378 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]);
379 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]);
380 }
381 }
382
383 if(A_p->row_block_size==3)
384 for (i=0; i<N; i++) {
385 for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
386 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]);
387 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]);
388 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]);
389 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]);
390 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]);
391 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]);
392 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]);
393 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]);
394 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]);
395 }
396 }
397
398 }
399
400 /* close the file */
401 fclose(fileHandle_p);
402
403 return;
404 }
405
406 index_t* Paso_SparseMatrix_borrowMainDiagonalPointer(Paso_SparseMatrix * A_p)
407 {
408 return Paso_Pattern_borrowMainDiagonalPointer(A_p->pattern);
409 }
410
411 dim_t Paso_SparseMatrix_getNumColors(Paso_SparseMatrix* A_p)
412 {
413 return Paso_Pattern_getNumColors(A_p->pattern);
414 }
415 index_t* Paso_SparseMatrix_borrowColoringPointer(Paso_SparseMatrix* A_p)
416 {
417 return Paso_Pattern_borrowColoringPointer(A_p->pattern);
418 }
419 dim_t Paso_SparseMatrix_maxDeg(Paso_SparseMatrix * A_p)
420 {
421 return Paso_Pattern_maxDeg(A_p->pattern);
422 }
423
424 Paso_SparseMatrix* Paso_SparseMatrix_unroll(Paso_SparseMatrix* A) {
425 Paso_SparseMatrix *out = NULL;
426 /*Paso_Pattern* mainPattern=NULL;*/
427 int i;
428 dim_t block_size=A->row_block_size;
429 index_t iptr, optr1=0, optr2=0, optr3=0;
430
431 /*mainPattern= Paso_Pattern_unrollBlocks(A->pattern, PATTERN_FORMAT_DEFAULT, A->row_block_size,A->col_block_size);*/
432 out = Paso_SparseMatrix_alloc(MATRIX_FORMAT_BLK1, A->pattern, A->row_block_size,A->col_block_size, FALSE);
433
434 if (block_size==1) {
435 #pragma omp parallel for private(i,iptr) schedule(static)
436 for (i=0;i<A->numRows;++i) {
437 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
438 out->val[iptr]=A->val[iptr];
439 }
440 }
441 } else if (block_size==2) {
442 #pragma omp parallel for private(i,iptr,optr1,optr2) schedule(static)
443 for (i=0;i<A->numRows;++i) {
444 optr1=out->pattern->ptr[i*block_size];
445 optr2=out->pattern->ptr[i*block_size+1];
446 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
447 out->val[optr1]=A->val[iptr*block_size*block_size];
448 optr1++;
449 out->val[optr1]=A->val[iptr*block_size*block_size+2];
450 optr1++;
451 out->val[optr2]=A->val[iptr*block_size*block_size+1];
452 optr2++;
453 out->val[optr2]=A->val[iptr*block_size*block_size+3];
454 optr2++;
455 }
456 }
457
458 } else if (block_size==3) {
459 #pragma omp parallel for private(i,iptr,optr1,optr2,optr3) schedule(static)
460 for (i=0;i<A->numRows;++i) {
461 optr1=out->pattern->ptr[i*block_size];
462 optr2=out->pattern->ptr[i*block_size+1];
463 optr3=out->pattern->ptr[i*block_size+2];
464 for (iptr=A->pattern->ptr[i];iptr<A->pattern->ptr[i+1]; ++iptr) {
465 out->val[optr1]=A->val[iptr*block_size*block_size];
466 optr1++;
467 out->val[optr1]=A->val[iptr*block_size*block_size+3];
468 optr1++;
469 out->val[optr1]=A->val[iptr*block_size*block_size+6];
470 optr1++;
471 out->val[optr2]=A->val[iptr*block_size*block_size+1];
472 optr2++;
473 out->val[optr2]=A->val[iptr*block_size*block_size+4];
474 optr2++;
475 out->val[optr2]=A->val[iptr*block_size*block_size+7];
476 optr2++;
477 out->val[optr3]=A->val[iptr*block_size*block_size+2];
478 optr3++;
479 out->val[optr3]=A->val[iptr*block_size*block_size+5];
480 optr3++;
481 out->val[optr3]=A->val[iptr*block_size*block_size+8];
482 optr3++;
483 }
484 }
485
486
487 }
488
489 return out;
490 }

  ViewVC Help
Powered by ViewVC 1.1.26