1 /*
2
3 Copyright (C) 2008-2021 Michele Martone
4
5 This file is part of librsb.
6
7 librsb is free software; you can redistribute it and/or modify it
8 under the terms of the GNU Lesser General Public License as published
9 by the Free Software Foundation; either version 3 of the License, or
10 (at your option) any later version.
11
12 librsb is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
16
17 You should have received a copy of the GNU Lesser General Public
18 License along with librsb; see the file COPYING.
19 If not, see <http://www.gnu.org/licenses/>.
20
21 */
22 /* @cond INNERDOC */
23 /*!
24 * @file
25 * @author Michele Martone
26 * @brief
27 * This source file contains matrix I/O functions.
28 * */
29
30 // FIXME: this code is messy and unclean
31
32 #include "rsb_internals.h"
33 #ifdef RSB_HAVE_SYS_STAT_H
34 #include <sys/stat.h>
35 #endif /* RSB_HAVE_SYS_STAT_H */
36 #if RSB_WANT_ZLIB_SUPPORT
37 #include <zlib.h>
38 #endif /* RSB_WANT_ZLIB_SUPPORT */
39 #include <stdio.h>
40 #define RSB_MMIOH_CL 4
41 #define RSB_20120321_IOBUFFERING 0
42 #if RSB_WANT_OMPIO_SUPPORT
43 #include "rsb_ompio.h"
44 #endif /* RSB_WANT_OMPIO_SUPPORT */
45 #define RSB_PATCH_GILLES_20130906 0
46 #if RSB_PATCH_GILLES_20130906
47 //#define _GNU_SOURCE
48 #include <asm/fcntl.h>
49 #include <unistd.h> /* O_DIRECT */
50 #include <sys/types.h> /* See NOTES */
51 #include <sys/stat.h>
52 #endif
53
54 #define rsb_util_sort_column_major(VA,IA,JA,nnz,nr,nc,typecode,flags) rsb__util_sort_row_major_inner(VA,JA,IA,nnz,nc,nr,typecode,flags)
55 #if 1
56 #define RSB_IO_VERBOSE_MSG(NZI,NNZ)
57 #else
58 #define RSB_IO_VERBOSE_GRNLRT RSB_MILLION_I
59 // #define RSB_IO_VERBOSE_GRNLRT 100000
60 #define RSB_IO_VERBOSE_MSG(NZI,NNZ) if((NZI)%RSB_IO_VERBOSE_GRNLRT==0) RSB_STDERR("%s%dM/%dM\n",RSB_CLEARTERM_STRING,(NZI)/RSB_IO_VERBOSE_GRNLRT,(NNZ)/RSB_IO_VERBOSE_GRNLRT )
61 #endif
62
63 #define RSB_FILE_ALLOW_LOAD_EMPTY_PATTERN 1 /* 20140324 */
64
65 #ifdef RSB_WITH_MM
rsb_util_mm_load_coo_matrix(const char * filename,struct rsb_coo_matrix_t * cmp)66 rsb_err_t rsb_util_mm_load_coo_matrix(const char *filename, struct rsb_coo_matrix_t * cmp)
67 {
68 /**
69 * \ingroup gr_internals
70 *
71 * Loads in a matrix in unsorted COO format. (new)
72 *
73 * FIXME : UNTESTED
74 *
75 * \note used by experiment.c files
76 */
77 rsb_err_t errval = RSB_ERR_NO_ERROR;
78 struct rsb_coo_matrix_t cm;
79 rsb_flags_t flags = RSB_FLAG_NOFLAGS;
80
81 if(!cmp || !filename)
82 {
83 errval = RSB_ERR_BADARGS;
84 RSB_PERR_GOTO(err,RSB_ERRM_ES);
85 }
86
87 RSB_BZERO_P(&cm);
88 cm.typecode = cmp->typecode; // should be specified
89 errval = rsb__util_mm_load_matrix_f(filename, &cm.IA, &cm.JA,&cm.VA , &cm.nr, &cm.nc, &cm.nnz , cm.typecode, flags, NULL, NULL);
90 if(RSB_SOME_ERROR(errval))
91 goto err;
92 rsb__memcpy(cmp,&cm,sizeof(cm));
93 err:
94 rsb__do_perror(NULL,errval);
95 RSB_DO_ERR_RETURN(errval)
96 }
97
rsb__util_mm_info_matrix_f(const char * fn,rsb_coo_idx_t * m,rsb_coo_idx_t * k,rsb_nnz_idx_t * nnz,rsb_type_t * typecode,rsb_bool_t * is_symmetric,rsb_bool_t * is_hermitian,rsb_bool_t * is_pattern,rsb_bool_t * is_lower,rsb_bool_t * is_upper,rsb_bool_t * is_vector)98 rsb_err_t rsb__util_mm_info_matrix_f(const char *fn, rsb_coo_idx_t *m, rsb_coo_idx_t *k , rsb_nnz_idx_t *nnz, rsb_type_t *typecode, rsb_bool_t * is_symmetric, rsb_bool_t * is_hermitian, rsb_bool_t * is_pattern, rsb_bool_t * is_lower, rsb_bool_t * is_upper , rsb_bool_t * is_vector)
99 {
100 /*!
101 * \ingroup internals
102 * FIXME : does not return cleanly in case of errors
103 * */
104
105 FILE * fd = NULL;
106 int innz = 0; /* FIXME */
107 int _m = 0,_k = 0;
108 char matcode[RSB_MMIOH_CL]; // !?
109 rsb_bool_t is_vector_ = RSB_BOOL_FALSE;
110
111 if(nnz)*nnz = RSB_MARKER_NNZ_VALUE ;
112 if(m)*m = RSB_MARKER_COO_VALUE;
113 if(k)*k = RSB_MARKER_COO_VALUE;
114 /* TODO: needs to define some RSB_NUMERICAL_COMPLEX_TYPE_DEFAULT macro and use it here, in case of a complex matrix */
115 if(typecode && RSB_MATRIX_UNSUPPORTED_TYPE(*typecode))*typecode = RSB_NUMERICAL_TYPE_DEFAULT;
116
117 if(!fn)
118 {
119 return RSB_ERR_BADARGS;
120 }
121
122 if(fd==NULL)
123 if ((fd = RSB_FOPEN(fn, "r")) == NULL)
124 {
125 RSB_STDERR("Failed opening file: %s\n",fn);
126 return RSB_ERR_GENERIC_ERROR;
127 }
128
129 if (rsb__mm_read_banner(fd,NULL,&(matcode)) != 0)
130 {
131 RSB_STDERR("Could not process Matrix Market banner.\n");
132 RSB_FCLOSE(fd);
133 return RSB_ERR_GENERIC_ERROR;
134 }
135
136 /* This is how one can screen matrix types if their application */
137 /* only supports a subset of the Matrix Market data types. */
138
139 is_vector_ = (rsb_mm_is_sparse(matcode))?RSB_BOOL_FALSE:RSB_BOOL_TRUE;
140 if ( !rsb_mm_is_matrix(matcode) /*|| !rsb_mm_is_sparse(matcode)*/ )
141 // if (!rsb_mm_is_real(matcode) || !rsb_mm_is_matrix(matcode) || !rsb_mm_is_sparse(matcode) )
142 {
143 RSB_STDERR("%s","Sorry, this application does not support ");
144 RSB_STDERR("Matrix Market type: [%s]\n", rsb__mm_typecode_to_str(matcode));
145 RSB_FCLOSE(fd);
146 return RSB_ERR_UNSUPPORTED_TYPE;
147 }
148
149 /* find out size of sparse matrix .... */
150
151
152 if( ((is_vector_) && (rsb__mm_read_mtx_array_size(fd,NULL,&_m,&_k) !=0)) || ((!is_vector_) && (rsb__mm_read_mtx_crd_size(fd,NULL,&_m,&_k,&innz)) !=0) )
153 {
154 RSB_FCLOSE(fd);
155 return RSB_ERR_GENERIC_ERROR;
156 }
157 if(m)*m = (rsb_coo_idx_t)_m;
158 if(k)*k = (rsb_coo_idx_t)_k;
159
160 if(is_vector)
161 {
162 *is_vector = is_vector_;
163 }
164
165 if(is_pattern)
166 {
167 if(rsb_mm_is_pattern(matcode))
168 *is_pattern = RSB_BOOL_TRUE;
169 else
170 *is_pattern = RSB_BOOL_FALSE;
171 }
172
173 if(is_symmetric)
174 {
175 if(rsb_mm_is_symmetric(matcode))
176 *is_symmetric = RSB_BOOL_TRUE;
177 else
178 *is_symmetric = RSB_BOOL_FALSE;
179 }
180
181 if(is_hermitian)
182 {
183 if(rsb_mm_is_hermitian(matcode))
184 *is_hermitian = RSB_BOOL_TRUE;
185 else
186 *is_hermitian = RSB_BOOL_FALSE;
187 }
188
189 if(is_lower)
190 {
191 // TODO
192 }
193
194 if(is_upper)
195 // TODO
196
197 if(m && k)
198 if (((int)*m != _m)||((int)*k != _k))
199 {
200 /* overflow */
201 RSB_IO_ERROR("overflow error while reading matrix dimensions.\n");
202 return RSB_ERR_INTERNAL_ERROR;
203 }
204
205 if(is_vector_)
206 innz = _m*_k; /* 20120904 */
207 if(nnz)
208 {
209 *nnz = innz;
210 #if RSB_FILE_ALLOW_LOAD_EMPTY_PATTERN
211 if(*nnz<0)
212 #else
213 if(*nnz<1)
214 #endif
215 return RSB_ERR_GENERIC_ERROR;
216 }
217
218 RSB_FCLOSE(fd);
219 return RSB_ERR_NO_ERROR;
220 }
221
rsb_zfscanf(FILE * fd,const char * fs,rsb_coo_idx_t * IV,rsb_coo_idx_t * JV,void * VAR,void * VAI,void * ngzfd)222 static int rsb_zfscanf(FILE * fd,const char * fs,rsb_coo_idx_t *IV, rsb_coo_idx_t *JV, void * VAR, void * VAI, void * ngzfd)
223 {
224 /**
225 * \ingroup internals
226 * FIXME: error handling is missing
227 * */
228 #if RSB_WANT_ZLIB_SUPPORT
229 if(ngzfd)
230 {
231 if((!IV) && (!JV))
232 {
233 if(VAI)
234 return fscanf(ngzfd,fs,VAR,VAI);
235 else
236 return fscanf(ngzfd,fs,VAR);
237 }
238 else
239 {
240 if(VAI)
241 return fscanf(ngzfd,fs,IV,JV,VAR,VAI);
242 if(VAR)
243 return fscanf(ngzfd,fs,IV,JV,VAR);
244 else
245 return fscanf(ngzfd,fs,IV,JV);
246 }
247 }
248 else
249 #endif /* RSB_WANT_ZLIB_SUPPORT */
250 return rsb_fscanf(fd,fs,IV,JV,VAR,VAI);
251 }
252
rsb_fscanf(FILE * fd,const char * fs,rsb_coo_idx_t * IV,rsb_coo_idx_t * JV,void * VAR,void * VAI)253 int rsb_fscanf(FILE * fd,const char * fs,rsb_coo_idx_t *IV, rsb_coo_idx_t *JV, void * VAR, void * VAI)
254 {
255 /**
256 * \ingroup internals
257 * FIXME: error handling is missing
258 * */
259 #if RSB_WANT_ZLIB_SUPPORT
260 char line[MM_MAX_LINE_LENGTH];
261 gzgets(fd,line,MM_MAX_LINE_LENGTH);
262
263 if((!IV) && (!JV))
264 {
265 if(VAI)
266 return sscanf(line,fs,VAR,VAI);
267 if(VAR)
268 return sscanf(line,fs,VAR);
269 else
270 return 0;
271 }
272 else
273 {
274 if(VAI)
275 return sscanf(line,fs,IV,JV,VAR,VAI);
276 if(VAR)
277 return sscanf(line,fs,IV,JV,VAR);
278 else
279 return sscanf(line,fs,IV,JV);
280 }
281 #else /* RSB_WANT_ZLIB_SUPPORT */
282 if((!IV) && (!JV))
283 {
284 if(VAI)
285 return fscanf(fd,fs,VAR,VAI);
286 if(VAR)
287 return fscanf(fd,fs,VAR);
288 else
289 return 0;
290 }
291 else
292 {
293 if(VAI)
294 return fscanf(fd,fs,IV,JV,VAR,VAI);
295 if(VAR)
296 return fscanf(fd,fs,IV,JV,VAR);
297 else
298 return fscanf(fd,fs,IV,JV);
299 }
300 #endif /* RSB_WANT_ZLIB_SUPPORT */
301 }
302
rsb_fgets(char * RSB_RESTRICT buf,int len,FILE * RSB_RESTRICT fd)303 char * rsb_fgets(char* RSB_RESTRICT buf, int len, FILE * RSB_RESTRICT fd)
304 {
305 /**
306 * \ingroup internals
307 * */
308 #if RSB_WANT_ZLIB_SUPPORT
309 return gzgets(fd,buf,len);
310 #else /* RSB_WANT_ZLIB_SUPPORT */
311 return fgets(buf,len,fd);
312 #endif /* RSB_WANT_ZLIB_SUPPORT */
313 }
314
rsb__util_mm_load_vector_f(const char * fn,void ** VA,rsb_nnz_idx_t * nnz,rsb_type_t typecode)315 rsb_err_t rsb__util_mm_load_vector_f(const char *fn, void **VA, rsb_nnz_idx_t *nnz, rsb_type_t typecode)
316 {
317 /* FIXME: in perpective, need stride, C/Fortran order, etc ... */
318 rsb_err_t errval = RSB_ERR_NO_ERROR;
319 errval = rsb__util_mm_load_matrix_f(fn,NULL,NULL,VA,NULL,NULL,nnz,typecode,RSB_FLAG_NOFLAGS,NULL,NULL);
320 return errval;
321 }
322
rsb__util_mm_load_matrix_f(const char * fn,rsb_coo_idx_t ** IA,rsb_coo_idx_t ** JA,void ** VA,rsb_coo_idx_t * m,rsb_coo_idx_t * k,rsb_nnz_idx_t * nnz,rsb_type_t typecode,rsb_flags_t flags,rsb_bool_t * is_lowerp,rsb_bool_t * is_upperp)323 rsb_err_t rsb__util_mm_load_matrix_f(const char *fn, rsb_coo_idx_t ** IA, rsb_coo_idx_t ** JA, void **VA, rsb_coo_idx_t *m, rsb_coo_idx_t *k , rsb_nnz_idx_t *nnz, rsb_type_t typecode, rsb_flags_t flags, rsb_bool_t *is_lowerp, rsb_bool_t *is_upperp)
324 {
325 /**
326 * \ingroup internals
327 * This function reads in a single Matrix Market format matrix drom a specified filename.
328 * The matrix will be loaded in coordinate storage format, and IA,JA, VA arrays will
329 * be allocated for this purpose, here.
330 *
331 * \param fn should contain a valid matrix file name
332 * \param IA should point to a pointer which will be allocated here to contain the elements row values
333 * \param JA should point to a pointer which will be allocated here to contain the elements column values
334 * \param VA should point to a pointer which will be allocated here to contain the elements column values
335 * \param m should point to the matrix rows count variable, which will be set in this function
336 * \param k should point to the matrix columns count variable, which will be set in this function
337 * \param nnz should point to the matrix nonzero count variable, which will be set in this function (shall be initialized to zero in advance!)
338 * \param type should specify a valid numerical type to convert the read data in. see rsb.h for this.
339 * \return RSB_ERR_NO_ERROR on correct operation, an error code (see \ref errors_section) otherwise.
340 *
341 * Notes:
342 *
343 * There is no guarantee that the loaded matrix will be free from duplicate values.
344 * The specified numerical type should be enabled at library generaton time.
345 *
346 * FIXME : lots of ths function's code should be generated from macros.
347 * FIXME : error handling is awful
348 * FIXME : otype stands for 'original' or 'output' ?
349 * TODO : some option to specify a pattern-only load
350 * FIXME : in a future version, should not allocate if pointers not NULL
351 * TODO: may print (at least for rsbench) out how many discarded duplicates, how many zeroes, etc etc
352 * TODO: may support different styles for tolerating / detecting e.g. reading a vector file instead a matrix one.
353 * */
354 rsb_err_t errval = RSB_ERR_NO_ERROR;
355 size_t re = 0;/* read elements */
356 FILE *fd = NULL;
357 //#if RSB_WANT_ZLIB_SUPPORT
358 FILE *ngzfd = NULL;
359 //#endif
360 rsb_nnz_idx_t i = 0;
361 rsb_nnz_idx_t annz = 0;/* allocated nnz */
362 rsb_bool_t is_symmetric = RSB_BOOL_FALSE,is_lower = RSB_BOOL_FALSE,is_upper = RSB_BOOL_FALSE,is_pattern = RSB_BOOL_FALSE, is_hermitian = RSB_BOOL_FALSE, is_vector = RSB_BOOL_FALSE;/* FIXME : no expansion support for hermitian */
363 rsb_bool_t aja = RSB_BOOL_FALSE, ava = RSB_BOOL_FALSE, aia = RSB_BOOL_FALSE;
364 rsb_flags_t otype = typecode;/* original type */
365 rsb_time_t frt = 0;
366 char matcode[RSB_MMIOH_CL]; // !?
367 rsb_bool_t is_gz = RSB_BOOL_FALSE;
368
369 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
370 double **dval = NULL;
371 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
372 #ifdef RSB_NUMERICAL_TYPE_FLOAT
373 float **fval = NULL;
374 #endif /* RSB_NUMERICAL_TYPE_FLOAT */
375 #ifdef RSB_NUMERICAL_TYPE_CHAR
376 char **cval = NULL;
377 #endif /* RSB_NUMERICAL_TYPE_CHAR */
378 #ifdef RSB_NUMERICAL_TYPE_INT
379 int **ival = NULL;
380 #endif /* RSB_NUMERICAL_TYPE_INT */
381 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
382 float complex **zval = NULL;
383 #endif /* RSB_NUMERICAL_TYPE_FLOAT_COMPLEX */
384 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
385 double complex **Zval = NULL;
386 #endif /* RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX */
387 int innz = 0; /* FIXME */
388 int _m = 0,_k = 0;
389 #if RSB_20120321_IOBUFFERING
390 char*iobuf = NULL;
391 size_t iobbs = 16*1024*1024;
392 #endif /* RSB_20120321_IOBUFFERING */
393
394 frt = - rsb_time();
395
396 if ( RSB_MATRIX_UNSUPPORTED_TYPE(typecode) )return RSB_ERR_UNSUPPORTED_TYPE ;
397
398 //if(!VA || !JA || !IA) return RSB_ERR_BADARGS;
399 if( (!(VA && JA && IA)) && (!(VA && (!JA) && (!IA))) ) return RSB_ERR_BADARGS;
400 if(IA)aia = RSB_BOOL_IS_POINTER_NON_NULL(*IA);
401 if(JA)aja = RSB_BOOL_IS_POINTER_NON_NULL(*JA);
402 if(VA)ava = RSB_BOOL_IS_POINTER_NON_NULL(*VA);
403 //if(!nnz || !k || !m) return RSB_ERR_BADARGS;
404 if((!(nnz && k && m)) && (!(nnz && (!k) && (!m)))) return RSB_ERR_BADARGS;
405 if(*nnz)
406 annz = *nnz;// if user set
407
408 errval = rsb__util_mm_info_matrix_f(fn,m,k,nnz,&typecode,&is_symmetric,&is_hermitian,&is_pattern,&is_lower,&is_upper,&is_vector);
409 if(RSB_SOME_ERROR(errval))
410 goto prerr;
411 if(annz==0)
412 annz = *nnz;
413 if(annz<*nnz)
414 {
415 RSB_ERROR("user-set array size (%d) does not fit actual input (%d)\n",(int)*nnz,(int)annz);
416 errval = RSB_ERR_BADARGS;
417 goto prerr;
418 }
419
420 if(is_pattern)
421 #ifdef RSB_NUMERICAL_TYPE_PATTERN
422 typecode = RSB_NUMERICAL_TYPE_PATTERN;
423 #else /* RSB_NUMERICAL_TYPE_PATTERN */
424 return RSB_ERR_UNSUPPORTED_FEATURE;
425 #endif /* RSB_NUMERICAL_TYPE_PATTERN */
426
427 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
428 if (otype == RSB_NUMERICAL_TYPE_DOUBLE)dval = (double**)(VA);
429 else
430 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
431 #ifdef RSB_NUMERICAL_TYPE_FLOAT
432 if (otype == RSB_NUMERICAL_TYPE_FLOAT)fval = (float **)(VA);
433 else
434 #endif /* RSB_NUMERICAL_TYPE_FLOAT */
435 #ifdef RSB_NUMERICAL_TYPE_INT
436 if (otype == RSB_NUMERICAL_TYPE_INT)ival = (int **)(VA);
437 else
438 #endif /* RSB_NUMERICAL_TYPE_INT */
439 #ifdef RSB_NUMERICAL_TYPE_CHAR
440 if (otype == RSB_NUMERICAL_TYPE_CHAR)cval = (char **)(VA);
441 else
442 #endif /* RSB_NUMERICAL_TYPE_CHAR */
443 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
444 if (otype == RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX)Zval = (double complex **)(VA);
445 else
446 #endif /* RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX */
447 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
448 if (otype == RSB_NUMERICAL_TYPE_FLOAT_COMPLEX)zval = (float complex **)(VA);
449 else
450 #endif /* RSB_NUMERICAL_TYPE_FLOAT_COMPLEX */
451 #ifdef RSB_NUMERICAL_TYPE_PATTERN
452 if (otype == RSB_NUMERICAL_TYPE_PATTERN){/* nothing to do */}
453 else
454 #endif /* RSB_NUMERICAL_TYPE_PATTERN */
455 return RSB_ERR_UNSUPPORTED_TYPE ;
456
457 if(IA)
458 *IA = NULL;
459 if(JA)
460 *JA = NULL;
461
462 {
463 // THIS IS A CODE DUPLICATION ..
464 is_gz = RSB_BOOL_FALSE;
465 #if RSB_WANT_ZLIB_SUPPORT
466 if ((fd = gzopen(fn,"r")) != NULL)
467 {
468 if(gzdirect(fd))
469 ;
470 else
471 is_gz = RSB_BOOL_TRUE;
472 gzclose(fd);
473 fd = NULL;
474 }
475 #endif /* RSB_WANT_ZLIB_SUPPORT */
476
477 #if RSB_WANT_ZLIB_SUPPORT
478 if(is_gz)
479 {
480 if((fd = gzopen(fn,"r")) == NULL)
481 {
482 /* TODO: the following code is not robust, shall fix it. */
483 #ifndef RSB_HAVE_DUP
484 #error Functions 'dup/fileno' is not present ? Reconfigure without Z library then!
485 int fnum = 0;/* */
486 ngzfd = NULL;
487 #else /* RSB_HAVE_DUP */
488 int fnum = dup( rsb__fileno(fd) );
489 ngzfd = gzdopen(fnum,"r");
490 #endif /* RSB_HAVE_DUP */
491 if(!ngzfd)
492 {
493 gzclose(fd);
494 RSB_ERROR(RSB_ERRMSG_FILEOPENPGZ"\n");
495 return RSB_ERR_GENERIC_ERROR;
496 }
497 else
498 ;
499 }
500 }
501 else
502 #endif /* RSB_WANT_ZLIB_SUPPORT */
503 #if !RSB_PATCH_GILLES_20130906
504 if ((fd = fopen(fn,"r")) == NULL)
505 #else /* RSB_PATCH_GILLES_20130906 */
506 {
507 int _fd;
508 if ((_fd = open(fn,O_RDONLY|O_DIRECT)) < 0)
509 {
510 RSB_STDERR(RSB_ERRMSG_FILEOPENP" %s\n",fn);
511 return RSB_ERR_GENERIC_ERROR;
512 }
513 else if ((fd = fdopen(_fd,"r")) == NULL)
514 #endif /* RSB_PATCH_GILLES_20130906 */
515 {
516 RSB_STDERR(RSB_ERRMSG_FILEOPENP" %s\n",fn);
517 return RSB_ERR_GENERIC_ERROR;
518 }
519 else
520 ngzfd = fd;
521 #if RSB_20120321_IOBUFFERING
522 //iobbs = BUFSIZ;
523 if(iobbs>0)
524 if(((iobuf = rsb__malloc(iobbs))==NULL)
525 //|| (0!=setvbuf(ngzfd,NULL,_IOLBF,0))// line buffering: slow
526 //|| (0!=setvbuf(ngzfd,NULL,_IONBF,0))// no buffering: super-slow
527 || (0!=setvbuf(ngzfd,iobuf,_IOFBF,iobbs)) // seem to not work
528 )
529 {
530 RSB_STDERR("problems setting up a buffer for file ""%s\n",fn);
531 }
532 //setbuf(ngzfd,iobbs);
533 //setbuffer(ngzfd,iobuf,iobbs);
534 #endif /* RSB_20120321_IOBUFFERING */
535 #if RSB_PATCH_GILLES_20130906
536 }
537 #endif /* RSB_PATCH_GILLES_20130906 */
538
539 if (rsb__mm_read_banner(fd,ngzfd,&(matcode)) != 0)
540 {
541 RSB_STDERR(RSB_ERRMSG_TMXMKTBANNER".\n");
542 RSB_FCLOSE(fd);
543 return RSB_ERR_GENERIC_ERROR;
544 }
545
546 if( ((is_vector) && (rsb__mm_read_mtx_array_size(fd,ngzfd,&_m,&_k) !=0)) || ((!is_vector) && (rsb__mm_read_mtx_crd_size(fd,ngzfd,&_m,&_k,&innz)) !=0) )
547 {
548 RSB_FCLOSE(fd);
549 return RSB_ERR_GENERIC_ERROR;
550 }
551 if(m)
552 *m = (rsb_coo_idx_t)_m;
553 if(k)
554 *k = (rsb_coo_idx_t)_k;
555 }
556
557 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
558 if (otype == RSB_NUMERICAL_TYPE_DOUBLE)*dval = *VA?*VA: rsb__calloc(sizeof(double)*annz);
559 else
560 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
561 #ifdef RSB_NUMERICAL_TYPE_FLOAT
562 if (otype == RSB_NUMERICAL_TYPE_FLOAT) *fval =*VA?*VA: rsb__calloc(sizeof(float) *annz);
563 else
564 #endif /* RSB_NUMERICAL_TYPE_FLOAT */
565 #ifdef RSB_NUMERICAL_TYPE_INT
566 if (otype == RSB_NUMERICAL_TYPE_INT) *ival =*VA?*VA: rsb__calloc(sizeof(int) *annz);
567 else
568 #endif /* RSB_NUMERICAL_TYPE_INT */
569 #ifdef RSB_NUMERICAL_TYPE_CHAR
570 if (otype == RSB_NUMERICAL_TYPE_CHAR) *cval =*VA?*VA: rsb__calloc(sizeof(char) *annz);
571 else
572 #endif /* RSB_NUMERICAL_TYPE_CHAR */
573 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
574 if (otype == RSB_NUMERICAL_TYPE_FLOAT_COMPLEX) *zval =*VA?*VA: rsb__calloc(sizeof(float complex) *annz);
575 else
576 #endif /* RSB_NUMERICAL_TYPE_FLOAT_COMPLEX */
577 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
578 if (otype == RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX) *Zval =*VA?*VA: rsb__calloc(sizeof(double complex) *annz);
579 else
580 #endif /* RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX */
581 #ifdef RSB_NUMERICAL_TYPE_PATTERN
582 if (otype == RSB_NUMERICAL_TYPE_PATTERN){/* no zval allocation (TODO : document this) */}
583 else
584 #endif /* RSB_NUMERICAL_TYPE_PATTERN */
585 {errval = RSB_ERR_UNSUPPORTED_TYPE;RSB_PERR_GOTO(err,RSB_ERRM_ES);}
586
587 if(IA && JA)
588 {
589 if( flags & RSB_FLAG_EXPERIMENTAL_IN_PLACE_CSR )
590 {
591 RSB_WARN("RSB_FLAG_EXPERIMENTAL_IN_PLACE_CSR will break with non matching index types\n");
592 /* Allocating slightly oversized arrays.
593 FIXME : potential type/size mismatches here !
594 */
595 if(!*IA)
596 *IA = rsb__calloc(sizeof(rsb_nnz_idx_t)*((annz>*m?annz:*m)+1));
597 if(!*JA)
598 *JA = rsb__calloc(sizeof(rsb_nnz_idx_t)*(annz+1));
599 }
600 else
601 {
602 if(!*IA)
603 *IA = rsb__calloc(sizeof(rsb_coo_idx_t)*annz);
604 if(!*JA)
605 *JA = rsb__calloc(sizeof(rsb_coo_idx_t)*annz);
606 }
607 if( !*IA || !*JA)
608 goto err;
609 }
610
611 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
612 if (otype == RSB_NUMERICAL_TYPE_DOUBLE) if(!dval) {RSB_PERR_GOTO(err,RSB_ERRM_ES)}
613 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
614 #ifdef RSB_NUMERICAL_TYPE_FLOAT
615 if (otype == RSB_NUMERICAL_TYPE_FLOAT ) if(!fval) {RSB_PERR_GOTO(err,RSB_ERRM_ES)}
616 #endif /* RSB_NUMERICAL_TYPE_FLOAT */
617 #ifdef RSB_NUMERICAL_TYPE_INT
618 if (otype == RSB_NUMERICAL_TYPE_INT ) if(!ival) {RSB_PERR_GOTO(err,RSB_ERRM_ES)}
619 #endif /* RSB_NUMERICAL_TYPE_INT */
620 #ifdef RSB_NUMERICAL_TYPE_CHAR
621 if (otype == RSB_NUMERICAL_TYPE_CHAR ) if(!cval) {RSB_PERR_GOTO(err,RSB_ERRM_ES)}
622 #endif /* RSB_NUMERICAL_TYPE_CHAR */
623 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
624 if (otype == RSB_NUMERICAL_TYPE_FLOAT_COMPLEX ) if(!zval) {RSB_PERR_GOTO(err,RSB_ERRM_ES)}
625 #endif /* RSB_NUMERICAL_TYPE_FLOAT_COMPLEX */
626 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
627 if (otype == RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX ) if(!Zval) {RSB_PERR_GOTO(err,RSB_ERRM_ES)}
628 #endif /* RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX */
629 #ifdef RSB_NUMERICAL_TYPE_PATTERN
630 /* :) */
631 #endif /* RSB_NUMERICAL_TYPE_PATTERN */
632
633 if( IA && JA)
634 goto full_scan;
635
636 #if RSB_WANT_OMPIO_SUPPORT
637 {errval = RSB_ERR_UNIMPLEMENTED_YET;RSB_PERR_GOTO(err,RSB_ERRM_ES);}
638 #endif /* RSB_WANT_OMPIO_SUPPORT */
639
640 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
641 if (typecode == RSB_NUMERICAL_TYPE_DOUBLE)
642 {
643 double iv;
644 if(rsb_mm_is_complex(matcode))
645 for (i=0; i<*nnz; i++)
646 {
647 re += (rsb_zfscanf(fd,"%lg %lg\n",NULL,NULL,*dval+i,&(iv),ngzfd)==2);
648 RSB_IO_VERBOSE_MSG(i,*nnz);
649 }
650 else
651 for (i=0; i<*nnz; i++)
652 {
653 re += (rsb_zfscanf(fd,"%lg\n",NULL,NULL,*dval+i,NULL,ngzfd)==1);
654 RSB_IO_VERBOSE_MSG(i,*nnz);
655 }
656 }
657 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
658
659 #ifdef RSB_NUMERICAL_TYPE_FLOAT
660 if (typecode == RSB_NUMERICAL_TYPE_FLOAT)
661 {
662 float iv;
663 if(rsb_mm_is_complex(matcode))
664 for (i=0; i<*nnz; i++)
665 {
666 re += (rsb_zfscanf(fd,"%g %g\n",NULL,NULL,*fval+i,&(iv),ngzfd)==2);
667 RSB_IO_VERBOSE_MSG(i,*nnz);
668 }
669 else
670 for (i=0; i<*nnz; i++)
671 {
672 re += (rsb_zfscanf(fd, "%g\n",NULL,NULL,*fval+i,NULL,ngzfd)==1);
673 RSB_IO_VERBOSE_MSG(i,*nnz);
674 }
675 }
676 #endif /* RSB_NUMERICAL_TYPE_FLOAT */
677
678 #ifdef RSB_NUMERICAL_TYPE_INT
679 if (typecode == RSB_NUMERICAL_TYPE_INT)
680 for (i=0; i<*nnz; i++)
681 {
682 double fv;
683 re += (rsb_zfscanf(fd,"%lg\n",NULL,NULL,&fv,NULL,ngzfd)==1);
684 (*ival)[i] = (int)fv;
685 RSB_IO_VERBOSE_MSG(i,*nnz);
686 }
687 #endif /* RSB_NUMERICAL_TYPE_INT */
688
689 #ifdef RSB_NUMERICAL_TYPE_CHAR
690 if (typecode == RSB_NUMERICAL_TYPE_CHAR)
691 for (i=0; i<*nnz; i++)
692 {
693 double fv;
694 re += (rsb_zfscanf(fd,"%g\n",NULL,NULL,&fv,NULL,ngzfd)==1);
695 (*cval)[i] = (char)fv;
696 RSB_IO_VERBOSE_MSG(i,*nnz);
697 }
698 #endif /* RSB_NUMERICAL_TYPE_CHAR */
699
700 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
701 if (typecode == RSB_NUMERICAL_TYPE_FLOAT_COMPLEX)
702 {
703 if(rsb_mm_is_complex(matcode))
704 for (i=0; i<*nnz; i++)
705 {
706 float rv,iv;
707 re += (rsb_zfscanf(fd,"%g %g\n",NULL,NULL,&(rv),&(iv),ngzfd)==2);
708 (*zval)[i] = (rv + I * iv);
709 RSB_IO_VERBOSE_MSG(i,*nnz);
710 }
711 else
712 for (i=0; i<*nnz; i++)
713 {
714 float rv;
715 re += (rsb_zfscanf(fd,"%g\n",NULL,NULL,&(rv),NULL,ngzfd)==1);
716 (*zval)[i] = (rv + I * 0);
717 RSB_IO_VERBOSE_MSG(i,*nnz);
718 }
719 }
720 #endif /* RSB_NUMERICAL_TYPE_FLOAT_COMPLEX */
721
722 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
723 if (typecode == RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX)
724 {
725 if(rsb_mm_is_complex(matcode))
726 for (i=0; i<*nnz; i++)
727 {
728 double rv,iv;
729 re += (rsb_zfscanf(fd,"%lg %lg\n",NULL,NULL,&(rv),&(iv),ngzfd)==2);
730 (*Zval)[i] = (rv + I * iv);
731 RSB_IO_VERBOSE_MSG(i,*nnz);
732 }
733 else
734 for (i=0; i<*nnz; i++)
735 {
736 double rv;
737 re += (rsb_zfscanf(fd,"%lg\n",NULL,NULL,&(rv),NULL,ngzfd)==1);
738 (*Zval)[i] = (rv + I * 0);
739 RSB_IO_VERBOSE_MSG(i,*nnz);
740 }
741 }
742 #endif /* RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX */
743 goto scan_done;
744 full_scan:
745 /* NOTE: when reading in doubles, ANSI C requires the use of the "l" */
746 /* specifier as in "%lg", "%lf", "%le", otherwise errors will occur */
747 /* (ANSI C X3.159-1989, Sec. 4.9.6.2, p. 136 lines 13-15) */
748
749 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
750 if (typecode == RSB_NUMERICAL_TYPE_DOUBLE)
751 {
752 if(rsb_mm_is_complex(matcode))
753 #if RSB_WANT_OMPIO_SUPPORT
754 {errval = RSB_ERR_UNIMPLEMENTED_YET;RSB_PERR_GOTO(err,RSB_ERRM_ES);}
755 #else /* RSB_WANT_OMPIO_SUPPORT */
756 for (i=0; i<*nnz; i++)
757 {
758 int iI,iJ;
759 double iv;
760 re += (rsb_zfscanf(fd,"%d %d %lg %lg\n",&iI,&iJ,*dval+i,&(iv),ngzfd)==4);
761 (*IA)[i] = (rsb_coo_idx_t)iI;
762 (*JA)[i] = (rsb_coo_idx_t)iJ;
763 (*IA)[i]--; /* adjust from 1-based to 0-based */
764 (*JA)[i]--;
765 RSB_IO_VERBOSE_MSG(i,*nnz);
766 }
767 #endif /* RSB_WANT_OMPIO_SUPPORT */
768 else
769 #if RSB_WANT_OMPIO_SUPPORT
770 rsb__ompio_DOUBLE(nnz,fd,ngzfd,dval,IA,JA,&re);
771 #else /* RSB_WANT_OMPIO_SUPPORT */
772 for (i=0; i<*nnz; i++)
773 {
774 int iI,iJ;
775 re += (rsb_zfscanf(fd,"%d %d %lg\n",&iI,&iJ,*dval+i,NULL,ngzfd)==3);
776 (*IA)[i] = (rsb_coo_idx_t)iI;
777 (*JA)[i] = (rsb_coo_idx_t)iJ;
778 (*IA)[i]--; /* adjust from 1-based to 0-based */
779 (*JA)[i]--;
780 RSB_IO_VERBOSE_MSG(i,*nnz);
781 }
782 #endif /* RSB_WANT_OMPIO_SUPPORT */
783 }
784 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
785
786 #ifdef RSB_NUMERICAL_TYPE_FLOAT
787 if (typecode == RSB_NUMERICAL_TYPE_FLOAT)
788 {
789 if(rsb_mm_is_complex(matcode))
790 #if RSB_WANT_OMPIO_SUPPORT
791 {errval = RSB_ERR_UNIMPLEMENTED_YET;RSB_PERR_GOTO(err,RSB_ERRM_ES);}
792 #else /* RSB_WANT_OMPIO_SUPPORT */
793 for (i=0; i<*nnz; i++)
794 {
795 int iI,iJ;
796 float iv;
797 re += (rsb_zfscanf(fd,"%d %d %g %g\n",&iI,&iJ,*fval+i,&(iv),ngzfd)==4);
798 (*IA)[i] = (rsb_coo_idx_t)iI;
799 (*JA)[i] = (rsb_coo_idx_t)iJ;
800 (*IA)[i]--; /* adjust from 1-based to 0-based */
801 (*JA)[i]--;
802 RSB_IO_VERBOSE_MSG(i,*nnz);
803 }
804 #endif /* RSB_WANT_OMPIO_SUPPORT */
805 else
806 #if RSB_WANT_OMPIO_SUPPORT
807 rsb_ompio_FLOAT(nnz,fd,ngzfd,dval,IA,JA,&re);
808 #else /* RSB_WANT_OMPIO_SUPPORT */
809 for (i=0; i<*nnz; i++)
810 {
811 int iI,iJ;
812 re += (rsb_zfscanf(fd, "%d %d %g\n",&iI,&iJ,*fval+i,NULL,ngzfd)==3);
813 (*IA)[i] = (rsb_coo_idx_t)iI;
814 (*JA)[i] = (rsb_coo_idx_t)iJ;
815 (*IA)[i]--; /* adjust from 1-based to 0-based */
816 (*JA)[i]--;
817 RSB_IO_VERBOSE_MSG(i,*nnz);
818 }
819 #endif /* RSB_WANT_OMPIO_SUPPORT */
820 }
821 #endif /* RSB_NUMERICAL_TYPE_FLOAT */
822
823 #ifdef RSB_NUMERICAL_TYPE_INT
824 if (typecode == RSB_NUMERICAL_TYPE_INT)
825 #if RSB_WANT_OMPIO_SUPPORT
826 rsb_ompio_INT(nnz,fd,ngzfd,dval,IA,JA,&re);
827 #else /* RSB_WANT_OMPIO_SUPPORT */
828 for (i=0; i<*nnz; i++)
829 {
830 int iI,iJ;
831 double fv;
832 re += (rsb_zfscanf(fd,"%d %d %lg\n",&iI,&iJ,&fv,NULL,ngzfd)==3);
833 (*IA)[i] = (rsb_coo_idx_t)iI;
834 (*JA)[i] = (rsb_coo_idx_t)iJ;
835 (*IA)[i]--; /* adjust from 1-based to 0-based */
836 (*JA)[i]--;
837 (*ival)[i] = (int)fv;
838 RSB_IO_VERBOSE_MSG(i,*nnz);
839 }
840 #endif /* RSB_WANT_OMPIO_SUPPORT */
841 #endif /* RSB_NUMERICAL_TYPE_INT */
842
843 #ifdef RSB_NUMERICAL_TYPE_CHAR
844 if (typecode == RSB_NUMERICAL_TYPE_CHAR)
845 #if RSB_WANT_OMPIO_SUPPORT
846 rsb_ompio_CHAR(nnz,fd,ngzfd,dval,IA,JA,&re);
847 #else /* RSB_WANT_OMPIO_SUPPORT */
848 for (i=0; i<*nnz; i++)
849 {
850 int iI,iJ;
851 double fv;
852 re += (rsb_zfscanf(fd,"%d %d %g\n",&iI,&iJ,&fv,NULL,ngzfd)==3);
853 (*IA)[i] = (rsb_coo_idx_t)iI;
854 (*JA)[i] = (rsb_coo_idx_t)iJ;
855 (*IA)[i]--; /* adjust from 1-based to 0-based */
856 (*JA)[i]--;
857 (*cval)[i] = (char)fv;
858 RSB_IO_VERBOSE_MSG(i,*nnz);
859 }
860 #endif /* RSB_WANT_OMPIO_SUPPORT */
861 #endif /* RSB_NUMERICAL_TYPE_INT */
862
863 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
864 if (typecode == RSB_NUMERICAL_TYPE_FLOAT_COMPLEX)
865 {
866 if(rsb_mm_is_complex(matcode))
867 #if RSB_WANT_OMPIO_SUPPORT
868 rsb_ompio_FLOAT_COMPLEX(nnz,fd,ngzfd,dval,IA,JA,&re);
869 #else /* RSB_WANT_OMPIO_SUPPORT */
870 for (i=0; i<*nnz; i++)
871 {
872 int iI,iJ;
873 float rv,iv;
874 re += (rsb_zfscanf(fd,"%d %d %g %g\n",&iI,&iJ,&(rv),&(iv),ngzfd)==4);
875 (*IA)[i] = (rsb_coo_idx_t)iI;
876 (*JA)[i] = (rsb_coo_idx_t)iJ;
877 (*IA)[i]--; /* adjust from 1-based to 0-based */
878 (*JA)[i]--;
879 (*zval)[i] = (rv + I * iv);
880 RSB_IO_VERBOSE_MSG(i,*nnz);
881 }
882 #endif /* RSB_WANT_OMPIO_SUPPORT */
883 else
884 #if RSB_WANT_OMPIO_SUPPORT
885 {errval = RSB_ERR_UNIMPLEMENTED_YET;RSB_PERR_GOTO(err,RSB_ERRM_ES);}
886 #else /* RSB_WANT_OMPIO_SUPPORT */
887 for (i=0; i<*nnz; i++)
888 {
889 int iI,iJ;
890 float rv;
891 re += (rsb_zfscanf(fd,"%d %d %g\n",&iI,&iJ,&(rv),NULL,ngzfd)==3);
892 (*IA)[i] = (rsb_coo_idx_t)iI;
893 (*JA)[i] = (rsb_coo_idx_t)iJ;
894 (*IA)[i]--; /* adjust from 1-based to 0-based */
895 (*JA)[i]--;
896 (*zval)[i] = (rv + I * 0);
897 RSB_IO_VERBOSE_MSG(i,*nnz);
898 }
899 #endif /* RSB_WANT_OMPIO_SUPPORT */
900 }
901 #endif /* RSB_NUMERICAL_TYPE_FLOAT_COMPLEX */
902
903 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
904 if (typecode == RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX)
905 {
906 if(rsb_mm_is_complex(matcode))
907 #if RSB_WANT_OMPIO_SUPPORT
908 rsb__ompio_DOUBLE_COMPLEX(nnz,fd,ngzfd,dval,IA,JA,&re);
909 #else /* RSB_WANT_OMPIO_SUPPORT */
910 for (i=0; i<*nnz; i++)
911 {
912 int iI,iJ;
913 double rv,iv;
914 re += (rsb_zfscanf(fd,"%d %d %lg %lg\n",&iI,&iJ,&(rv),&(iv),ngzfd)==4);
915 (*IA)[i] = (rsb_coo_idx_t)iI;
916 (*JA)[i] = (rsb_coo_idx_t)iJ;
917 (*IA)[i]--; /* adjust from 1-based to 0-based */
918 (*JA)[i]--;
919 (*Zval)[i] = (rv + I * iv);
920 RSB_IO_VERBOSE_MSG(i,*nnz);
921 }
922 #endif /* RSB_WANT_OMPIO_SUPPORT */
923 else
924 for (i=0; i<*nnz; i++)
925 {
926 int iI,iJ;
927 double rv;
928 re += (rsb_zfscanf(fd,"%d %d %lg\n",&iI,&iJ,&(rv),NULL,ngzfd)==3);
929 (*IA)[i] = (rsb_coo_idx_t)iI;
930 (*JA)[i] = (rsb_coo_idx_t)iJ;
931 (*IA)[i]--; /* adjust from 1-based to 0-based */
932 (*JA)[i]--;
933 (*Zval)[i] = (rv + I * 0);
934 RSB_IO_VERBOSE_MSG(i,*nnz);
935 }
936 }
937 #endif /* RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX */
938
939 #ifdef RSB_NUMERICAL_TYPE_PATTERN
940 if (typecode == RSB_NUMERICAL_TYPE_PATTERN)
941 #if RSB_WANT_OMPIO_SUPPORT
942 rsb__ompio_PATTERN(nnz,fd,ngzfd,IA,JA,&re);
943 #else /* RSB_WANT_OMPIO_SUPPORT */
944 for (i=0; i<*nnz; i++)
945 {
946 int iI,iJ;
947 re += (rsb_zfscanf(fd,"%d %d\n",&iI,&iJ,NULL,NULL,ngzfd)==2);
948 (*IA)[i] = (rsb_coo_idx_t)iI;
949 (*JA)[i] = (rsb_coo_idx_t)iJ;
950 (*IA)[i]--; /* adjust from 1-based to 0-based */
951 (*JA)[i]--;
952 RSB_IO_VERBOSE_MSG(i,*nnz);
953 }
954 #endif /* RSB_WANT_OMPIO_SUPPORT */
955 #endif /* RSB_NUMERICAL_TYPE_PATTERN */
956
957 if( is_lowerp || is_upperp )
958 {
959 rsb_flags_t flags = RSB_FLAG_NOFLAGS;
960
961 flags = rsb__do_detect_and_add_triangular_flags(*IA,*JA,*nnz,flags);
962 if( is_lowerp )
963 *is_lowerp = RSB_DO_FLAG_HAS(flags,RSB_FLAG_LOWER);
964 if( is_upperp )
965 *is_upperp = RSB_DO_FLAG_HAS(flags,RSB_FLAG_UPPER);
966 }
967
968 #ifdef RSB_NUMERICAL_TYPE_PATTERN
969 if(typecode == RSB_NUMERICAL_TYPE_PATTERN && otype!=typecode && VA)
970 rsb__fill_with_ones(*VA,otype,*nnz,1);
971 #endif /* RSB_NUMERICAL_TYPE_PATTERN */
972 scan_done:
973 if ( (fd !=stdin ) || ngzfd )
974 {
975 if(ngzfd)
976 RSB_FCLOSE(ngzfd);
977 else
978 RSB_FCLOSE(fd);
979 }
980 if(re!=*nnz)
981 {
982 /* FIXME : this can happen when reading as double a complex matrix, now. */
983 RSB_STDERR("read only %zu out of %d matrix elements (incomplete or not a matrix file?)!\n",re,*nnz);
984 goto err;
985 }
986
987 if( _m != _k && is_symmetric )
988 {
989 RSB_STDERR("matrix declared as symmetric but not square!\n");
990 goto err;
991 }
992
993 if(!(IA && JA))
994 goto afterpmtxchecks;
995 if(is_symmetric && RSB_EXPERIMENTAL_EXPAND_SYMMETRIC_MATRICES_BY_DEFAULT)
996 {
997 /* FIXME: this breaks reallocation tricks !! */
998 if( rsb__reallocate_with_symmetry(IA,JA,VA,nnz,(otype)) )
999 {
1000 RSB_STDERR("problems handling matrix symmetry!\n");
1001 goto err;
1002 }
1003 }
1004
1005 if(is_symmetric && !RSB_EXPERIMENTAL_EXPAND_SYMMETRIC_MATRICES_BY_DEFAULT)
1006 {
1007 rsb_bool_t has_diagonal_elements = RSB_BOOL_FALSE;
1008
1009 if(rsb__util_coo_check_if_triangle_non_empty(*IA,*JA,*nnz,RSB_FLAG_UPPER))
1010 {
1011 RSB_STDERR("#converting upper to lower triangle..\n");
1012 rsb__util_coo_upper_to_lower_symmetric(*IA,*JA,*nnz);
1013 if(is_lower)
1014 is_lower = RSB_BOOL_TRUE;
1015 if( is_lowerp )
1016 *is_lowerp = RSB_BOOL_TRUE;
1017 if( is_upperp )
1018 *is_upperp = RSB_BOOL_FALSE;
1019 }
1020
1021 if(rsb__util_coo_check_if_triangle_non_empty(*IA,*JA,*nnz,RSB_FLAG_UPPER))
1022 {
1023 RSB_STDERR("input declared as symmetric, but it is unsymmetric!\n");
1024 goto err;
1025 }
1026
1027 errval = rsb__util_coo_check_if_has_diagonal_elements(*IA,*JA,*nnz,*m,&has_diagonal_elements);
1028 if(RSB_SOME_ERROR(errval))
1029 {
1030 RSB_STDERR("error while checking diagonal elements!\n");
1031 goto err;
1032 }
1033
1034 if(!has_diagonal_elements)
1035 {
1036 RSB_STDERR("Input has missing elements on the diagonal.\n"); /* FIXME: emit this in a verbose mode only */
1037 }
1038 }
1039 afterpmtxchecks:
1040 frt += rsb_time();
1041 /* this should be a disk read only routine, not an output reporting one */
1042 if(0)
1043 RSB_IO_NOTICE("file I/O took %lf s (%lf nnz, %lf nnz/s ) \n",frt,((double)*nnz),(((double)*nnz)/frt));
1044
1045 if(IA && JA)
1046 if(RSB_DO_FLAG_HAS(flags,RSB_FLAG_UPPER))
1047 RSB_SWAP(rsb_coo_idx_t*,*IA,*JA);
1048
1049 errval = RSB_ERR_NO_ERROR;
1050 prerr:
1051 goto ret;
1052 err:
1053 #if RSB_20120321_IOBUFFERING
1054 RSB_CONDITIONAL_FREE(iobuf);
1055 #endif /* RSB_20120321_IOBUFFERING */
1056 /* FIXME: this will free also already allocated arrays */
1057 if(!aia) if(IA)
1058 RSB_CONDITIONAL_FREE(*IA);
1059 if(!aja) if(JA)
1060 RSB_CONDITIONAL_FREE(*JA);
1061 if(!ava){
1062 #ifdef RSB_NUMERICAL_TYPE_DOUBLE
1063 if(dval) RSB_CONDITIONAL_FREE(*dval);
1064 #endif /* RSB_NUMERICAL_TYPE_DOUBLE */
1065 #ifdef RSB_NUMERICAL_TYPE_FLOAT
1066 if(fval) RSB_CONDITIONAL_FREE(*fval);
1067 #endif /* RSB_NUMERICAL_TYPE_FLOAT */
1068 #ifdef RSB_NUMERICAL_TYPE_INT
1069 if(ival) RSB_CONDITIONAL_FREE(*ival);
1070 #endif /* RSB_NUMERICAL_TYPE_INT */
1071 #ifdef RSB_NUMERICAL_TYPE_CHAR
1072 if(cval) RSB_CONDITIONAL_FREE(*cval);
1073 #endif /* RSB_NUMERICAL_TYPE_CHAR */
1074 #ifdef RSB_NUMERICAL_TYPE_FLOAT_COMPLEX
1075 if(zval) RSB_CONDITIONAL_FREE(*zval);
1076 #endif /* RSB_NUMERICAL_TYPE_FLOAT_COMPLEX */
1077 #ifdef RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX
1078 if(Zval) RSB_CONDITIONAL_FREE(*Zval);
1079 #endif /* RSB_NUMERICAL_TYPE_DOUBLE_COMPLEX */
1080 }
1081 errval = RSB_ERR_GENERIC_ERROR;
1082 ret:
1083 RSB_DO_ERR_RETURN(errval);
1084 }
1085 #endif /* RSB_WITH_MM */
1086
1087 rsb_err_t rsb__do_util_get_matrix_dimensions(const char * filename, size_t * cols, size_t * rows, size_t * nnzp, rsb_flags_t*flagsp)
1088 {
1089 /**
1090 * \ingroup gr_internals
1091 *
1092 * FIXME : needs error handling
1093 */
1094 rsb_coo_idx_t m,k;
1095 rsb_nnz_idx_t nnz;
1096 rsb_flags_t flags = RSB_FLAG_NOFLAGS;
1097 rsb_type_t typecode = RSB_NUMERICAL_TYPE_DEFAULT;
1098 rsb_bool_t is_symmetric = RSB_BOOL_FALSE;
1099 rsb_bool_t is_hermitian = RSB_BOOL_FALSE;
1100 rsb_bool_t is_pattern = RSB_BOOL_FALSE;
1101 rsb_bool_t is_lower = RSB_BOOL_FALSE;
1102 rsb_bool_t is_upper = RSB_BOOL_FALSE;
1103 rsb_bool_t is_vector = RSB_BOOL_FALSE;
1104 rsb_err_t errval = RSB_ERR_NO_ERROR;
1105
1106 errval = rsb__util_mm_info_matrix_f(filename,&m,&k,&nnz,&typecode,&is_symmetric,&is_hermitian,&is_pattern,&is_lower,&is_upper,&is_vector);
1107 if(cols)
1108 *cols = k;
1109 if(rows)
1110 *rows = m;
1111 if(nnzp)
1112 *nnzp = nnz;
1113 if(is_symmetric)
1114 RSB_DO_FLAG_ADD(flags,RSB_FLAG_SYMMETRIC);
1115 if(is_hermitian)
1116 RSB_DO_FLAG_ADD(flags,RSB_FLAG_HERMITIAN);
1117 /* if(is_pattern) ... */
1118 if(is_lower)
1119 RSB_DO_FLAG_ADD(flags,RSB_FLAG_LOWER);
1120 if(is_upper)
1121 RSB_DO_FLAG_ADD(flags,RSB_FLAG_UPPER);
1122 if(flagsp)
1123 *(flagsp) = flags;
1124 RSB_DO_ERR_RETURN(errval)
1125 }
1126
1127 rsb_err_t rsb__util_mm_load_matrix_f_as_csr(const char *filename, rsb_nnz_idx_t ** INDX, rsb_coo_idx_t ** JA, void **VA, rsb_coo_idx_t *m, rsb_coo_idx_t *k , rsb_nnz_idx_t *nnz, rsb_type_t typecode, rsb_flags_t flags/*, rsb_bool_t *is_lowerp, rsb_bool_t *is_upperp*/)
1128 {
1129 /**
1130 * \ingroup gr_internals
1131 * FIXME : should optimize
1132 * */
1133 rsb_err_t errval = RSB_ERR_NO_ERROR;
1134 struct rsb_coo_matrix_t coo;
1135
1136 RSB_DO_FLAG_ADD(flags,RSB_FLAG_WANT_BCSS_STORAGE);
1137 *INDX = NULL;
1138 *JA = NULL;
1139 *VA = NULL;
1140 RSB_BZERO_P(&coo);
1141 coo.typecode = typecode;
1142 errval = rsb_util_mm_load_coo_matrix(filename,&coo);
1143 if(RSB_SOME_ERROR(errval))
1144 goto err;
1145 errval = rsb__util_sort_row_major_inner(coo.VA,coo.IA,coo.JA,coo.nnz,coo.nr,coo.nc,coo.typecode,flags);
1146 if(RSB_SOME_ERROR(errval))
1147 goto err;
1148 errval = rsb__allocate_csr_arrays_from_coo_sorted(coo.VA,coo.IA,coo.JA,coo.nnz,coo.nr,coo.nc,coo.typecode,VA,JA,INDX);
1149 if(RSB_SOME_ERROR(errval))
1150 goto err;
1151 *m = coo.nr; *k = coo.nc; *nnz = coo.nnz;
1152 err:
1153 rsb__destroy_coo_matrix_t(&coo);
1154 RSB_DO_ERR_RETURN(errval)
1155 }
1156
1157 rsb_err_t rsb__util_mm_load_matrix_f_as_csc(const char *filename, rsb_nnz_idx_t ** INDX, rsb_coo_idx_t ** IA, void **VA, rsb_coo_idx_t *m, rsb_coo_idx_t *k , rsb_nnz_idx_t *nnz, rsb_type_t typecode, rsb_flags_t flags/*, rsb_bool_t *is_lowerp, rsb_bool_t *is_upperp*/)
1158 {
1159 /**
1160 * \ingroup gr_internals
1161 * FIXME : should optimize
1162 * */
1163 rsb_err_t errval = RSB_ERR_NO_ERROR;
1164 struct rsb_coo_matrix_t coo;
1165
1166 RSB_DO_FLAG_ADD(flags,RSB_FLAG_WANT_BCSS_STORAGE);
1167 *INDX = NULL;*IA = NULL;*VA = NULL;
1168 RSB_BZERO_P(&coo);
1169 coo.typecode = typecode;
1170 errval = rsb_util_mm_load_coo_matrix(filename,&coo);
1171 if(RSB_SOME_ERROR(errval))
1172 goto err;
1173 errval = rsb_util_sort_column_major(coo.VA,coo.IA,coo.JA,coo.nnz,coo.nr,coo.nc,coo.typecode,flags);
1174 if(RSB_SOME_ERROR(errval))
1175 goto err;
1176 errval = rsb__allocate_csc_arrays_from_coo_sorted(coo.VA,coo.IA,coo.JA,coo.nnz,coo.nr,coo.nc,coo.typecode,VA,IA,INDX);
1177 if(RSB_SOME_ERROR(errval))
1178 goto err;
1179 *m = coo.nr; *k = coo.nc; *nnz = coo.nnz;
1180 err:
1181 rsb__destroy_coo_matrix_t(&coo);
1182 RSB_DO_ERR_RETURN(errval)
1183 }
1184
1185 rsb_err_t rsb_util_mm_fill_arrays_for_csc(const char *filename, rsb_nnz_idx_t * INDX, rsb_coo_idx_t * IA, void *VA, rsb_type_t typecode, rsb_flags_t flags)
1186 {
1187 /**
1188 * \ingroup gr_internals
1189 * FIXME : should optimize
1190 * */
1191 rsb_err_t errval = RSB_ERR_NO_ERROR;
1192 struct rsb_coo_matrix_t coo;
1193 rsb_nnz_idx_t *iINDX = NULL;
1194 rsb_coo_idx_t *iIA = NULL;
1195 void *iVA = NULL;
1196
1197 RSB_DO_FLAG_ADD(flags,RSB_FLAG_WANT_BCSS_STORAGE);
1198 if(!INDX || !IA || !VA)
1199 return RSB_ERR_BADARGS;
1200
1201 RSB_BZERO_P(&coo);
1202 coo.typecode = typecode;
1203 errval = rsb_util_mm_load_coo_matrix(filename,&coo);
1204 if(RSB_SOME_ERROR(errval))
1205 goto err;
1206 errval = rsb_util_sort_column_major(coo.VA,coo.IA,coo.JA,coo.nnz,coo.nr,coo.nc,coo.typecode,flags);
1207 if(RSB_SOME_ERROR(errval))
1208 goto err;
1209 errval = rsb__allocate_csc_arrays_from_coo_sorted(coo.VA,coo.IA,coo.JA,coo.nnz,coo.nr,coo.nc,coo.typecode,&iVA,&iIA,&iINDX);
1210 if(RSB_SOME_ERROR(errval))
1211 goto err;
1212 errval = rsb__copy_css_arrays(iVA,iINDX,iIA,coo.nnz,coo.nc,typecode,VA,INDX,IA);
1213 if(RSB_SOME_ERROR(errval))
1214 goto err;
1215 err:
1216 rsb__destroy_coo_matrix_t(&coo);
1217 RSB_DO_ERR_RETURN(errval)
1218 }
1219
1220 size_t rsb__sys_filesize(const char *filename)
1221 {
1222 /**
1223 * file size, in bytes
1224 * FIXME
1225 * TODO : to sys.c
1226 * TODO: rsb__do_util_get_matrix_dimensions shall invoke this.
1227 * */
1228 #ifdef RSB_HAVE_SYS_STAT_H
1229 struct stat ss;
1230 #endif /* RSB_HAVE_SYS_STAT_H */
1231 if(!filename)
1232 goto err;
1233 #ifdef RSB_HAVE_SYS_STAT_H
1234 stat(filename,&ss);
1235 return ss.st_size;
1236 #else /* RSB_HAVE_SYS_STAT_H */
1237 #endif /* RSB_HAVE_SYS_STAT_H */
1238 err:
1239 return 0;
1240 }
1241
1242 rsb_err_t rsb__do_file_mtx_get_dims(const char * filename, rsb_coo_idx_t* nrp, rsb_coo_idx_t *ncp, rsb_coo_idx_t *nzp, rsb_flags_t*flagsp)
1243 {
1244 rsb_err_t errval = RSB_ERR_NO_ERROR;
1245 size_t nrA = 0,ncA = 0,nzA = 0;
1246 rsb_flags_t flags = RSB_FLAG_NOFLAGS;
1247
1248 errval = rsb__do_util_get_matrix_dimensions(filename,&ncA,&nrA,&nzA,&flags);
1249
1250 if(RSB_SOME_ERROR(errval))
1251 goto ret;
1252
1253 /* RSB_STDOUT("%d / %d %d / %d %d / %d\n",nrA,RSB_MAX_MATRIX_DIM,ncA,RSB_MAX_MATRIX_DIM,nnzA,RSB_MAX_MATRIX_NNZ); */
1254 if( nrp && RSB_INVALID_COO_INDEX(nrA) )
1255 errval |= RSB_ERR_LIMITS;
1256 if( ncp && RSB_INVALID_COO_INDEX(ncA) )
1257 errval |= RSB_ERR_LIMITS;
1258 if( nzp && RSB_INVALID_NNZ_INDEX(nzA) )
1259 errval |= RSB_ERR_LIMITS;
1260 if(nrp)
1261 *nrp = (rsb_coo_idx_t)nrA;
1262 if(ncp)
1263 *ncp = (rsb_coo_idx_t)ncA;
1264 if(nzp)
1265 *nzp = (rsb_nnz_idx_t)nzA;
1266 if(flagsp)
1267 *flagsp = (rsb_flags_t)flags;
1268 /* FIXME: need overflow check here */
1269 ret:
1270 return errval;
1271 }
1272
1273 /* @endcond */
1274