/[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 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years ago) by ksteube
File MIME type: text/plain
File size: 9691 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

1
2 /* $Id$ */
3
4 /*******************************************************
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
16 /**************************************************************/
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 index_t dist[2];
100 Paso_Distribution* input_dist=NULL, *output_dist=NULL;
101 index_t *col_ind = NULL;
102 index_t *row_ind = NULL;
103 index_t *row_ptr = NULL;
104 double *val = NULL;
105 FILE *fileHandle_p = NULL;
106 Paso_Pattern* mainPattern=NULL, *couplePattern=NULL;
107 Paso_SystemMatrixPattern *pattern = NULL;
108 Paso_SystemMatrix *out = NULL;
109 Paso_SharedComponents *send =NULL;
110 Paso_Coupler *coupler=NULL;
111 int i, curr_row;
112 MM_typecode matrixCode;
113 Paso_MPIInfo* mpi_info=Paso_MPIInfo_alloc( MPI_COMM_WORLD);
114 Paso_resetError();
115 if (mpi_info->size >1) {
116 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: support single processor only");
117 return NULL;
118 }
119 /* open the file */
120 fileHandle_p = fopen( fileName_p, "r" );
121 if( fileHandle_p == NULL )
122 {
123 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Cannot read file for reading.");
124 Paso_MPIInfo_free(mpi_info);
125 return NULL;
126 }
127
128 /* process banner */
129 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
130 {
131 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Error processing MM banner.");
132 Paso_MPIInfo_free(mpi_info);
133 fclose( fileHandle_p );
134 return NULL;
135 }
136 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
137 {
138
139 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSR: found Matrix Market type is not supported.");
140 Paso_MPIInfo_free(mpi_info);
141 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 Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not parse matrix size");
149 Paso_MPIInfo_free(mpi_info);
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 Paso_setError(MEMORY_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not allocate memory" );
164
165 Paso_MPIInfo_free(mpi_info);
166 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 /* create return value */
193 /* create F_SMP and F_SM */
194 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 send=Paso_SharedComponents_alloc(0,NULL,NULL,NULL,1,0,mpi_info);
202 coupler=Paso_Coupler_alloc(send,send);
203 pattern=Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,output_dist,input_dist,
204 mainPattern,couplePattern,coupler);
205
206 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_DEFAULT, pattern, 1, 1);
207 /* copy values and cleanup temps */
208 for( i=0; i<nz; i++ ) out->mainBlock->val[i] = val[i];
209
210 Paso_SystemMatrixPattern_free(pattern);
211 Paso_Pattern_free(mainPattern);
212 Paso_Pattern_free(couplePattern);
213 Paso_Coupler_free(coupler);
214 Paso_Distribution_free(output_dist);
215 Paso_Distribution_free(input_dist);
216 Paso_SharedComponents_free(send);
217 Paso_MPIInfo_free(mpi_info);
218 MEMFREE( val );
219 MEMFREE( row_ind );
220
221 return out;
222 }
223
224 Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSC( char *fileName_p )
225 {
226 index_t dist[2];
227 Paso_Distribution* input_dist=NULL, *output_dist=NULL;
228 FILE *fileHandle_p = NULL;
229 Paso_Pattern* mainPattern=NULL, *couplePattern=NULL;
230 Paso_SystemMatrixPattern *pattern = NULL;
231 Paso_SystemMatrix *out = NULL;
232 Paso_SharedComponents *send =NULL;
233 Paso_Coupler *coupler=NULL;
234 index_t *col_ind = NULL;
235 index_t *row_ind = NULL;
236 index_t *col_ptr = NULL;
237 double *val = NULL;
238 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
246 Paso_resetError();
247
248 /* open the file */
249 fileHandle_p = fopen( fileName_p, "r" );
250 if( fileHandle_p == NULL )
251 {
252 Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: File could not be opened for reading");
253 Paso_MPIInfo_free(mpi_info);
254 return NULL;
255 }
256
257 /* process banner */
258 if( mm_read_banner(fileHandle_p, &matrixCode) != 0 )
259 {
260 Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: Error processing MM banner");
261 fclose( fileHandle_p );
262 Paso_MPIInfo_free(mpi_info);
263 return NULL;
264 }
265 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) )
266 {
267 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
268 fclose( fileHandle_p );
269 Paso_MPIInfo_free(mpi_info);
270 return NULL;
271 }
272
273 /* get matrix size */
274 if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 )
275 {
276 Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported.");
277 fclose( fileHandle_p );
278 Paso_MPIInfo_free(mpi_info);
279 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 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 send=Paso_SharedComponents_alloc(0,NULL,NULL,NULL,1,0,mpi_info);
320 coupler=Paso_Coupler_alloc(send,send);
321 pattern=Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,output_dist,input_dist,
322 mainPattern,couplePattern,coupler);
323 out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_CSC, pattern, 1, 1);
324 /* copy values and cleanup temps */
325 for( i=0; i<nz; i++ )
326 out->mainBlock->val[i] = val[i];
327
328 Paso_SystemMatrixPattern_free(pattern);
329 Paso_Pattern_free(mainPattern);
330 Paso_Pattern_free(couplePattern);
331 Paso_Coupler_free(coupler);
332 Paso_Distribution_free(output_dist);
333 Paso_Distribution_free(input_dist);
334 Paso_SharedComponents_free(send);
335 Paso_MPIInfo_free(mpi_info);
336 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