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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 969 - (show annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 8 months ago) by ksteube
File MIME type: text/plain
File size: 8279 byte(s)
Parallelization using MPI for solution of implicit problems.

Parallelization for explicit problems has already been accomplished in
the main SVN branch.

This is incomplete and is not ready for use.


1 /* $Id$ */
2
3
4 /*
5 ********************************************************************************
6 * Copyright 2006 by ACcESS MNRF *
7 * *
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 /**************************************************************/
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 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Cannot read file for reading.");
115 return NULL;
116 }
117
118 /* process banner */
119 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
120 {
121 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Error processing MM banner.");
122 fclose( fileHandle_p );
123 return NULL;
124 }
125 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
126 {
127
128 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
129 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 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not parse matrix size");
137 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 Paso_setError(MEMORY_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not allocate memory" );
151
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 loc_pattern = Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT, M, row_ptr, col_ind, NULL ); /* NULL should be an MPIInfo if this is ever to be used */
180 if(! Paso_noError() )
181 return NULL;
182
183 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_DEFAULT, loc_pattern, 1, 1 );
184 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 Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: File could not be opened for reading");
220 return NULL;
221 }
222
223 /* process banner */
224 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
225 {
226 Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: Error processing MM banner");
227 fclose( fileHandle_p );
228 return NULL;
229 }
230 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
231 {
232 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: 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(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
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 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 // loc_pattern = Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT, M, row_ptr, col_ind, NULL ); /* NULL should be an MPIInfo if this is ever to be used */
292 loc_pattern = Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT, N, col_ptr, row_ind, NULL ); /* NULL should be an MPIInfo if this is ever to be used */
293 if(! Paso_noError() )
294 return NULL;
295
296 // out = Paso_SystemMatrix_alloc( MATRIX_FORMAT_DEFAULT, loc_pattern, 1, 1 );
297 out = Paso_SystemMatrix_alloc( MATRIX_FORMAT_CSC, loc_pattern, 1, 1 );
298 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