/[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 1736 - (show annotations)
Fri Aug 29 02:23:16 2008 UTC (11 years, 2 months ago) by gross
File MIME type: text/plain
File size: 9779 byte(s)
This fixes a problem which is typically arising when using reduced order
with MPI and a "small" number of elements per processor. In this case it
can happen that the couple matrix is not using all entries sent to the
processor. The old implementations assumed that the indices will cover
the entire input. This assumption has been removed.


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_Connector *connector=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,N,row_ptr,row_ind);
200 couplePattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,N,NULL,NULL);
201 send=Paso_SharedComponents_alloc(M,0,NULL,NULL,NULL,1,0,mpi_info);
202 connector=Paso_Connector_alloc(send,send);
203 pattern=Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,output_dist,input_dist,
204 mainPattern,couplePattern,couplePattern,connector,connector);
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_Connector_free(connector);
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_Connector *connector=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,M,col_ptr,col_ind);
318 couplePattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,N,M,NULL,NULL);
319 send=Paso_SharedComponents_alloc(N,0,NULL,NULL,NULL,1,0,mpi_info);
320 connector=Paso_Connector_alloc(send,send);
321 pattern=Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,output_dist,input_dist,
322 mainPattern,couplePattern,couplePattern,connector,connector);
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_Connector_free(connector);
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