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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26