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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1562 - (hide annotations)
Wed May 21 13:04:40 2008 UTC (11 years, 4 months ago) by gross
File MIME type: text/plain
File size: 9771 byte(s)
The algebraic upwinding with MPI. The case of boundary constraint needs still some attention. 


1 ksteube 1312
2 jgs 150 /* $Id$ */
3    
4 ksteube 1312 /*******************************************************
5     *
6     * Copyright 2003-2007 by ACceSS MNRF
7     * Copyright 2007 by University of Queensland
8     *
9     * http://esscc.uq.edu.au
10     * Primary Business: Queensland, Australia
11     * Licensed under the Open Software License version 3.0
12     * http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15 dhawcroft 631
16 jgs 150 /**************************************************************/
17    
18     /* Paso: Matrix Market format is loaded to a SystemMatrix */
19    
20     /**************************************************************/
21    
22     /* Copyrights by ACcESS Australia 2003,2004,2005 */
23     /* Author: imran@access.edu.au */
24    
25     /**************************************************************/
26    
27     #include "Paso.h"
28     #include "mmio.h"
29     #include "SystemMatrix.h"
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     void print_entries( index_t *r, index_t *c, double *v )
40     {
41     int i;
42    
43     for( i=0; i<nz; i++ )
44     {
45     printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] );
46     }
47     }
48    
49     /* swap function */
50     void swap( index_t *r, index_t *c, double *v, int left, int right )
51     {
52     double v_temp;
53     index_t temp;
54    
55     temp = r[left];
56     r[left] = r[right];
57     r[right] = temp;
58    
59     temp = c[left];
60     c[left] = c[right];
61     c[right] = temp;
62    
63     v_temp = v[left];
64     v[left] = v[right];
65     v[right] = v_temp;
66     }
67    
68     void q_sort( index_t *row, index_t *col, double *val, int begin, int end )
69     {
70     int l, r;
71     index_t pivot, lval;
72    
73     if( end > begin )
74     {
75     pivot = N * row[begin] + col[begin];
76     l = begin + 1;
77     r = end;
78    
79     while( l < r )
80     {
81     lval = N * row[l] + col[l];
82     if( lval < pivot )
83     l++;
84     else
85     {
86     r--;
87     swap( row, col, val, l, r );
88     }
89     }
90     l--;
91     swap( row, col, val, begin, l );
92     q_sort( row, col, val, begin, l );
93     q_sort( row, col, val, r, end );
94     }
95     }
96    
97     Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSR( char *fileName_p )
98     {
99 ksteube 1312 index_t dist[2];
100     Paso_Distribution* input_dist=NULL, *output_dist=NULL;
101 jgs 150 index_t *col_ind = NULL;
102     index_t *row_ind = NULL;
103     index_t *row_ptr = NULL;
104     double *val = NULL;
105 gross 1028 FILE *fileHandle_p = NULL;
106 ksteube 1312 Paso_Pattern* mainPattern=NULL, *couplePattern=NULL;
107     Paso_SystemMatrixPattern *pattern = NULL;
108 jgs 150 Paso_SystemMatrix *out = NULL;
109 ksteube 1312 Paso_SharedComponents *send =NULL;
110 gross 1552 Paso_Connector *connector=NULL;
111 ksteube 1312 int i, curr_row;
112     MM_typecode matrixCode;
113     Paso_MPIInfo* mpi_info=Paso_MPIInfo_alloc( MPI_COMM_WORLD);
114 gross 1028 Paso_resetError();
115 ksteube 1312 if (mpi_info->size >1) {
116     Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: support single processor only");
117     return NULL;
118     }
119 jgs 150 /* open the file */
120 gross 1028 fileHandle_p = fopen( fileName_p, "r" );
121 jgs 150 if( fileHandle_p == NULL )
122     {
123 gross 415 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Cannot read file for reading.");
124 ksteube 1312 Paso_MPIInfo_free(mpi_info);
125 jgs 150 return NULL;
126     }
127    
128     /* process banner */
129     if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
130     {
131 gross 415 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Error processing MM banner.");
132 ksteube 1312 Paso_MPIInfo_free(mpi_info);
133 jgs 150 fclose( fileHandle_p );
134     return NULL;
135     }
136     if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
137     {
138 gross 415
139     Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
140 ksteube 1312 Paso_MPIInfo_free(mpi_info);
141 jgs 150 fclose( fileHandle_p );
142     return NULL;
143     }
144    
145     /* get matrix size */
146     if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
147     {
148 gross 415 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not parse matrix size");
149 ksteube 1312 Paso_MPIInfo_free(mpi_info);
150 jgs 150 fclose( fileHandle_p );
151     return NULL;
152     }
153    
154     /* prepare storage */
155     col_ind = MEMALLOC( nz, index_t );
156     row_ind = MEMALLOC( nz, index_t );
157     val = MEMALLOC( nz, double );
158    
159     row_ptr = MEMALLOC( (M+1), index_t );
160    
161     if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
162     {
163 gross 415 Paso_setError(MEMORY_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not allocate memory" );
164 jgs 150
165 ksteube 1312 Paso_MPIInfo_free(mpi_info);
166 jgs 150 fclose( fileHandle_p );
167     return NULL;
168     }
169    
170     /* perform actual read of elements */
171     for( i=0; i<nz; i++ )
172     {
173     fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
174     row_ind[i]--;
175     col_ind[i]--;
176     }
177     fclose( fileHandle_p );
178    
179     /* sort the entries */
180     q_sort( row_ind, col_ind, val, 0, nz );
181    
182     /* setup row_ptr */
183     curr_row = 0;
184     for( i=0; (i<nz && curr_row<M); curr_row++ )
185     {
186     while( row_ind[i] != curr_row )
187     i++;
188     row_ptr[curr_row] = i;
189     }
190     row_ptr[M] = nz;
191    
192 ksteube 1312 /* create return value */
193 jgs 150 /* create F_SMP and F_SM */
194 ksteube 1312 dist[0]=0;
195     dist[1]=M;
196     output_dist=Paso_Distribution_alloc(mpi_info, dist,1,0);
197     dist[1]=N;
198     input_dist=Paso_Distribution_alloc(mpi_info, dist,1,0);
199     mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,row_ptr,row_ind);
200     couplePattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,NULL,NULL);
201 gross 1562 send=Paso_SharedComponents_alloc(M,0,NULL,NULL,NULL,1,0,mpi_info);
202 gross 1552 connector=Paso_Connector_alloc(send,send);
203 ksteube 1312 pattern=Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,output_dist,input_dist,
204 gross 1552 mainPattern,couplePattern,couplePattern,connector,connector);
205 jgs 150
206 ksteube 1312 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_DEFAULT, pattern, 1, 1);
207 jgs 150 /* copy values and cleanup temps */
208 ksteube 1312 for( i=0; i<nz; i++ ) out->mainBlock->val[i] = val[i];
209 jgs 150
210 ksteube 1312 Paso_SystemMatrixPattern_free(pattern);
211     Paso_Pattern_free(mainPattern);
212     Paso_Pattern_free(couplePattern);
213 gross 1552 Paso_Connector_free(connector);
214 ksteube 1312 Paso_Distribution_free(output_dist);
215     Paso_Distribution_free(input_dist);
216     Paso_SharedComponents_free(send);
217     Paso_MPIInfo_free(mpi_info);
218 jgs 150 MEMFREE( val );
219     MEMFREE( row_ind );
220    
221     return out;
222     }
223    
224     Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSC( char *fileName_p )
225     {
226 ksteube 1312 index_t dist[2];
227     Paso_Distribution* input_dist=NULL, *output_dist=NULL;
228 gross 1028 FILE *fileHandle_p = NULL;
229 ksteube 1312 Paso_Pattern* mainPattern=NULL, *couplePattern=NULL;
230     Paso_SystemMatrixPattern *pattern = NULL;
231     Paso_SystemMatrix *out = NULL;
232     Paso_SharedComponents *send =NULL;
233 gross 1552 Paso_Connector *connector=NULL;
234 jgs 150 index_t *col_ind = NULL;
235     index_t *row_ind = NULL;
236     index_t *col_ptr = NULL;
237     double *val = NULL;
238 ksteube 1312 int i, curr_col=0;
239     MM_typecode matrixCode;
240     Paso_MPIInfo* mpi_info=Paso_MPIInfo_alloc( MPI_COMM_WORLD);
241     if (mpi_info->size >1) {
242     Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSC: support single processor only");
243     return NULL;
244     }
245 jgs 150
246     Paso_resetError();
247    
248     /* open the file */
249 gross 1028 fileHandle_p = fopen( fileName_p, "r" );
250 jgs 150 if( fileHandle_p == NULL )
251     {
252 gross 415 Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: File could not be opened for reading");
253 ksteube 1312 Paso_MPIInfo_free(mpi_info);
254 jgs 150 return NULL;
255     }
256    
257     /* process banner */
258     if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
259     {
260 gross 415 Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: Error processing MM banner");
261 jgs 150 fclose( fileHandle_p );
262 ksteube 1312 Paso_MPIInfo_free(mpi_info);
263 jgs 150 return NULL;
264     }
265     if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
266     {
267 gross 415 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
268 jgs 150 fclose( fileHandle_p );
269 ksteube 1312 Paso_MPIInfo_free(mpi_info);
270 jgs 150 return NULL;
271     }
272    
273     /* get matrix size */
274     if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
275     {
276 gross 415 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
277 jgs 150 fclose( fileHandle_p );
278 ksteube 1312 Paso_MPIInfo_free(mpi_info);
279 jgs 150 return NULL;
280     }
281    
282     /* prepare storage */
283     col_ind = MEMALLOC( nz, index_t );
284     row_ind = MEMALLOC( nz, index_t );
285     val = MEMALLOC( nz, double );
286    
287     col_ptr = MEMALLOC( (N+1), index_t );
288    
289    
290     /* perform actual read of elements */
291     for( i=0; i<nz; i++ )
292     {
293     fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
294     row_ind[i]--;
295     col_ind[i]--;
296     }
297     fclose( fileHandle_p );
298    
299     /* sort the entries */
300     q_sort( col_ind, row_ind, val, 0, nz );
301    
302     /* setup row_ptr */
303     for( i=0; (i<nz && curr_col<N); curr_col++ )
304     {
305     while( col_ind[i] != curr_col )
306     i++;
307     col_ptr[curr_col] = i;
308     }
309     col_ptr[N] = nz;
310    
311     /* create F_SMP and F_SM */
312 ksteube 1312 dist[0]=0;
313     dist[1]=N;
314     output_dist=Paso_Distribution_alloc(mpi_info, dist,1,0);
315     dist[1]=M;
316     input_dist=Paso_Distribution_alloc(mpi_info, dist,1,0);
317     mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,N,col_ptr,col_ind);
318     couplePattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,N,NULL,NULL);
319 gross 1562 send=Paso_SharedComponents_alloc(N,0,NULL,NULL,NULL,1,0,mpi_info);
320 gross 1552 connector=Paso_Connector_alloc(send,send);
321 ksteube 1312 pattern=Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,output_dist,input_dist,
322 gross 1552 mainPattern,couplePattern,couplePattern,connector,connector);
323 ksteube 1312 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_CSC, pattern, 1, 1);
324 jgs 150 /* copy values and cleanup temps */
325     for( i=0; i<nz; i++ )
326 ksteube 1312 out->mainBlock->val[i] = val[i];
327 jgs 150
328 ksteube 1312 Paso_SystemMatrixPattern_free(pattern);
329     Paso_Pattern_free(mainPattern);
330     Paso_Pattern_free(couplePattern);
331 gross 1552 Paso_Connector_free(connector);
332 ksteube 1312 Paso_Distribution_free(output_dist);
333     Paso_Distribution_free(input_dist);
334     Paso_SharedComponents_free(send);
335     Paso_MPIInfo_free(mpi_info);
336 jgs 150 MEMFREE( val );
337     MEMFREE( col_ind );
338     return out;
339     }

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26