/[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 415 - (hide annotations)
Wed Jan 4 05:37:33 2006 UTC (13 years, 9 months ago) by gross
File MIME type: text/plain
File size: 7438 byte(s)
a better way of representing the matrix format type is introduced. this is needed for the Paradiso and UMFPACK interface
1 jgs 150 /* $Id$ */
2    
3     /**************************************************************/
4    
5     /* Paso: Matrix Market format is loaded to a SystemMatrix */
6    
7     /**************************************************************/
8    
9     /* Copyrights by ACcESS Australia 2003,2004,2005 */
10     /* Author: imran@access.edu.au */
11    
12     /**************************************************************/
13    
14     #include "Paso.h"
15     #include "mmio.h"
16     #include "SystemMatrix.h"
17    
18     static void swap( index_t*, index_t*, double*, int, int );
19     static void q_sort( index_t*, index_t*, double*, int, int );
20     static void print_entries( index_t*, index_t*, double* );
21    
22     static int M, N, nz;
23    
24    
25     /* debug: print the entries */
26     void print_entries( index_t *r, index_t *c, double *v )
27     {
28     int i;
29    
30     for( i=0; i<nz; i++ )
31     {
32     printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] );
33     }
34     }
35    
36     /* swap function */
37     void swap( index_t *r, index_t *c, double *v, int left, int right )
38     {
39     double v_temp;
40     index_t temp;
41    
42     temp = r[left];
43     r[left] = r[right];
44     r[right] = temp;
45    
46     temp = c[left];
47     c[left] = c[right];
48     c[right] = temp;
49    
50     v_temp = v[left];
51     v[left] = v[right];
52     v[right] = v_temp;
53     }
54    
55     void q_sort( index_t *row, index_t *col, double *val, int begin, int end )
56     {
57     int l, r;
58     index_t pivot, lval;
59    
60     if( end > begin )
61     {
62     pivot = N * row[begin] + col[begin];
63     l = begin + 1;
64     r = end;
65    
66     while( l < r )
67     {
68     lval = N * row[l] + col[l];
69     if( lval < pivot )
70     l++;
71     else
72     {
73     r--;
74     swap( row, col, val, l, r );
75     }
76     }
77     l--;
78     swap( row, col, val, begin, l );
79     q_sort( row, col, val, begin, l );
80     q_sort( row, col, val, r, end );
81     }
82     }
83    
84     Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSR( char *fileName_p )
85     {
86     int i, curr_row;
87     MM_typecode matrixCode;
88     Paso_resetError();
89    
90     index_t *col_ind = NULL;
91     index_t *row_ind = NULL;
92     index_t *row_ptr = NULL;
93     double *val = NULL;
94    
95     Paso_SystemMatrixPattern *loc_pattern = NULL;
96     Paso_SystemMatrix *out = NULL;
97    
98     /* open the file */
99     FILE *fileHandle_p = fopen( fileName_p, "r" );
100     if( fileHandle_p == NULL )
101     {
102 gross 415 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Cannot read file for reading.");
103 jgs 150 return NULL;
104     }
105    
106     /* process banner */
107     if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
108     {
109 gross 415 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Error processing MM banner.");
110 jgs 150 fclose( fileHandle_p );
111     return NULL;
112     }
113     if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
114     {
115 gross 415
116     Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
117 jgs 150 fclose( fileHandle_p );
118     return NULL;
119     }
120    
121     /* get matrix size */
122     if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
123     {
124 gross 415 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not parse matrix size");
125 jgs 150 fclose( fileHandle_p );
126     return NULL;
127     }
128    
129     /* prepare storage */
130     col_ind = MEMALLOC( nz, index_t );
131     row_ind = MEMALLOC( nz, index_t );
132     val = MEMALLOC( nz, double );
133    
134     row_ptr = MEMALLOC( (M+1), index_t );
135    
136     if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
137     {
138 gross 415 Paso_setError(MEMORY_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not allocate memory" );
139 jgs 150
140     fclose( fileHandle_p );
141     return NULL;
142     }
143    
144     /* perform actual read of elements */
145     for( i=0; i<nz; i++ )
146     {
147     fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
148     row_ind[i]--;
149     col_ind[i]--;
150     }
151     fclose( fileHandle_p );
152    
153     /* sort the entries */
154     q_sort( row_ind, col_ind, val, 0, nz );
155    
156     /* setup row_ptr */
157     curr_row = 0;
158     for( i=0; (i<nz && curr_row<M); curr_row++ )
159     {
160     while( row_ind[i] != curr_row )
161     i++;
162     row_ptr[curr_row] = i;
163     }
164     row_ptr[M] = nz;
165    
166     /* create F_SMP and F_SM */
167 gross 415 loc_pattern = Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT, M, row_ptr, col_ind );
168 jgs 150 if(! Paso_noError() )
169     return NULL;
170    
171 gross 415 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_DEFAULT, loc_pattern, 1, 1 );
172 jgs 150 if(! Paso_noError() )
173     return NULL;
174    
175     /* copy values and cleanup temps */
176     for( i=0; i<nz; i++ )
177     out->val[i] = val[i];
178    
179     MEMFREE( val );
180     MEMFREE( row_ind );
181    
182     return out;
183     }
184    
185     Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSC( char *fileName_p )
186     {
187     int i;
188     // int curr_row;
189     int curr_col;
190     MM_typecode matrixCode;
191    
192     index_t *col_ind = NULL;
193     index_t *row_ind = NULL;
194     // index_t *row_ptr = NULL;
195     index_t *col_ptr = NULL;
196     double *val = NULL;
197    
198     Paso_SystemMatrixPattern *loc_pattern = NULL;
199     Paso_SystemMatrix *out = NULL;
200    
201     Paso_resetError();
202    
203     /* open the file */
204     FILE *fileHandle_p = fopen( fileName_p, "r" );
205     if( fileHandle_p == NULL )
206     {
207 gross 415 Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: File could not be opened for reading");
208 jgs 150 return NULL;
209     }
210    
211     /* process banner */
212     if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
213     {
214 gross 415 Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: Error processing MM banner");
215 jgs 150 fclose( fileHandle_p );
216     return NULL;
217     }
218     if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
219     {
220 gross 415 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
221 jgs 150 fclose( fileHandle_p );
222     return NULL;
223     }
224    
225     /* get matrix size */
226     if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
227     {
228 gross 415 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
229 jgs 150 fclose( fileHandle_p );
230     return NULL;
231     }
232    
233     /* prepare storage */
234     col_ind = MEMALLOC( nz, index_t );
235     row_ind = MEMALLOC( nz, index_t );
236     val = MEMALLOC( nz, double );
237    
238     // row_ptr = MEMALLOC( (M+1), index_t );
239     col_ptr = MEMALLOC( (N+1), index_t );
240    
241    
242     // if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
243     if( col_ind == NULL || row_ind == NULL || val == NULL || col_ptr == NULL )
244     {
245     Paso_setError(MEMORY_ERROR,"Could not allocate memory" );
246     fclose( fileHandle_p );
247     return NULL;
248     }
249    
250     /* perform actual read of elements */
251     for( i=0; i<nz; i++ )
252     {
253     fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
254     row_ind[i]--;
255     col_ind[i]--;
256     }
257     fclose( fileHandle_p );
258    
259     /* sort the entries */
260     // q_sort( row_ind, col_ind, val, 0, nz );
261     q_sort( col_ind, row_ind, val, 0, nz );
262    
263     /* setup row_ptr */
264     // curr_row = 0;
265     curr_col = 0;
266     // for( i=0; (i<nz && curr_row<M); curr_row++ )
267     for( i=0; (i<nz && curr_col<N); curr_col++ )
268     {
269     // while( row_ind[i] != curr_row )
270     while( col_ind[i] != curr_col )
271     i++;
272     // row_ptr[curr_row] = i;
273     col_ptr[curr_col] = i;
274     }
275     // row_ptr[M] = nz;
276     col_ptr[N] = nz;
277    
278     /* create F_SMP and F_SM */
279 gross 415 // loc_pattern = Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT, M, row_ptr, col_ind );
280     loc_pattern = Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT, N, col_ptr, row_ind );
281 jgs 150 if(! Paso_noError() )
282     return NULL;
283    
284 gross 415 // out = Paso_SystemMatrix_alloc( MATRIX_FORMAT_DEFAULT, loc_pattern, 1, 1 );
285     out = Paso_SystemMatrix_alloc( MATRIX_FORMAT_CSC, loc_pattern, 1, 1 );
286 jgs 150 if(! Paso_noError() )
287     return NULL;
288    
289     /* copy values and cleanup temps */
290     for( i=0; i<nz; i++ )
291     out->val[i] = val[i];
292    
293     MEMFREE( val );
294     // MEMFREE( row_ind );
295     MEMFREE( col_ind );
296    
297     return out;
298     }
299    
300    
301     /*
302     * $Log$
303     * Revision 1.2 2005/09/15 03:44:39 jgs
304     * Merge of development branch dev-02 back to main trunk on 2005-09-15
305     *
306     * Revision 1.1.2.1 2005/09/05 06:29:48 gross
307     * These files have been extracted from finley to define a stand alone libray for iterative
308     * linear solvers on the ALTIX. main entry through Paso_solve. this version compiles but
309     * has not been tested yet.
310     *
311     *
312     */

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26