/[escript]/branches/doubleplusgood/paso/src/SystemMatrix_loadMM.cpp
ViewVC logotype

Contents of /branches/doubleplusgood/paso/src/SystemMatrix_loadMM.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4280 - (show annotations)
Wed Mar 6 06:45:32 2013 UTC (6 years, 1 month ago) by jfenwick
File size: 12236 byte(s)
some memory management replacement
1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2013 by University of Queensland
5 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
10 *
11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 * Development since 2012 by School of Earth Sciences
13 *
14 *****************************************************************************/
15
16
17 /************************************************************************************/
18
19 /* Paso: Matrix Market format is loaded to a SystemMatrix */
20
21 /************************************************************************************/
22
23 /* Copyrights by ACcESS Australia 2003,2004,2005 */
24 /* Author: imran@access.edu.au */
25
26 /************************************************************************************/
27
28 #include "Paso.h"
29 #include "mmio.h"
30 #include "SystemMatrix.h"
31
32 #include "limits.h"
33
34 #define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) perror(reason); return NULL; }
35
36 static void swap( index_t*, index_t*, double*, int, int );
37 static void q_sort( index_t*, index_t*, double*, int, int );
38 /*static void print_entries( index_t*, index_t*, double* );*/
39
40 static int M, N, nz;
41
42
43 /* debug: print the entries */
44 /*
45 void print_entries( index_t *r, index_t *c, double *v )
46 {
47 int i;
48
49 for( i=0; i<nz; i++ )
50 {
51 printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] );
52 }
53 }
54 */
55
56 /* swap function */
57 void swap( index_t *r, index_t *c, double *v, int left, int right )
58 {
59 double v_temp;
60 index_t temp;
61
62 temp = r[left];
63 r[left] = r[right];
64 r[right] = temp;
65
66 temp = c[left];
67 c[left] = c[right];
68 c[right] = temp;
69
70 v_temp = v[left];
71 v[left] = v[right];
72 v[right] = v_temp;
73 }
74
75 void q_sort( index_t *row, index_t *col, double *val, int begin, int end )
76 {
77 int l, r;
78 int flag;
79
80 if( end > begin )
81 {
82 l = begin + 1;
83 r = end;
84
85 while( l < r )
86 {
87 /* This whole section is for checking lval<pivot, where
88 pivot=N*row[begin]+col[begin] and lval=N*row[l]+col[l]. */
89 if (row[l]<row[begin])
90 {
91 if (ABS(row[l]-row[begin])==1 && ABS(col[l]-col[begin])==N)
92 flag=0;
93 else
94 flag=1;
95 }
96 else if (row[l]==row[begin])
97 if (col[l]<col[begin])
98 flag=1;
99 else
100 flag=0;
101 else {
102 if (ABS(row[l]-row[begin])==1 && ABS(col[l]-col[begin])==N)
103 flag=1;
104 else
105 flag=0;
106 }
107
108
109 if(flag==1)
110 l++;
111 else {
112 r--;
113 swap( row, col, val, l, r );
114 }
115
116 }
117 l--;
118 swap( row, col, val, begin, l );
119 q_sort( row, col, val, begin, l );
120 q_sort( row, col, val, r, end );
121 }
122 }
123
124 Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSR( char *fileName_p )
125 {
126 index_t dist[2];
127 Paso_Distribution* input_dist=NULL, *output_dist=NULL;
128 index_t *col_ind = NULL;
129 index_t *row_ind = NULL;
130 index_t *row_ptr = NULL;
131 double *val = NULL;
132 FILE *fileHandle_p = NULL;
133 Paso_Pattern* mainPattern=NULL, *couplePattern=NULL;
134 Paso_SystemMatrixPattern *pattern = NULL;
135 Paso_SystemMatrix *out = NULL;
136 Paso_SharedComponents *send =NULL;
137 Paso_Connector *connector=NULL;
138 int i, curr_row, scan_ret;
139 MM_typecode matrixCode;
140 Esys_MPIInfo* mpi_info=Esys_MPIInfo_alloc( MPI_COMM_WORLD);
141 Esys_resetError();
142 if (mpi_info->size >1) {
143 Esys_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: supports single processor only");
144 return NULL;
145 }
146 /* open the file */
147 fileHandle_p = fopen( fileName_p, "r" );
148 if( fileHandle_p == NULL )
149 {
150 Esys_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Cannot open file for reading.");
151 Esys_MPIInfo_free(mpi_info);
152 return NULL;
153 }
154
155 /* process banner */
156 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
157 {
158 Esys_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Error processing MM banner.");
159 Esys_MPIInfo_free(mpi_info);
160 fclose( fileHandle_p );
161 return NULL;
162 }
163 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
164 {
165
166 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
167 Esys_MPIInfo_free(mpi_info);
168 fclose( fileHandle_p );
169 return NULL;
170 }
171
172 /* get matrix size */
173 if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
174 {
175 Esys_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not read sparse matrix size.");
176 Esys_MPIInfo_free(mpi_info);
177 fclose( fileHandle_p );
178 return NULL;
179 }
180
181 /* prepare storage */
182 col_ind = new index_t [ nz];
183 row_ind = new index_t [ nz];
184 val = new double [ nz];
185
186 row_ptr = new index_t [ (M+1)];
187
188 if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL )
189 {
190 Esys_setError(MEMORY_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not allocate memory.");
191
192 Esys_MPIInfo_free(mpi_info);
193 fclose( fileHandle_p );
194 return NULL;
195 }
196
197 /* perform actual read of elements */
198 for( i=0; i<nz; i++ )
199 {
200 scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
201 if (scan_ret!=3)
202 {
203 delete[] val ;
204 delete[] row_ind ;
205 delete[] col_ind ;
206 delete[] row_ptr ;
207 Esys_MPIInfo_free(mpi_info);
208 fclose(fileHandle_p);
209 return NULL;
210 }
211 row_ind[i]--;
212 col_ind[i]--;
213 }
214 fclose( fileHandle_p );
215 /* sort the entries */
216 q_sort( row_ind, col_ind, val, 0, nz );
217
218 /* setup row_ptr */
219 curr_row = 0;
220 for( i=0; (i<nz && curr_row<M); curr_row++ )
221 {
222 while( row_ind[i] != curr_row ){
223 i++;
224 }
225 row_ptr[curr_row] = i;
226 }
227 row_ptr[M] = nz;
228
229 /* create return value */
230 /* create F_SMP and F_SM */
231 dist[0]=0;
232 dist[1]=M;
233 output_dist=Paso_Distribution_alloc(mpi_info, dist,1,0);
234 dist[1]=N;
235 input_dist=Paso_Distribution_alloc(mpi_info, dist,1,0);
236 mainPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,M,N,row_ptr,col_ind);
237 couplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,M,N,NULL,NULL);
238 dist[0]=M;
239 send=Paso_SharedComponents_alloc(M,0,NULL,NULL,dist,1,0,mpi_info);
240 dist[0]=0;
241 connector=Paso_Connector_alloc(send,send);
242 pattern=Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,output_dist,input_dist,
243 mainPattern,couplePattern,couplePattern,connector,connector);
244
245 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_DEFAULT, pattern, 1, 1, TRUE);
246 /* copy values and cleanup temps */
247 for( i=0; i<nz; i++ ) out->mainBlock->val[i] = val[i];
248
249 Paso_SystemMatrixPattern_free(pattern);
250 Paso_Pattern_free(mainPattern);
251 Paso_Pattern_free(couplePattern);
252 Paso_Connector_free(connector);
253 Paso_Distribution_free(output_dist);
254 Paso_Distribution_free(input_dist);
255 Paso_SharedComponents_free(send);
256 Esys_MPIInfo_free(mpi_info);
257 delete[] val ;
258 delete[] row_ind ;
259 return out;
260 }
261
262 Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSC( char *fileName_p )
263 {
264 index_t dist[2];
265 Paso_Distribution* input_dist=NULL, *output_dist=NULL;
266 FILE *fileHandle_p = NULL;
267 Paso_Pattern* mainPattern=NULL, *couplePattern=NULL;
268 Paso_SystemMatrixPattern *pattern = NULL;
269 Paso_SystemMatrix *out = NULL;
270 Paso_SharedComponents *send =NULL;
271 Paso_Connector *connector=NULL;
272 index_t *col_ind = NULL;
273 index_t *row_ind = NULL;
274 index_t *col_ptr = NULL;
275 double *val = NULL;
276 int i, curr_col=0, scan_ret;
277 MM_typecode matrixCode;
278 Esys_MPIInfo* mpi_info=Esys_MPIInfo_alloc( MPI_COMM_WORLD);
279 if (mpi_info->size >1) {
280 Esys_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSC: supports single processor only");
281 return NULL;
282 }
283
284 Esys_resetError();
285
286 /* open the file */
287 fileHandle_p = fopen( fileName_p, "r" );
288 if( fileHandle_p == NULL )
289 {
290 Esys_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: File could not be opened for reading.");
291 Esys_MPIInfo_free(mpi_info);
292 return NULL;
293 }
294
295 /* process banner */
296 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
297 {
298 Esys_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: Error processing MM banner.");
299 fclose( fileHandle_p );
300 Esys_MPIInfo_free(mpi_info);
301 return NULL;
302 }
303 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
304 {
305 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
306 fclose( fileHandle_p );
307 Esys_MPIInfo_free(mpi_info);
308 return NULL;
309 }
310
311 /* get matrix size */
312 if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
313 {
314 Esys_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
315 fclose( fileHandle_p );
316 Esys_MPIInfo_free(mpi_info);
317 return NULL;
318 }
319
320 /* prepare storage */
321 col_ind = new index_t [ nz];
322 row_ind = new index_t [ nz];
323 val = new double [ nz];
324
325 col_ptr = new index_t [ (N+1)];
326
327
328 /* perform actual read of elements */
329 for( i=0; i<nz; i++ )
330 {
331 scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] );
332 if (scan_ret!=3)
333 {
334 delete[] val ;
335 delete[] row_ind ;
336 delete[] col_ind ;
337 delete[] col_ptr ;
338 Esys_MPIInfo_free(mpi_info);
339 fclose(fileHandle_p);
340 return NULL;
341 }
342 row_ind[i]--;
343 col_ind[i]--;
344 }
345 fclose( fileHandle_p );
346
347 /* sort the entries */
348 q_sort( col_ind, row_ind, val, 0, nz );
349
350 /* setup row_ptr */
351 for( i=0; (i<nz && curr_col<N); curr_col++ )
352 {
353 while( col_ind[i] != curr_col )
354 i++;
355 col_ptr[curr_col] = i;
356 }
357 col_ptr[N] = nz;
358
359 /* create F_SMP and F_SM */
360 dist[0]=0;
361 dist[1]=N;
362 output_dist=Paso_Distribution_alloc(mpi_info, dist,1,0);
363 dist[1]=M;
364 input_dist=Paso_Distribution_alloc(mpi_info, dist,1,0);
365 mainPattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,N,M,col_ptr,col_ind);
366 couplePattern=Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT,N,M,NULL,NULL);
367 send=Paso_SharedComponents_alloc(N,0,NULL,NULL,NULL,1,0,mpi_info);
368 connector=Paso_Connector_alloc(send,send);
369 pattern=Paso_SystemMatrixPattern_alloc(MATRIX_FORMAT_DEFAULT,output_dist,input_dist,
370 mainPattern,couplePattern,couplePattern,connector,connector);
371 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_CSC, pattern, 1, 1, TRUE);
372 /* copy values and cleanup temps */
373 for( i=0; i<nz; i++ )
374 out->mainBlock->val[i] = val[i];
375
376 Paso_SystemMatrixPattern_free(pattern);
377 Paso_Pattern_free(mainPattern);
378 Paso_Pattern_free(couplePattern);
379 Paso_Connector_free(connector);
380 Paso_Distribution_free(output_dist);
381 Paso_Distribution_free(input_dist);
382 Paso_SharedComponents_free(send);
383 Esys_MPIInfo_free(mpi_info);
384 delete[] val ;
385 delete[] row_ind ;
386 return out;
387 }
388
389 void Paso_RHS_loadMM_toCSR( char *fileName_p, double *b, dim_t size)
390 {
391 FILE *fileHandle_p = NULL;
392 int i, scan_ret;
393 MM_typecode matrixCode;
394 Esys_resetError();
395 /* open the file */
396 fileHandle_p = fopen( fileName_p, "r" );
397 if( fileHandle_p == NULL )
398 {
399 Esys_setError(IO_ERROR, "Paso_RHS_loadMM_toCSR: Cannot open file for reading.");
400 }
401
402 /* process banner */
403 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
404 {
405 Esys_setError(IO_ERROR, "Paso_RHS_loadMM_toCSR: Error processing MM banner.");
406 }
407 if( !(mm_is_real(matrixCode) && mm_is_general(matrixCode) && mm_is_array(matrixCode)) )
408 {
409
410 Esys_setError(TYPE_ERROR,"Paso_RHS_loadMM_toCSR: found Matrix Market type is not supported.");
411 }
412
413 /* get matrix size */
414 if( mm_read_mtx_array_size(fileHandle_p, &M, &N) != 0 )
415 {
416 Esys_setError(IO_ERROR, "Paso_RHS_loadMM_toCSR: Could not read sparse matrix size.");
417 }
418
419 if(M!=size){
420 Esys_setError(IO_ERROR, "Paso_RHS_loadMM_toCSR: Actual and provided sizes do not match.");
421 }
422
423 if (Esys_noError()) {
424 nz=M;
425 /* perform actual read of elements */
426 for( i=0; i<nz; i++ )
427 {
428 scan_ret = fscanf( fileHandle_p, "%le\n", &b[i] );
429 if (scan_ret!=1)
430 {
431 fclose(fileHandle_p);
432 Esys_setError(IO_ERROR, "Paso_RHS_loadMM_toCSR: Could not read some of the values.");
433 }
434 }
435 }
436 else {
437 fclose( fileHandle_p );
438 }
439
440 }
441

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26