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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3642 - (show annotations)
Thu Oct 27 03:41:51 2011 UTC (7 years, 5 months ago) by caltinay
File MIME type: text/plain
File size: 12909 byte(s)
Assorted spelling/comment fixes in paso.

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 given type using the given matrix pattern.
104 Values are initialized with zero.
105 If patternIsUnrolled and type & MATRIX_FORMAT_BLK1, it is assumed that the
106 pattern is already unrolled to match the requested block size
107 and offsets. Otherwise unrolling and offset adjustment will be performed.
108 */
109 Paso_SparseMatrix* Paso_SparseMatrix_alloc(Paso_SparseMatrixType type,Paso_Pattern *pattern, int row_block_size, int col_block_size, const bool_t patternIsUnrolled) {
110
111 Paso_SparseMatrix*out=NULL;
112 Paso_SparseMatrixType pattern_format_out;
113 bool_t unroll=FALSE;
114
115 if (patternIsUnrolled) {
116 if (! XNOR(type & MATRIX_FORMAT_OFFSET1, pattern->type & MATRIX_FORMAT_OFFSET1) ) {
117 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_alloc: requested offset and pattern offset do 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 #ifndef USE_LAPACK
126 /* or any block size bigger than 3 */
127 || (col_block_size>3)
128 # endif
129 /* or if block size one requested and the block size is not 1 */
130 || ((type & MATRIX_FORMAT_BLK1) && (col_block_size>1) )
131 /* offsets don't match */
132 || ( (type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & MATRIX_FORMAT_OFFSET1) ) ;
133
134 pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? MATRIX_FORMAT_OFFSET1: MATRIX_FORMAT_DEFAULT;
135
136 Esys_resetError();
137 out=MEMALLOC(1,Paso_SparseMatrix);
138 if (! Esys_checkPtr(out)) {
139 out->pattern=NULL;
140 out->val=NULL;
141 out->reference_counter=1;
142 out->type=type;
143 out->solver_package=PASO_PASO;
144 out->solver_p=NULL;
145
146 /* ====== compressed sparse columns === */
147 if (type & MATRIX_FORMAT_CSC) {
148 if (unroll) {
149 if (patternIsUnrolled) {
150 out->pattern=Paso_Pattern_getReference(pattern);
151 } else {
152 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
153 }
154 out->row_block_size=1;
155 out->col_block_size=1;
156 } else {
157 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
158 out->row_block_size=row_block_size;
159 out->col_block_size=col_block_size;
160 }
161 if (Esys_noError()) {
162 out->numRows = out->pattern->numInput;
163 out->numCols = out->pattern->numOutput;
164 }
165
166 } else {
167 /* ====== compressed sparse row === */
168 if (unroll) {
169 if (patternIsUnrolled) {
170 out->pattern=Paso_Pattern_getReference(pattern);
171 } else {
172 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
173 }
174 out->row_block_size=1;
175 out->col_block_size=1;
176 } else {
177 out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
178 out->row_block_size=row_block_size;
179 out->col_block_size=col_block_size;
180 }
181 if (Esys_noError()) {
182 out->numRows = out->pattern->numOutput;
183 out->numCols = out->pattern->numInput;
184 }
185 }
186 if (Esys_noError()) {
187 if (type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
188 out->block_size=MIN(out->row_block_size,out->col_block_size);
189 } else {
190 out->block_size=out->row_block_size*out->col_block_size;
191 }
192 out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
193
194 out->val=MEMALLOC(out->len,double);
195 if (! Esys_checkPtr(out->val)) Paso_SparseMatrix_setValues(out,DBLE(0));
196 }
197 }
198 /* all done: */
199 if (Esys_noError()) {
200 return out;
201 } else {
202 Paso_SparseMatrix_free(out);
203 return NULL;
204 }
205 }
206
207 /* returns a reference to Paso_SparseMatrix in */
208
209 Paso_SparseMatrix* Paso_SparseMatrix_getReference(Paso_SparseMatrix* in) {
210 if (in!=NULL) ++(in->reference_counter);
211 return in;
212 }
213
214 /* deallocates a SparseMatrix: */
215
216 void Paso_SparseMatrix_free(Paso_SparseMatrix* in) {
217 if (in!=NULL) {
218 in->reference_counter--;
219 if (in->reference_counter<=0) {
220 switch(in->solver_package) {
221
222 case PASO_SMOOTHER:
223 Paso_Preconditioner_LocalSmoother_free((Paso_Preconditioner_LocalSmoother*) in->solver_p);
224 break;
225
226 case PASO_MKL:
227 Paso_MKL_free(in); /* releases solver_p */
228 break;
229
230 case PASO_UMFPACK:
231 Paso_UMFPACK_free(in); /* releases solver_p */
232 break;
233 }
234 MEMFREE(in->val);
235 Paso_Pattern_free(in->pattern);
236 MEMFREE(in);
237 }
238 }
239 }
240
241 Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p )
242 {
243 index_t *col_ind = NULL;
244 index_t *row_ind = NULL;
245 index_t *row_ptr = NULL;
246 double *val = NULL;
247 FILE *fileHandle_p = NULL;
248 Paso_Pattern* mainPattern=NULL;
249 Paso_SparseMatrix *out = NULL;
250 int i, curr_row, scan_ret;
251 MM_typecode matrixCode;
252 Esys_resetError();
253
254 /* open the file */
255 fileHandle_p = fopen( fileName_p, "r" );
256 if( fileHandle_p == NULL )
257 {
258 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Cannot open file for reading.");
259 return NULL;
260 }
261
262 /* process banner */
263 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
264 {
265 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner.");
266 fclose( fileHandle_p );
267 return NULL;
268 }
269 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
270 {
271
272 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
273 fclose( fileHandle_p );
274 return NULL;
275 }
276
277 /* get matrix size */
278 if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
279 {
280 Esys_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size.");
281 fclose( fileHandle_p );
282 return NULL;
283 }
284
285 /* prepare storage */
286 col_ind = MEMALLOC( nz, index_t );
287 row_ind = MEMALLOC( nz, index_t );
288 val = MEMALLOC( nz, double );
289
290 row_ptr = MEMALLOC( (M+1), index_t );
291
292 if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
293 {
294 Esys_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory." );
295 fclose( fileHandle_p );
296 return NULL;
297 }
298
299 /* perform actual read of elements */
300 for( i=0; i<nz; i++ )
301 {
302 scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
303 if (scan_ret!=3)
304 {
305 MEMFREE( val );
306 MEMFREE( row_ind );
307 MEMFREE( col_ind );
308 MEMFREE( row_ptr );
309 fclose(fileHandle_p);
310 return NULL;
311 }
312 row_ind[i]--;
313 col_ind[i]--;
314 }
315 fclose( fileHandle_p );
316
317 /* sort the entries */
318 q_sort( row_ind, col_ind, val, 0, nz );
319
320 /* setup row_ptr */
321 curr_row = 0;
322 for( i=0; (i<nz && curr_row<M); curr_row++ )
323 {
324 while( row_ind[i] != curr_row )
325 i++;
326 row_ptr[curr_row] = i;
327 }
328 row_ptr[M] = nz;
329
330 mainPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,M,N,row_ptr,col_ind);
331 out = Paso_SparseMatrix_alloc(MATRIX_FORMAT_DEFAULT, mainPattern, 1, 1, TRUE);
332 /* copy values and cleanup temps */
333 for( i=0; i<nz; i++ ) out->val[i] = val[i];
334
335 Paso_Pattern_free(mainPattern);
336 MEMFREE( val );
337 MEMFREE( row_ind );
338 return out;
339 }
340
341
342 void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) {
343 FILE * fileHandle_p = NULL;
344 dim_t N,M,i,j, irow, icol, ib, irb, icb, iptr;
345 MM_typecode matcode;
346 const dim_t col_block_size = A_p->col_block_size;
347 const dim_t row_block_size = A_p->row_block_size;
348 const dim_t block_size = A_p->block_size;
349 const index_t offset=(A_p->type & MATRIX_FORMAT_OFFSET1 ? 1:0);
350
351 if (col_block_size !=row_block_size) {
352 Esys_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported.");
353 return;
354 }
355
356 /* open the file */
357 fileHandle_p = fopen(fileName_p, "w");
358 if (fileHandle_p==NULL) {
359 Esys_setError(IO_ERROR,"Paso_SparseMatrix_saveMM: File could not be opened for writing");
360 return;
361 }
362 if (A_p->type & MATRIX_FORMAT_CSC) {
363 Esys_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet.");
364 } else {
365 mm_initialize_typecode(&matcode);
366 mm_set_matrix(&matcode);
367 mm_set_coordinate(&matcode);
368 mm_set_real(&matcode);
369
370 N=Paso_SparseMatrix_getNumRows(A_p);
371 M=Paso_SparseMatrix_getNumCols(A_p);
372 mm_write_banner(fileHandle_p, matcode);
373 mm_write_mtx_crd_size(fileHandle_p, N*row_block_size, M*col_block_size, A_p->pattern->ptr[N]*block_size);
374
375 if (A_p->type & MATRIX_FORMAT_DIAGONAL_BLOCK) {
376 for (i=0; i<N; i++) {
377 for (iptr = A_p->pattern->ptr[i]-offset;iptr<A_p->pattern->ptr[i+1]-offset; ++iptr) {
378 j=A_p->pattern->index[iptr]-offset;
379 for (ib=0;ib<block_size;ib++) {
380 irow=ib+row_block_size*i;
381 icol=ib+col_block_size*j;
382 fprintf(fileHandle_p, "%d %d %25.15e\n", irow+1, icol+1, A_p->val[iptr*block_size+ib]);
383 }
384 }
385 }
386 } else {
387 for (i=0; i<N; i++) {
388 for (iptr = A_p->pattern->ptr[i]-offset;iptr<A_p->pattern->ptr[i+1]-offset; ++iptr) {
389 j=A_p->pattern->index[iptr]-offset;
390 for (irb=0;irb<row_block_size;irb++) {
391 irow=irb+row_block_size*i;
392 for (icb=0;icb<col_block_size;icb++) {
393 icol=icb+col_block_size*j;
394 fprintf(fileHandle_p, "%d %d %25.15e\n", irow+1, icol+1, A_p->val[iptr*block_size+irb+row_block_size*icb]);
395 }
396 }
397 }
398 }
399 }
400 }
401 /* close the file */
402 fclose(fileHandle_p);
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 dim_t Paso_SparseMatrix_getTotalNumRows(const Paso_SparseMatrix* A){
424 return A->numRows * A->row_block_size;
425 }
426
427 dim_t Paso_SparseMatrix_getTotalNumCols(const Paso_SparseMatrix* A){
428 return A->numCols * A->col_block_size;
429 }
430 dim_t Paso_SparseMatrix_getNumRows(const Paso_SparseMatrix* A) {
431 return A->numRows;
432 }
433 dim_t Paso_SparseMatrix_getNumCols(const Paso_SparseMatrix* A) {
434 return A->numCols;
435 }
436
437 double Paso_SparseMatrix_getSize(const Paso_SparseMatrix* A)
438 {
439 if (A!=NULL) {
440 return DBLE(A->len);
441 } else {
442 return DBLE(0.);
443 }
444 }
445
446 double Paso_SparseMatrix_getSparsity(const Paso_SparseMatrix* A){
447 return Paso_SparseMatrix_getSize(A)/(DBLE(Paso_SparseMatrix_getTotalNumRows(A))*DBLE(Paso_SparseMatrix_getTotalNumCols(A)));
448 }

  ViewVC Help
Powered by ViewVC 1.1.26