1 |
|
2 |
/******************************************************* |
3 |
* |
4 |
* Copyright (c) 2003-2008 by University of Queensland |
5 |
* Earth Systems Science Computational Center (ESSCC) |
6 |
* http://www.uq.edu.au/esscc |
7 |
* |
8 |
* Primary Business: Queensland, Australia |
9 |
* Licensed under the Open Software License version 3.0 |
10 |
* http://www.opensource.org/licenses/osl-3.0.php |
11 |
* |
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 |
#define FSCANF_CHECK(scan_ret, reason) { if (scan_ret == EOF) perror(reason); return NULL; } |
31 |
|
32 |
static void swap( index_t*, index_t*, double*, int, int ); |
33 |
static void q_sort( index_t*, index_t*, double*, int, int ); |
34 |
static void print_entries( index_t*, index_t*, double* ); |
35 |
|
36 |
static int M, N, nz; |
37 |
|
38 |
|
39 |
/* debug: print the entries */ |
40 |
void print_entries( index_t *r, index_t *c, double *v ) |
41 |
{ |
42 |
int i; |
43 |
|
44 |
for( i=0; i<nz; i++ ) |
45 |
{ |
46 |
printf( "(%ld, %ld) == %e\n", (long)r[i], (long)c[i], v[i] ); |
47 |
} |
48 |
} |
49 |
|
50 |
/* swap function */ |
51 |
void swap( index_t *r, index_t *c, double *v, int left, int right ) |
52 |
{ |
53 |
double v_temp; |
54 |
index_t temp; |
55 |
|
56 |
temp = r[left]; |
57 |
r[left] = r[right]; |
58 |
r[right] = temp; |
59 |
|
60 |
temp = c[left]; |
61 |
c[left] = c[right]; |
62 |
c[right] = temp; |
63 |
|
64 |
v_temp = v[left]; |
65 |
v[left] = v[right]; |
66 |
v[right] = v_temp; |
67 |
} |
68 |
|
69 |
void q_sort( index_t *row, index_t *col, double *val, int begin, int end ) |
70 |
{ |
71 |
int l, r; |
72 |
index_t pivot, lval; |
73 |
|
74 |
if( end > begin ) |
75 |
{ |
76 |
pivot = N * row[begin] + col[begin]; |
77 |
l = begin + 1; |
78 |
r = end; |
79 |
|
80 |
while( l < r ) |
81 |
{ |
82 |
lval = N * row[l] + col[l]; |
83 |
if( lval < pivot ) |
84 |
l++; |
85 |
else |
86 |
{ |
87 |
r--; |
88 |
swap( row, col, val, l, r ); |
89 |
} |
90 |
} |
91 |
l--; |
92 |
swap( row, col, val, begin, l ); |
93 |
q_sort( row, col, val, begin, l ); |
94 |
q_sort( row, col, val, r, end ); |
95 |
} |
96 |
} |
97 |
|
98 |
Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSR( char *fileName_p ) |
99 |
{ |
100 |
index_t dist[2]; |
101 |
Paso_Distribution* input_dist=NULL, *output_dist=NULL; |
102 |
index_t *col_ind = NULL; |
103 |
index_t *row_ind = NULL; |
104 |
index_t *row_ptr = NULL; |
105 |
double *val = NULL; |
106 |
FILE *fileHandle_p = NULL; |
107 |
Paso_Pattern* mainPattern=NULL, *couplePattern=NULL; |
108 |
Paso_SystemMatrixPattern *pattern = NULL; |
109 |
Paso_SystemMatrix *out = NULL; |
110 |
Paso_SharedComponents *send =NULL; |
111 |
Paso_Connector *connector=NULL; |
112 |
int i, curr_row, scan_ret; |
113 |
MM_typecode matrixCode; |
114 |
Paso_MPIInfo* mpi_info=Paso_MPIInfo_alloc( MPI_COMM_WORLD); |
115 |
Paso_resetError(); |
116 |
if (mpi_info->size >1) { |
117 |
Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: support single processor only"); |
118 |
return NULL; |
119 |
} |
120 |
/* open the file */ |
121 |
fileHandle_p = fopen( fileName_p, "r" ); |
122 |
if( fileHandle_p == NULL ) |
123 |
{ |
124 |
Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Cannot read file for reading."); |
125 |
Paso_MPIInfo_free(mpi_info); |
126 |
return NULL; |
127 |
} |
128 |
|
129 |
/* process banner */ |
130 |
if( mm_read_banner(fileHandle_p, &matrixCode) != 0 ) |
131 |
{ |
132 |
Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Error processing MM banner."); |
133 |
Paso_MPIInfo_free(mpi_info); |
134 |
fclose( fileHandle_p ); |
135 |
return NULL; |
136 |
} |
137 |
if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) ) |
138 |
{ |
139 |
|
140 |
Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSR: found Matrix Market type is not supported."); |
141 |
Paso_MPIInfo_free(mpi_info); |
142 |
fclose( fileHandle_p ); |
143 |
return NULL; |
144 |
} |
145 |
|
146 |
/* get matrix size */ |
147 |
if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 ) |
148 |
{ |
149 |
Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not parse matrix size"); |
150 |
Paso_MPIInfo_free(mpi_info); |
151 |
fclose( fileHandle_p ); |
152 |
return NULL; |
153 |
} |
154 |
|
155 |
/* prepare storage */ |
156 |
col_ind = MEMALLOC( nz, index_t ); |
157 |
row_ind = MEMALLOC( nz, index_t ); |
158 |
val = MEMALLOC( nz, double ); |
159 |
|
160 |
row_ptr = MEMALLOC( (M+1), index_t ); |
161 |
|
162 |
if( col_ind == NULL || row_ind == NULL || val == NULL || row_ptr == NULL ) |
163 |
{ |
164 |
Paso_setError(MEMORY_ERROR, "Paso_SystemMatrix_loadMM_toCSR: Could not allocate memory" ); |
165 |
|
166 |
Paso_MPIInfo_free(mpi_info); |
167 |
fclose( fileHandle_p ); |
168 |
return NULL; |
169 |
} |
170 |
|
171 |
/* perform actual read of elements */ |
172 |
for( i=0; i<nz; i++ ) |
173 |
{ |
174 |
scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] ); |
175 |
FSCANF_CHECK(scan_ret, "fscanf: read elements 1") |
176 |
row_ind[i]--; |
177 |
col_ind[i]--; |
178 |
} |
179 |
fclose( fileHandle_p ); |
180 |
|
181 |
/* sort the entries */ |
182 |
q_sort( row_ind, col_ind, val, 0, nz ); |
183 |
|
184 |
/* setup row_ptr */ |
185 |
curr_row = 0; |
186 |
for( i=0; (i<nz && curr_row<M); curr_row++ ) |
187 |
{ |
188 |
while( row_ind[i] != curr_row ) |
189 |
i++; |
190 |
row_ptr[curr_row] = i; |
191 |
} |
192 |
row_ptr[M] = nz; |
193 |
|
194 |
/* create return value */ |
195 |
/* create F_SMP and F_SM */ |
196 |
dist[0]=0; |
197 |
dist[1]=M; |
198 |
output_dist=Paso_Distribution_alloc(mpi_info, dist,1,0); |
199 |
dist[1]=N; |
200 |
input_dist=Paso_Distribution_alloc(mpi_info, dist,1,0); |
201 |
mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,N,row_ptr,row_ind); |
202 |
couplePattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,M,N,NULL,NULL); |
203 |
send=Paso_SharedComponents_alloc(M,0,NULL,NULL,NULL,1,0,mpi_info); |
204 |
connector=Paso_Connector_alloc(send,send); |
205 |
pattern=Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,output_dist,input_dist, |
206 |
mainPattern,couplePattern,couplePattern,connector,connector); |
207 |
|
208 |
out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_DEFAULT, pattern, 1, 1); |
209 |
/* copy values and cleanup temps */ |
210 |
for( i=0; i<nz; i++ ) out->mainBlock->val[i] = val[i]; |
211 |
|
212 |
Paso_SystemMatrixPattern_free(pattern); |
213 |
Paso_Pattern_free(mainPattern); |
214 |
Paso_Pattern_free(couplePattern); |
215 |
Paso_Connector_free(connector); |
216 |
Paso_Distribution_free(output_dist); |
217 |
Paso_Distribution_free(input_dist); |
218 |
Paso_SharedComponents_free(send); |
219 |
Paso_MPIInfo_free(mpi_info); |
220 |
MEMFREE( val ); |
221 |
MEMFREE( row_ind ); |
222 |
|
223 |
return out; |
224 |
} |
225 |
|
226 |
Paso_SystemMatrix* Paso_SystemMatrix_loadMM_toCSC( char *fileName_p ) |
227 |
{ |
228 |
index_t dist[2]; |
229 |
Paso_Distribution* input_dist=NULL, *output_dist=NULL; |
230 |
FILE *fileHandle_p = NULL; |
231 |
Paso_Pattern* mainPattern=NULL, *couplePattern=NULL; |
232 |
Paso_SystemMatrixPattern *pattern = NULL; |
233 |
Paso_SystemMatrix *out = NULL; |
234 |
Paso_SharedComponents *send =NULL; |
235 |
Paso_Connector *connector=NULL; |
236 |
index_t *col_ind = NULL; |
237 |
index_t *row_ind = NULL; |
238 |
index_t *col_ptr = NULL; |
239 |
double *val = NULL; |
240 |
int i, curr_col=0, scan_ret; |
241 |
MM_typecode matrixCode; |
242 |
Paso_MPIInfo* mpi_info=Paso_MPIInfo_alloc( MPI_COMM_WORLD); |
243 |
if (mpi_info->size >1) { |
244 |
Paso_setError(IO_ERROR, "Paso_SystemMatrix_loadMM_toCSC: support single processor only"); |
245 |
return NULL; |
246 |
} |
247 |
|
248 |
Paso_resetError(); |
249 |
|
250 |
/* open the file */ |
251 |
fileHandle_p = fopen( fileName_p, "r" ); |
252 |
if( fileHandle_p == NULL ) |
253 |
{ |
254 |
Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: File could not be opened for reading"); |
255 |
Paso_MPIInfo_free(mpi_info); |
256 |
return NULL; |
257 |
} |
258 |
|
259 |
/* process banner */ |
260 |
if( mm_read_banner(fileHandle_p, &matrixCode) != 0 ) |
261 |
{ |
262 |
Paso_setError(IO_ERROR,"Paso_SystemMatrix_loadMM_toCSC: Error processing MM banner"); |
263 |
fclose( fileHandle_p ); |
264 |
Paso_MPIInfo_free(mpi_info); |
265 |
return NULL; |
266 |
} |
267 |
if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) ) |
268 |
{ |
269 |
Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported."); |
270 |
fclose( fileHandle_p ); |
271 |
Paso_MPIInfo_free(mpi_info); |
272 |
return NULL; |
273 |
} |
274 |
|
275 |
/* get matrix size */ |
276 |
if( mm_read_mtx_crd_size(fileHandle_p, &M, &N, &nz) != 0 ) |
277 |
{ |
278 |
Paso_setError(TYPE_ERROR,"Paso_SystemMatrix_loadMM_toCSC: found Matrix Market type is not supported."); |
279 |
fclose( fileHandle_p ); |
280 |
Paso_MPIInfo_free(mpi_info); |
281 |
return NULL; |
282 |
} |
283 |
|
284 |
/* prepare storage */ |
285 |
col_ind = MEMALLOC( nz, index_t ); |
286 |
row_ind = MEMALLOC( nz, index_t ); |
287 |
val = MEMALLOC( nz, double ); |
288 |
|
289 |
col_ptr = MEMALLOC( (N+1), index_t ); |
290 |
|
291 |
|
292 |
/* perform actual read of elements */ |
293 |
for( i=0; i<nz; i++ ) |
294 |
{ |
295 |
scan_ret = fscanf( fileHandle_p, "%d %d %le\n", &row_ind[i], &col_ind[i], &val[i] ); |
296 |
FSCANF_CHECK(scan_ret, "fscanf: read elements 2") |
297 |
row_ind[i]--; |
298 |
col_ind[i]--; |
299 |
} |
300 |
fclose( fileHandle_p ); |
301 |
|
302 |
/* sort the entries */ |
303 |
q_sort( col_ind, row_ind, val, 0, nz ); |
304 |
|
305 |
/* setup row_ptr */ |
306 |
for( i=0; (i<nz && curr_col<N); curr_col++ ) |
307 |
{ |
308 |
while( col_ind[i] != curr_col ) |
309 |
i++; |
310 |
col_ptr[curr_col] = i; |
311 |
} |
312 |
col_ptr[N] = nz; |
313 |
|
314 |
/* create F_SMP and F_SM */ |
315 |
dist[0]=0; |
316 |
dist[1]=N; |
317 |
output_dist=Paso_Distribution_alloc(mpi_info, dist,1,0); |
318 |
dist[1]=M; |
319 |
input_dist=Paso_Distribution_alloc(mpi_info, dist,1,0); |
320 |
mainPattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,N,M,col_ptr,col_ind); |
321 |
couplePattern=Paso_Pattern_alloc(PATTERN_FORMAT_DEFAULT,1,1,N,M,NULL,NULL); |
322 |
send=Paso_SharedComponents_alloc(N,0,NULL,NULL,NULL,1,0,mpi_info); |
323 |
connector=Paso_Connector_alloc(send,send); |
324 |
pattern=Paso_SystemMatrixPattern_alloc(PATTERN_FORMAT_DEFAULT,output_dist,input_dist, |
325 |
mainPattern,couplePattern,couplePattern,connector,connector); |
326 |
out = Paso_SystemMatrix_alloc(MATRIX_FORMAT_CSC, pattern, 1, 1); |
327 |
/* copy values and cleanup temps */ |
328 |
for( i=0; i<nz; i++ ) |
329 |
out->mainBlock->val[i] = val[i]; |
330 |
|
331 |
Paso_SystemMatrixPattern_free(pattern); |
332 |
Paso_Pattern_free(mainPattern); |
333 |
Paso_Pattern_free(couplePattern); |
334 |
Paso_Connector_free(connector); |
335 |
Paso_Distribution_free(output_dist); |
336 |
Paso_Distribution_free(input_dist); |
337 |
Paso_SharedComponents_free(send); |
338 |
Paso_MPIInfo_free(mpi_info); |
339 |
MEMFREE( val ); |
340 |
MEMFREE( col_ind ); |
341 |
return out; |
342 |
} |