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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26