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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2554 - (hide annotations)
Fri Jul 24 05:38:54 2009 UTC (11 years, 2 months ago) by gross
File MIME type: text/plain
File size: 12493 byte(s)
bug with MKL call fixed.
1 ksteube 1313
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 ksteube 1313
14 ksteube 1811
15 ksteube 1313 /**************************************************************/
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 artak 2240 #include "mmio.h"
29 ksteube 1313
30     /**************************************************************/
31 artak 2240 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 ksteube 1313
35 artak 2240 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 ksteube 1313 /* allocates a SparseMatrix of type type using the given matrix pattern
101 gross 2551 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 ksteube 1313
107     Paso_SparseMatrix*out=NULL;
108     Paso_SparseMatrixType pattern_format_out;
109 gross 2551 bool_t unroll=FALSE;
110 ksteube 1313
111 gross 2551 if (type & MATRIX_FORMAT_SYM) {
112     Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_alloc: symmetric matrix pattern are not supported.");
113     return NULL;
114     }
115 gross 2554
116 gross 2551 if (patternIsUnrolled) {
117 gross 2554 if (! XNOR(type & MATRIX_FORMAT_OFFSET1, pattern->type & PATTERN_FORMAT_OFFSET1) ) {
118     Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_alloc: requested offset and pattern offset does not match.");
119 gross 2551 return NULL;
120     }
121     }
122     /* do we need to apply unrolling ? */
123     unroll
124     /* we don't like non-square blocks */
125     = (row_block_size!=col_block_size)
126     /* or any block size bigger than 3 */
127     || (col_block_size>3)
128     /* or if lock size one requested and the block size is not 1 */
129 gross 2554 || ((type & MATRIX_FORMAT_BLK1) && (col_block_size>1) )
130     /* offsets don't match */
131     || ( (type & MATRIX_FORMAT_OFFSET1) != ( pattern->type & PATTERN_FORMAT_OFFSET1) ) ;
132 gross 2551
133     pattern_format_out= (type & MATRIX_FORMAT_OFFSET1)? PATTERN_FORMAT_OFFSET1: PATTERN_FORMAT_DEFAULT;
134    
135 ksteube 1313 Paso_resetError();
136     out=MEMALLOC(1,Paso_SparseMatrix);
137     if (! Paso_checkPtr(out)) {
138     out->pattern=NULL;
139     out->val=NULL;
140     out->reference_counter=1;
141     out->type=type;
142 artak 2489 out->solver=NULL;
143 ksteube 1313
144     /* ====== compressed sparse columns === */
145     if (type & MATRIX_FORMAT_CSC) {
146 gross 2551 if (unroll) {
147     if (patternIsUnrolled) {
148     out->pattern=Paso_Pattern_getReference(pattern);
149     } else {
150     out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,col_block_size,row_block_size);
151     }
152 ksteube 1313 out->row_block_size=1;
153     out->col_block_size=1;
154     } else {
155     out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
156     out->row_block_size=row_block_size;
157     out->col_block_size=col_block_size;
158     }
159 gross 2551 if (Paso_noError()) {
160     out->numRows = out->pattern->numInput;
161     out->numCols = out->pattern->numOutput;
162     }
163 ksteube 1313 } else {
164     /* ====== compressed sparse row === */
165 gross 2551 if (unroll) {
166     if (patternIsUnrolled) {
167     out->pattern=Paso_Pattern_getReference(pattern);
168     } else {
169     out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,row_block_size,col_block_size);
170     }
171 ksteube 1313 out->row_block_size=1;
172     out->col_block_size=1;
173     } else {
174     out->pattern=Paso_Pattern_unrollBlocks(pattern,pattern_format_out,1,1);
175     out->row_block_size=row_block_size;
176     out->col_block_size=col_block_size;
177     }
178 gross 2551 if (Paso_noError()) {
179     out->numRows = out->pattern->numOutput;
180     out->numCols = out->pattern->numInput;
181     }
182 ksteube 1313 }
183 gross 2551 if (Paso_noError()) {
184     out->block_size=out->row_block_size*out->col_block_size;
185     out->len=(size_t)(out->pattern->len)*(size_t)(out->block_size);
186    
187     out->val=MEMALLOC(out->len,double);
188     if (! Paso_checkPtr(out->val)) Paso_SparseMatrix_setValues(out,DBLE(0));
189     }
190 ksteube 1313 }
191     /* all done: */
192 gross 2551 if (Paso_noError()) {
193     return out;
194     } else {
195 ksteube 1313 Paso_SparseMatrix_free(out);
196     return NULL;
197     }
198     }
199    
200     /* returns a reference to Paso_SparseMatrix in */
201    
202     Paso_SparseMatrix* Paso_SparseMatrix_getReference(Paso_SparseMatrix* in) {
203     if (in!=NULL) ++(in->reference_counter);
204     return in;
205     }
206    
207     /* deallocates a SparseMatrix: */
208    
209     void Paso_SparseMatrix_free(Paso_SparseMatrix* in) {
210     if (in!=NULL) {
211     in->reference_counter--;
212     if (in->reference_counter<=0) {
213     Paso_Pattern_free(in->pattern);
214     MEMFREE(in->val);
215     MEMFREE(in);
216     #ifdef Paso_TRACE
217     printf("Paso_SparseMatrix_free: system matrix as been deallocated.\n");
218     #endif
219     }
220     }
221     }
222 artak 2240
223     Paso_SparseMatrix* Paso_SparseMatrix_loadMM_toCSR( char *fileName_p )
224     {
225     index_t *col_ind = NULL;
226     index_t *row_ind = NULL;
227     index_t *row_ptr = NULL;
228     double *val = NULL;
229     FILE *fileHandle_p = NULL;
230     Paso_Pattern* mainPattern=NULL;
231     Paso_SparseMatrix *out = NULL;
232     int i, curr_row, scan_ret;
233     MM_typecode matrixCode;
234     Paso_resetError();
235    
236     /* open the file */
237     fileHandle_p = fopen( fileName_p, "r" );
238     if( fileHandle_p == NULL )
239     {
240     Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Cannot read file for reading.");
241     return NULL;
242     }
243    
244     /* process banner */
245     if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
246     {
247     Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Error processing MM banner.");
248     fclose( fileHandle_p );
249     return NULL;
250     }
251     if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
252     {
253    
254     Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
255     fclose( fileHandle_p );
256     return NULL;
257     }
258    
259     /* get matrix size */
260     if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
261     {
262     Paso_setError(IO_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not parse matrix size");
263     fclose( fileHandle_p );
264     return NULL;
265     }
266    
267     /* prepare storage */
268     col_ind = MEMALLOC( nz, index_t );
269     row_ind = MEMALLOC( nz, index_t );
270     val = MEMALLOC( nz, double );
271    
272     row_ptr = MEMALLOC( (M+1), index_t );
273    
274     if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
275     {
276     Paso_setError(MEMORY_ERROR, "Paso_SparseMatrix_loadMM_toCSR: Could not allocate memory" );
277     fclose( fileHandle_p );
278     return NULL;
279     }
280    
281     /* perform actual read of elements */
282     for( i=0; i<nz; i++ )
283     {
284     scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
285     if (scan_ret!=3)
286     {
287     MEMFREE( val );
288     MEMFREE( row_ind );
289     MEMFREE( col_ind );
290     MEMFREE( row_ptr );
291     fclose(fileHandle_p);
292     return NULL;
293     }
294     row_ind[i]--;
295     col_ind[i]--;
296     }
297     fclose( fileHandle_p );
298    
299     /* sort the entries */
300     q_sort( row_ind, col_ind, val, 0, nz );
301    
302     /* setup row_ptr */
303     curr_row = 0;
304     for( i=0; (i<nz && curr_row<M); curr_row++ )
305     {
306     while( row_ind[i] != curr_row )
307     i++;
308     row_ptr[curr_row] = i;
309     }
310     row_ptr[M] = nz;
311    
312 gross 2551 mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,M,N,row_ptr,col_ind);
313     out = Paso_SparseMatrix_alloc(MATRIX_FORMAT_DEFAULT, mainPattern, 1, 1, TRUE);
314 artak 2240 /* copy values and cleanup temps */
315     for( i=0; i<nz; i++ ) out->val[i] = val[i];
316    
317     Paso_Pattern_free(mainPattern);
318     MEMFREE( val );
319     MEMFREE( row_ind );
320     return out;
321     }
322    
323    
324     void Paso_SparseMatrix_saveMM(Paso_SparseMatrix * A_p, char * fileName_p) {
325     FILE * fileHandle_p = NULL;
326     dim_t N,M,i, iptr_ij;
327     MM_typecode matcode;
328    
329     if (A_p->col_block_size !=A_p->row_block_size) {
330     Paso_setError(TYPE_ERROR, "Paso_SparseMatrix_saveMM: currently only square blocks are supported.");
331     return;
332     }
333     if (A_p->row_block_size>3) {
334     Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM: currently only block size 3 is supported.\n");
335     return;
336     }
337    
338     if (A_p->type & MATRIX_FORMAT_SYM) {
339     Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support symmetric storage scheme");
340     return;
341     }
342    
343     /* open the file */
344     fileHandle_p = fopen(fileName_p, "w");
345     if (fileHandle_p==NULL) {
346     Paso_setError(IO_ERROR,"file could not be opened for writing");
347     return;
348     }
349    
350     if (A_p->type & MATRIX_FORMAT_CSC) {
351     Paso_setError(TYPE_ERROR,"Paso_SparseMatrix_saveMM does not support CSC yet.");
352     } else {
353     mm_initialize_typecode(&matcode);
354     mm_set_matrix(&matcode);
355     mm_set_coordinate(&matcode);
356     mm_set_real(&matcode);
357    
358     N= A_p->numRows;
359     M= A_p->numCols;
360     mm_write_banner(fileHandle_p, matcode);
361     mm_write_mtx_crd_size(fileHandle_p, N, M, A_p->pattern->ptr[N]);
362     if(A_p->row_block_size==1)
363     for (i=0; i<N; i++)
364     for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
365     fprintf(fileHandle_p, "%d %d %25.15e\n", i+1, A_p->pattern->index[iptr_ij]+1, A_p->val[iptr_ij]);
366     }
367    
368     if(A_p->row_block_size==2)
369     for (i=0; i<N; i++) {
370     for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
371     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]);
372     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]);
373     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]);
374     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]);
375     }
376     }
377    
378     if(A_p->row_block_size==3)
379     for (i=0; i<N; i++) {
380     for (iptr_ij=A_p->pattern->ptr[i];iptr_ij<A_p->pattern->ptr[i+1]; ++iptr_ij) {
381     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]);
382     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]);
383     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]);
384     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]);
385     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]);
386     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]);
387     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]);
388     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]);
389     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]);
390     }
391     }
392    
393     }
394    
395     /* close the file */
396     fclose(fileHandle_p);
397    
398     return;
399     }
400    

  ViewVC Help
Powered by ViewVC 1.1.26