/[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 971 - (hide annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 7 months ago) by ksteube
File MIME type: text/plain
File size: 8084 byte(s)
Had to undo commit to new MPI branch. The changes went into the original and
not the branch. The files committed here are exactly the same as revision 969.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26