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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5933 - (show annotations)
Wed Feb 17 23:53:30 2016 UTC (21 months, 4 weeks ago) by caltinay
File size: 10680 byte(s)
sync with trunk.

1
2 /*****************************************************************************
3 *
4 * Copyright (c) 2003-2016 by The 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 2012-2013 by School of Earth Sciences
13 * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 *
15 *****************************************************************************/
16
17
18 /****************************************************************************/
19
20 /* Paso: Matrix Market format is loaded to a SystemMatrix */
21
22 /****************************************************************************/
23
24 /* Copyrights by ACcESS Australia 2003,2004,2005 */
25 /* Author: imran@access.edu.au */
26
27 /****************************************************************************/
28
29 #include "Paso.h"
30 #include "mmio.h"
31 #include "SystemMatrix.h"
32
33 #include "limits.h"
34
35 namespace paso {
36
37 static void swap(index_t*, index_t*, double*, int, int);
38 static void q_sort(index_t*, index_t*, double*, int, int);
39
40 static int M, N, nz;
41
42 void swap(index_t *r, index_t *c, double *v, int left, int right)
43 {
44 double v_temp;
45 index_t temp;
46
47 temp = r[left];
48 r[left] = r[right];
49 r[right] = temp;
50
51 temp = c[left];
52 c[left] = c[right];
53 c[right] = temp;
54
55 v_temp = v[left];
56 v[left] = v[right];
57 v[right] = v_temp;
58 }
59
60 void q_sort(index_t *row, index_t *col, double *val, int begin, int end)
61 {
62 int l, r;
63 int flag;
64
65 if (end > begin) {
66 l = begin + 1;
67 r = end;
68
69 while(l < r) {
70 /* This whole section is for checking lval<pivot, where
71 pivot=N*row[begin]+col[begin] and lval=N*row[l]+col[l]. */
72 if (row[l]<row[begin]) {
73 if (ABS(row[l]-row[begin])==1 && ABS(col[l]-col[begin])==N)
74 flag=0;
75 else
76 flag=1;
77 } else if (row[l]==row[begin]) {
78 if (col[l]<col[begin])
79 flag=1;
80 else
81 flag=0;
82 } else {
83 if (ABS(row[l]-row[begin])==1 && ABS(col[l]-col[begin])==N)
84 flag=1;
85 else
86 flag=0;
87 }
88
89 if (flag==1) {
90 l++;
91 } else {
92 r--;
93 swap(row, col, val, l, r);
94 }
95 }
96 l--;
97 swap(row, col, val, begin, l);
98 q_sort(row, col, val, begin, l);
99 q_sort(row, col, val, r, end);
100 }
101 }
102
103 SystemMatrix_ptr SystemMatrix::loadMM_toCSR(const char *filename)
104 {
105 index_t *col_ind = NULL;
106 index_t *row_ind = NULL;
107 index_t *row_ptr = NULL;
108 double *val = NULL;
109 SystemMatrix_ptr out;
110 int curr_row;
111 MM_typecode matrixCode;
112 esysUtils::JMPI mpi_info=esysUtils::makeInfo(MPI_COMM_WORLD);
113 Esys_resetError();
114 if (mpi_info->size > 1) {
115 Esys_setError(IO_ERROR, "SystemMatrix::loadMM_toCSR: supports single processor only");
116 return out;
117 }
118
119 // open the file
120 std::ifstream f(filename);
121 if (!f.good()) {
122 Esys_setError(IO_ERROR, "SystemMatrix::loadMM_toCSR: Cannot open file for reading.");
123 return out;
124 }
125
126 // process banner
127 if (mm_read_banner(f, &matrixCode) != 0) {
128 Esys_setError(IO_ERROR, "SystemMatrix::loadMM_toCSR: Error processing MM banner.");
129 f.close();
130 return out;
131 }
132 if ( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) ) {
133 Esys_setError(TYPE_ERROR, "SystemMatrix::loadMM_toCSR: found Matrix Market type is not supported.");
134 f.close();
135 return out;
136 }
137
138 // get matrix size
139 if (mm_read_mtx_crd_size(f, &M, &N, &nz) != 0) {
140 Esys_setError(IO_ERROR, "SystemMatrix::loadMM_toCSR: Could not read sparse matrix size.");
141 f.close();
142 return out;
143 }
144
145 // prepare storage
146 col_ind = new index_t[nz];
147 row_ind = new index_t[nz];
148 val = new double[nz];
149
150 row_ptr = new index_t[M+1];
151
152 // perform actual read of elements
153 for (int i=0; i<nz; i++) {
154 f >> row_ind[i] >> col_ind[i] >> val[i];
155 if (!f.good()) {
156 delete[] val;
157 delete[] row_ind;
158 delete[] col_ind;
159 delete[] row_ptr;
160 f.close();
161 return out;
162 }
163 row_ind[i]--;
164 col_ind[i]--;
165 }
166 f.close();
167
168 // sort the entries
169 q_sort(row_ind, col_ind, val, 0, nz);
170
171 // setup row_ptr
172 curr_row = 0;
173 for (int i=0; (i<nz && curr_row<M); curr_row++) {
174 while (row_ind[i] != curr_row) {
175 i++;
176 }
177 row_ptr[curr_row] = i;
178 }
179 row_ptr[M] = nz;
180
181 // create return value
182 index_t dist[2];
183 dist[0]=0;
184 dist[1]=M;
185 Distribution_ptr output_dist(new Distribution(mpi_info, dist, 1, 0));
186 dist[1]=N;
187 Distribution_ptr input_dist(new Distribution(mpi_info, dist, 1, 0));
188 Pattern_ptr mainPattern(new Pattern(MATRIX_FORMAT_DEFAULT, M, N, row_ptr, col_ind));
189 Pattern_ptr couplePattern(new Pattern(MATRIX_FORMAT_DEFAULT, M, N, NULL, NULL));
190 dist[0]=M;
191 SharedComponents_ptr send(new SharedComponents(
192 M, 0, NULL, NULL, dist, 1, 0, mpi_info));
193 Connector_ptr connector(new Connector(send, send));
194 SystemMatrixPattern_ptr pattern(new SystemMatrixPattern(
195 MATRIX_FORMAT_DEFAULT, output_dist, input_dist, mainPattern,
196 couplePattern, couplePattern, connector, connector));
197 out.reset(new SystemMatrix(MATRIX_FORMAT_DEFAULT, pattern, 1, 1, true,
198 escript::FunctionSpace(), escript::FunctionSpace()));
199
200 // copy values
201 #pragma omp parallel for
202 for (int i=0; i<nz; i++)
203 out->mainBlock->val[i] = val[i];
204
205 delete[] val;
206 delete[] row_ind;
207 return out;
208 }
209
210 SystemMatrix_ptr SystemMatrix::loadMM_toCSC(const char* filename)
211 {
212 Pattern_ptr mainPattern, couplePattern;
213 SystemMatrixPattern_ptr pattern;
214 SystemMatrix_ptr out;
215 Connector_ptr connector;
216 index_t *col_ind = NULL;
217 index_t *row_ind = NULL;
218 index_t *col_ptr = NULL;
219 double *val = NULL;
220 int curr_col=0;
221 MM_typecode matrixCode;
222 esysUtils::JMPI mpi_info=esysUtils::makeInfo( MPI_COMM_WORLD);
223 if (mpi_info->size > 1) {
224 Esys_setError(IO_ERROR, "SystemMatrix::loadMM_toCSC: supports single processor only");
225 return out;
226 }
227
228 Esys_resetError();
229
230 // open the file
231 std::ifstream f(filename);
232 if (!f.good()) {
233 Esys_setError(IO_ERROR, "SystemMatrix::loadMM_toCSC: File could not be opened for reading.");
234 return out;
235 }
236
237 // process banner
238 if (mm_read_banner(f, &matrixCode) != 0) {
239 Esys_setError(IO_ERROR,"SystemMatrix::loadMM_toCSC: Error processing MM banner.");
240 f.close();
241 return out;
242 }
243 if( !(mm_is_real(matrixCode) && mm_is_sparse(matrixCode) && mm_is_general(matrixCode)) ) {
244 Esys_setError(TYPE_ERROR, "SystemMatrix::loadMM_toCSC: found Matrix Market type is not supported.");
245 f.close();
246 return out;
247 }
248
249 // get matrix size
250 if (mm_read_mtx_crd_size(f, &M, &N, &nz) != 0) {
251 Esys_setError(TYPE_ERROR, "SystemMatrix::loadMM_toCSC: found Matrix Market type is not supported.");
252 f.close();
253 return out;
254 }
255
256 // prepare storage
257 col_ind = new index_t[nz];
258 row_ind = new index_t[nz];
259 val = new double [nz];
260 col_ptr = new index_t[N+1];
261
262 // perform actual read of elements
263 for (int i=0; i<nz; i++) {
264 f >> row_ind[i] >> col_ind[i] >> val[i];
265 if (!f.good()) {
266 delete[] val;
267 delete[] row_ind;
268 delete[] col_ind;
269 delete[] col_ptr;
270 f.close();
271 return out;
272 }
273 row_ind[i]--;
274 col_ind[i]--;
275 }
276 f.close();
277
278 // sort the entries
279 q_sort( col_ind, row_ind, val, 0, nz );
280
281 // setup row_ptr
282 for(int i=0; (i<nz && curr_col<N); curr_col++) {
283 while( col_ind[i] != curr_col )
284 i++;
285 col_ptr[curr_col] = i;
286 }
287 col_ptr[N] = nz;
288
289 index_t dist[2];
290 dist[0]=0;
291 dist[1]=N;
292 Distribution_ptr output_dist(new Distribution(mpi_info, dist,1,0));
293 dist[1]=M;
294 Distribution_ptr input_dist(new Distribution(mpi_info, dist,1,0));
295 mainPattern.reset(new Pattern(MATRIX_FORMAT_DEFAULT,N,M,col_ptr,col_ind));
296 couplePattern.reset(new Pattern(MATRIX_FORMAT_DEFAULT,N,M,NULL,NULL));
297 SharedComponents_ptr send(new SharedComponents(
298 N, 0, NULL, NULL, NULL, 1, 0, mpi_info));
299 connector.reset(new Connector(send,send));
300 pattern.reset(new SystemMatrixPattern(MATRIX_FORMAT_DEFAULT,
301 output_dist, input_dist, mainPattern, couplePattern,
302 couplePattern, connector, connector));
303 out.reset(new SystemMatrix(MATRIX_FORMAT_CSC, pattern, 1, 1, true,
304 escript::FunctionSpace(),
305 escript::FunctionSpace()));
306
307 // copy values
308 #pragma omp parallel for
309 for (int i=0; i<nz; i++)
310 out->mainBlock->val[i] = val[i];
311
312 delete[] val;
313 delete[] row_ind;
314 return out;
315 }
316
317 void RHS_loadMM_toCSR(const char *filename, double *b, dim_t size)
318 {
319 MM_typecode matrixCode;
320 Esys_resetError();
321 // open the file
322 std::ifstream f(filename);
323 if (!f.good()) {
324 Esys_setError(IO_ERROR, "RHS_loadMM_toCSR: Cannot open file for reading.");
325 }
326
327 // process banner
328 if (mm_read_banner(f, &matrixCode) != 0) {
329 Esys_setError(IO_ERROR, "RHS_loadMM_toCSR: Error processing MM banner.");
330 }
331 if( !(mm_is_real(matrixCode) && mm_is_general(matrixCode) && mm_is_array(matrixCode)) ) {
332 Esys_setError(TYPE_ERROR,"RHS_loadMM_toCSR: found Matrix Market type is not supported.");
333 }
334
335 // get matrix size
336 if (mm_read_mtx_array_size(f, &M, &N) != 0) {
337 Esys_setError(IO_ERROR, "RHS_loadMM_toCSR: Could not read sparse matrix size.");
338 }
339
340 if (M != size) {
341 Esys_setError(IO_ERROR, "RHS_loadMM_toCSR: Actual and provided sizes do not match.");
342 }
343
344 if (Esys_noError()) {
345 nz=M;
346 // perform actual read of elements
347 for (int i=0; i<nz; i++) {
348 f >> b[i];
349 if (!f.good()) {
350 f.close();
351 Esys_setError(IO_ERROR, "RHS_loadMM_toCSR: Could not read some of the values.");
352 }
353 }
354 }
355 f.close();
356 }
357
358 } // namespace paso
359

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/paso/src/SystemMatrix_loadMM.cpp:5567-5588 /branches/amg_from_3530/paso/src/SystemMatrix_loadMM.cpp:3531-3826 /branches/lapack2681/paso/src/SystemMatrix_loadMM.cpp:2682-2741 /branches/pasowrap/paso/src/SystemMatrix_loadMM.cpp:3661-3674 /branches/py3_attempt2/paso/src/SystemMatrix_loadMM.cpp:3871-3891 /branches/restext/paso/src/SystemMatrix_loadMM.cpp:2610-2624 /branches/ripleygmg_from_3668/paso/src/SystemMatrix_loadMM.cpp:3669-3791 /branches/stage3.0/paso/src/SystemMatrix_loadMM.cpp:2569-2590 /branches/symbolic_from_3470/paso/src/SystemMatrix_loadMM.cpp:3471-3974 /branches/symbolic_from_3470/ripley/test/python/paso/src/SystemMatrix_loadMM.cpp:3517-3974 /release/3.0/paso/src/SystemMatrix_loadMM.cpp:2591-2601 /release/4.0/paso/src/SystemMatrix_loadMM.cpp:5380-5406 /trunk/paso/src/SystemMatrix_loadMM.cpp:4257-4344,5898-5932 /trunk/ripley/test/python/paso/src/SystemMatrix_loadMM.cpp:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26