GetFEM  5.5
gmm_inoutput.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2003-2026 Yves Renard, Julien Pommier
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program. If not, see https://www.gnu.org/licenses/.
19 
20  As a special exception, you may use this file as it is a part of a free
21  software library without restriction. Specifically, if other files
22  instantiate templates or use macros or inline functions from this file,
23  or you compile this file and link it with other files to produce an
24  executable, this file does not by itself cause the resulting executable
25  to be covered by the GNU Lesser General Public License. This exception
26  does not however invalidate any other reasons why the executable file
27  might be covered by the GNU Lesser General Public License.
28 
29 ===========================================================================*/
30 
31 /**@file gmm_inoutput.h
32  @author Yves Renard <[email protected]>
33  @author Julien Pommier <[email protected]>
34  @date July 8, 2003.
35  @brief Input/output on sparse matrices
36 
37  Support Harwell-Boeing and Matrix-Market formats.
38 */
39 #ifndef GMM_INOUTPUT_H
40 #define GMM_INOUTPUT_H
41 
42 #include <stdio.h>
43 #include "gmm_kernel.h"
44 namespace gmm {
45 
46  /*************************************************************************/
47  /* */
48  /* Functions to read and write Harwell Boeing format. */
49  /* */
50  /*************************************************************************/
51 
52  // Fri Aug 15 16:29:47 EDT 1997
53  //
54  // Harwell-Boeing File I/O in C
55  // V. 1.0
56  //
57  // National Institute of Standards and Technology, MD.
58  // K.A. Remington
59  // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
60  // NOTICE
61  //
62  // Permission to use, copy, modify, and distribute this software and
63  // its documentation for any purpose and without fee is hereby granted
64  // provided that the above copyright notice appear in all copies and
65  // that both the copyright notice and this permission notice appear in
66  // supporting documentation.
67  //
68  // Neither the Author nor the Institution (National Institute of Standards
69  // and Technology) make any representations about the suitability of this
70  // software for any purpose. This software is provided "as is" without
71  // expressed or implied warranty.
72  // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
73 
74  inline void IOHBTerminate(const char *a) { GMM_ASSERT1(false, a);}
75 
76  inline bool is_complex_double__(std::complex<double>) { return true; }
77  inline bool is_complex_double__(double) { return false; }
78 
79  inline int ParseIfmt(const char *fmt, int* perline, int* width) {
80  if (SECURE_NONCHAR_SSCANF(fmt, " (%dI%d)", perline, width) != 2) {
81  *perline = 1;
82  int s = SECURE_NONCHAR_SSCANF(fmt, " (I%d)", width);
83  GMM_ASSERT1(s == 1, "invalid HB I-format: " << fmt);
84  }
85  return *width;
86  }
87 
88  inline int ParseRfmt(const char *fmt, int* perline, int* width,
89  int* prec, int* flag) {
90  char p;
91  *perline = *width = *flag = *prec = 0;
92 #ifdef GMM_SECURE_CRT
93  if (sscanf_s(fmt, " (%d%c%d.%d)", perline, &p, sizeof(char), width, prec)
94  < 3 || !strchr("PEDF", p))
95 #else
96  if (sscanf(fmt, " (%d%c%d.%d)", perline, &p, width, prec) < 3
97  || !strchr("PEDF", p))
98 #endif
99  {
100  *perline = 1;
101 #ifdef GMM_SECURE_CRT
102  int s = sscanf_s(fmt, " (%c%d.%d)", &p, sizeof(char), width, prec);
103 #else
104  int s = sscanf(fmt, " (%c%d.%d)", &p, width, prec);
105 #endif
106  GMM_ASSERT1(s>=2 && strchr("PEDF",p), "invalid HB REAL format: " << fmt);
107  }
108  *flag = p;
109  return *width;
110  }
111 
112  /** matrix input/output for Harwell-Boeing format */
113  struct HarwellBoeing_IO {
114  int nrows() const { return Nrow; }
115  int ncols() const { return Ncol; }
116  int nnz() const { return Nnzero; }
117  int is_complex() const { return Type[0] == 'C'; }
118  int is_symmetric() const { return Type[1] == 'S'; }
119  int is_hermitian() const { return Type[1] == 'H'; }
120  HarwellBoeing_IO() { clear(); }
121  HarwellBoeing_IO(const char *filename) { clear(); open(filename); }
122  ~HarwellBoeing_IO() { close(); }
123  /** open filename and reads header */
124  void open(const char *filename);
125  /** read the opened file */
126  template <typename T, typename IND_TYPE, int shift> void read(csc_matrix<T, IND_TYPE, shift>& A);
127  template <typename MAT> void read(MAT &M) IS_DEPRECATED;
128  template <typename T, typename IND_TYPE, int shift>
129  static void write(const char *filename, const csc_matrix<T, IND_TYPE, shift>& A);
130  template <typename T, typename IND_TYPE, int shift>
131  static void write(const char *filename, const csc_matrix<T, IND_TYPE, shift>& A,
132  const std::vector<T> &rhs);
133  template <typename T, typename INDI, typename INDJ, int shift>
134  static void write(const char *filename,
135  const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A);
136  template <typename T, typename INDI, typename INDJ, int shift>
137  static void write(const char *filename,
138  const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A,
139  const std::vector<T> &rhs);
140 
141  /** static method for saving the matrix */
142  template <typename MAT> static void write(const char *filename,
143  const MAT& A) IS_DEPRECATED;
144  private:
145  FILE *f;
146  char Title[73], Key[9], Rhstype[4], Type[4];
147  int Nrow, Ncol, Nnzero, Nrhs;
148  char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
149  int Ptrcrd, Indcrd, Valcrd, Rhscrd;
150  int lcount;
151 
152 
153  void close() { if (f) fclose(f); clear(); }
154  void clear() {
155  Nrow = Ncol = Nnzero = Nrhs = 0; f = 0; lcount = 0;
156  memset(Type, 0, sizeof Type);
157  memset(Key, 0, sizeof Key);
158  memset(Title, 0, sizeof Title);
159  }
160  char *getline(char *buf) {
161  char *p = fgets(buf, BUFSIZ, f); ++lcount;
162  int s = SECURE_NONCHAR_SSCANF(buf,"%*s");
163  GMM_ASSERT1(s >= 0 && p != 0,
164  "blank line in HB file at line " << lcount);
165  return buf;
166  }
167 
168  int substrtoi(const char *p, size_type len) {
169  char s[100]; len = std::min(len, sizeof s - 1);
170  SECURE_STRNCPY(s, 100, p, len); s[len] = 0; return atoi(s);
171  }
172  double substrtod(const char *p, size_type len, int Valflag) {
173  char s[100]; len = std::min(len, sizeof s - 1);
174  SECURE_STRNCPY(s, 100, p, len); s[len] = 0;
175  if ( Valflag != 'F' && !strchr(s,'E')) {
176  /* insert a char prefix for exp */
177  int last = int(strlen(s));
178  for (int j=last+1;j>=0;j--) {
179  s[j] = s[j-1];
180  if ( s[j] == '+' || s[j] == '-' ) {
181  s[j-1] = char(Valflag);
182  break;
183  }
184  }
185  }
186  return atof(s);
187  }
188  template <typename IND_TYPE>
189  int readHB_data(IND_TYPE colptr[], IND_TYPE rowind[],
190  double val[]) {
191  /***********************************************************************/
192  /* This function opens and reads the specified file, interpreting its */
193  /* contents as a sparse matrix stored in the Harwell/Boeing standard */
194  /* format and creating compressed column storage scheme vectors to */
195  /* hold the index and nonzero value information. */
196  /* */
197  /* ---------- */
198  /* **CAVEAT** */
199  /* ---------- */
200  /* Parsing real formats from Fortran is tricky, and this file reader */
201  /* does not claim to be foolproof. It has been tested for cases */
202  /* when the real values are printed consistently and evenly spaced on */
203  /* each line, with Fixed (F), and Exponential (E or D) formats. */
204  /* */
205  /* ** If the input file does not adhere to the H/B format, the ** */
206  /* ** results will be unpredictable. ** */
207  /* */
208  /***********************************************************************/
209  int i,ind,col,offset,count;
210  int Ptrperline, Ptrwidth, Indperline, Indwidth;
211  int Valperline, Valwidth, Valprec, Nentries;
212  int Valflag = 'D'; /* Indicates 'E','D', or 'F' float format */
213  char line[BUFSIZ];
214  gmm::standard_locale sl;
215 
216 
217  /* Parse the array input formats from Line 3 of HB file */
218  ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
219  ParseIfmt(Indfmt,&Indperline,&Indwidth);
220  if ( Type[0] != 'P' ) { /* Skip if pattern only */
221  ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
222  }
223 
224  /* Read column pointer array: */
225  offset = 0; /* if base 0 storage is declared (via macro def), */
226  /* then storage entries are offset by 1 */
227 
228  for (count = 0, i=0;i<Ptrcrd;i++) {
229  getline(line);
230  for (col = 0, ind = 0;ind<Ptrperline;ind++) {
231  if (count > Ncol) break;
232  colptr[count] = substrtoi(line+col,Ptrwidth)-offset;
233  count++; col += Ptrwidth;
234  }
235  }
236 
237  /* Read row index array: */
238  for (count = 0, i=0;i<Indcrd;i++) {
239  getline(line);
240  for (col = 0, ind = 0;ind<Indperline;ind++) {
241  if (count == Nnzero) break;
242  rowind[count] = substrtoi(line+col,Indwidth)-offset;
243  count++; col += Indwidth;
244  }
245  }
246 
247  /* Read array of values: */
248  if ( Type[0] != 'P' ) { /* Skip if pattern only */
249  if ( Type[0] == 'C' ) Nentries = 2*Nnzero;
250  else Nentries = Nnzero;
251 
252  count = 0;
253  for (i=0;i<Valcrd;i++) {
254  getline(line);
255  if (Valflag == 'D') {
256  // const_cast Due to aCC excentricity
257  char *p;
258  while( (p = const_cast<char *>(strchr(line,'D')) )) *p = 'E';
259  }
260  for (col = 0, ind = 0;ind<Valperline;ind++) {
261  if (count == Nentries) break;
262  val[count] = substrtod(line+col, Valwidth, Valflag);
263  count++; col += Valwidth;
264  }
265  }
266  }
267  return 1;
268  }
269  };
270 
271  inline void HarwellBoeing_IO::open(const char *filename) {
272  int Totcrd,Neltvl,Nrhsix;
273  char line[BUFSIZ];
274  close();
275  SECURE_FOPEN(&f, filename, "r");
276  GMM_ASSERT1(f, "could not open " << filename);
277  /* First line: */
278 #ifdef GMM_SECURE_CRT
279  sscanf_s(getline(line), "%c%s", Title, 72, Key, 8);
280 #else
281  sscanf(getline(line), "%72c%8s", Title, Key);
282 #endif
283  Key[8] = Title[72] = 0;
284  /* Second line: */
285  Totcrd = Ptrcrd = Indcrd = Valcrd = Rhscrd = 0;
286  SECURE_NONCHAR_SSCANF(getline(line), "%d%d%d%d%d", &Totcrd, &Ptrcrd,
287  &Indcrd, &Valcrd, &Rhscrd);
288 
289  /* Third line: */
290  Nrow = Ncol = Nnzero = Neltvl = 0;
291 #ifdef GMM_SECURE_CRT
292  if (sscanf_s(getline(line), "%c%d%d%d%d", Type, 3, &Nrow, &Ncol, &Nnzero,
293  &Neltvl) < 1)
294 #else
295  if (sscanf(getline(line), "%3c%d%d%d%d", Type, &Nrow, &Ncol, &Nnzero,
296  &Neltvl) < 1)
297 #endif
298  IOHBTerminate("Invalid Type info, line 3 of Harwell-Boeing file.\n");
299  for (size_type i = 0; i < 3; ++i) Type[i] = char(toupper(Type[i]));
300 
301  /* Fourth line: */
302 #ifdef GMM_SECURE_CRT
303  if ( sscanf_s(getline(line), "%c%c%c%c",Ptrfmt, 16,Indfmt, 16,Valfmt,
304  20,Rhsfmt, 20) < 3)
305 #else
306  if ( sscanf(getline(line), "%16c%16c%20c%20c",Ptrfmt,Indfmt,Valfmt,
307  Rhsfmt) < 3)
308 #endif
309  IOHBTerminate("Invalid format info, line 4 of Harwell-Boeing file.\n");
310  Ptrfmt[16] = Indfmt[16] = Valfmt[20] = Rhsfmt[20] = 0;
311 
312  /* (Optional) Fifth line: */
313  if (Rhscrd != 0 ) {
314  Nrhs = Nrhsix = 0;
315 #ifdef GMM_SECURE_CRT
316  if ( sscanf_s(getline(line), "%c%d%d", Rhstype, 3, &Nrhs, &Nrhsix) != 1)
317 #else
318  if ( sscanf(getline(line), "%3c%d%d", Rhstype, &Nrhs, &Nrhsix) != 1)
319 #endif
320  IOHBTerminate("Invalid RHS type information, line 5 of"
321  " Harwell-Boeing file.\n");
322  }
323  }
324 
325  /* only valid for double and complex<double> csc matrices */
326  template <typename T, typename IND_TYPE, int shift> void
327  HarwellBoeing_IO::read(csc_matrix<T, IND_TYPE, shift>& A) {
328 
329  // typedef typename csc_matrix<T, shift>::IND_TYPE IND_TYPE;
330 
331  GMM_ASSERT1(f, "no file opened!");
332  GMM_ASSERT1(Type[0] != 'P',
333  "Bad HB matrix format (pattern matrices not supported)");
334  GMM_ASSERT1(!is_complex_double__(T()) || Type[0] != 'R',
335  "Bad HB matrix format (file contains a REAL matrix)");
336  GMM_ASSERT1(is_complex_double__(T()) || Type[0] != 'C',
337  "Bad HB matrix format (file contains a COMPLEX matrix)");
338  A.nc = ncols(); A.nr = nrows();
339  A.jc.resize(ncols()+1);
340  A.ir.resize(nnz());
341  A.pr.resize(nnz());
342  readHB_data(&A.jc[0], &A.ir[0], (double*)&A.pr[0]);
343  for (int i = 0; i <= ncols(); ++i) { A.jc[i] += shift; A.jc[i] -= 1; }
344  for (int i = 0; i < nnz(); ++i) { A.ir[i] += shift; A.ir[i] -= 1; }
345  }
346 
347  template <typename MAT> void
348  HarwellBoeing_IO::read(MAT &M) {
349  csc_matrix<typename gmm::linalg_traits<MAT>::value_type> csc;
350  read(csc);
351  resize(M, mat_nrows(csc), mat_ncols(csc));
352  copy(csc, M);
353  }
354 
355  template <typename IND_TYPE>
356  inline int writeHB_mat_double(const char* filename, int M, int N, int nz,
357  const IND_TYPE colptr[],
358  const IND_TYPE rowind[],
359  const double val[], int Nrhs,
360  const double rhs[], const double guess[],
361  const double exact[], const char* Title,
362  const char* Key, const char* Type,
363  const char* Ptrfmt, const char* Indfmt,
364  const char* Valfmt, const char* Rhsfmt,
365  const char* Rhstype, int shift) {
366  /************************************************************************/
367  /* The writeHB function opens the named file and writes the specified */
368  /* matrix and optional right-hand-side(s) to that file in */
369  /* Harwell-Boeing format. */
370  /* */
371  /* For a description of the Harwell Boeing standard, see: */
372  /* Duff, et al., ACM TOMS Vol.15, No.1, March 1989 */
373  /* */
374  /************************************************************************/
375  FILE *out_file;
376  int i, entry, offset, j, acount, linemod;
377  int totcrd, ptrcrd, indcrd, valcrd, rhscrd;
378  int nvalentries, nrhsentries;
379  int Ptrperline, Ptrwidth, Indperline, Indwidth;
380  int Rhsperline, Rhswidth, Rhsprec, Rhsflag;
381  int Valperline, Valwidth, Valprec;
382  int Valflag; /* Indicates 'E','D', or 'F' float format */
383  char pformat[16],iformat[16],vformat[19],rformat[19];
384  // char *pValflag, *pRhsflag;
385  gmm::standard_locale sl;
386 
387  if ( Type[0] == 'C' )
388  { nvalentries = 2*nz; nrhsentries = 2*M; }
389  else
390  { nvalentries = nz; nrhsentries = M; }
391 
392  if ( filename != NULL ) {
393  SECURE_FOPEN(&out_file, filename, "w");
394  GMM_ASSERT1(out_file != NULL, "Error: Cannot open file: " << filename);
395  } else out_file = stdout;
396 
397  if ( Ptrfmt == NULL ) Ptrfmt = "(8I10)";
398  ParseIfmt(Ptrfmt, &Ptrperline, &Ptrwidth);
399  SECURE_SPRINTF1(pformat,sizeof(pformat),"%%%dd",Ptrwidth);
400  ptrcrd = (N+1)/Ptrperline;
401  if ( (N+1)%Ptrperline != 0) ptrcrd++;
402 
403  if ( Indfmt == NULL ) Indfmt = Ptrfmt;
404  ParseIfmt(Indfmt, &Indperline, &Indwidth);
405  SECURE_SPRINTF1(iformat,sizeof(iformat), "%%%dd",Indwidth);
406  indcrd = nz/Indperline;
407  if ( nz%Indperline != 0) indcrd++;
408 
409  if ( Type[0] != 'P' ) { /* Skip if pattern only */
410  if ( Valfmt == NULL ) Valfmt = "(4E21.13)";
411  ParseRfmt(Valfmt, &Valperline, &Valwidth, &Valprec, &Valflag);
412 // if (Valflag == 'D') {
413 // pValflag = (char *) strchr(Valfmt,'D');
414 // *pValflag = 'E';
415 // }
416  if (Valflag == 'F')
417  SECURE_SPRINTF2(vformat, sizeof(vformat), "%% %d.%df", Valwidth,
418  Valprec);
419  else
420  SECURE_SPRINTF2(vformat, sizeof(vformat), "%% %d.%dE", Valwidth,
421  Valprec);
422  valcrd = nvalentries/Valperline;
423  if ( nvalentries%Valperline != 0) valcrd++;
424  } else valcrd = 0;
425 
426  if ( Nrhs > 0 ) {
427  if ( Rhsfmt == NULL ) Rhsfmt = Valfmt;
428  ParseRfmt(Rhsfmt,&Rhsperline,&Rhswidth,&Rhsprec, &Rhsflag);
429  if (Rhsflag == 'F')
430  SECURE_SPRINTF2(rformat,sizeof(rformat), "%% %d.%df",Rhswidth,Rhsprec);
431  else
432  SECURE_SPRINTF2(rformat,sizeof(rformat), "%% %d.%dE",Rhswidth,Rhsprec);
433 // if (Valflag == 'D') {
434 // pRhsflag = (char *) strchr(Rhsfmt,'D');
435 // *pRhsflag = 'E';
436 // }
437  rhscrd = nrhsentries/Rhsperline;
438  if ( nrhsentries%Rhsperline != 0) rhscrd++;
439  if ( Rhstype[1] == 'G' ) rhscrd+=rhscrd;
440  if ( Rhstype[2] == 'X' ) rhscrd+=rhscrd;
441  rhscrd*=Nrhs;
442  } else rhscrd = 0;
443 
444  totcrd = 4+ptrcrd+indcrd+valcrd+rhscrd;
445 
446 
447  /* Print header information: */
448 
449  fprintf(out_file,"%-72s%-8s\n%14d%14d%14d%14d%14d\n",Title, Key, totcrd,
450  ptrcrd, indcrd, valcrd, rhscrd);
451  fprintf(out_file,"%3s%11s%14d%14d%14d%14d\n",Type," ", M, N, nz, 0);
452  fprintf(out_file,"%-16s%-16s%-20s", Ptrfmt, Indfmt, Valfmt);
453  if ( Nrhs != 0 ) {
454  /* Print Rhsfmt on fourth line and */
455  /* optional fifth header line for auxillary vector information:*/
456  fprintf(out_file,"%-20s\n%-14s%d\n",Rhsfmt,Rhstype,Nrhs);
457  }
458  else
459  fprintf(out_file,"\n");
460 
461  offset = 1 - shift; /* if base 0 storage is declared (via macro def), */
462  /* then storage entries are offset by 1 */
463 
464  /* Print column pointers: */
465  for (i = 0; i < N+1; i++) {
466  entry = colptr[i]+offset;
467  fprintf(out_file,pformat,entry);
468  if ( (i+1)%Ptrperline == 0 ) fprintf(out_file,"\n");
469  }
470 
471  if ( (N+1) % Ptrperline != 0 ) fprintf(out_file,"\n");
472 
473  /* Print row indices: */
474  for (i=0;i<nz;i++) {
475  entry = rowind[i]+offset;
476  fprintf(out_file,iformat,entry);
477  if ( (i+1)%Indperline == 0 ) fprintf(out_file,"\n");
478  }
479 
480  if ( nz % Indperline != 0 ) fprintf(out_file,"\n");
481 
482  /* Print values: */
483 
484  if ( Type[0] != 'P' ) { /* Skip if pattern only */
485  for (i=0;i<nvalentries;i++) {
486  fprintf(out_file,vformat,val[i]);
487  if ( (i+1)%Valperline == 0 ) fprintf(out_file,"\n");
488  }
489 
490  if ( nvalentries % Valperline != 0 ) fprintf(out_file,"\n");
491 
492  /* Print right hand sides: */
493  acount = 1;
494  linemod=0;
495  if ( Nrhs > 0 ) {
496  for (j=0;j<Nrhs;j++) {
497  for (i=0;i<nrhsentries;i++) {
498  fprintf(out_file,rformat,rhs[i] /* *Rhswidth */);
499  if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
500  }
501  if ( acount%Rhsperline != linemod ) {
502  fprintf(out_file,"\n");
503  linemod = (acount-1)%Rhsperline;
504  }
505  if ( Rhstype[1] == 'G' ) {
506  for (i=0;i<nrhsentries;i++) {
507  fprintf(out_file,rformat,guess[i] /* *Rhswidth */);
508  if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
509  }
510  if ( acount%Rhsperline != linemod ) {
511  fprintf(out_file,"\n");
512  linemod = (acount-1)%Rhsperline;
513  }
514  }
515  if ( Rhstype[2] == 'X' ) {
516  for (i=0;i<nrhsentries;i++) {
517  fprintf(out_file,rformat,exact[i] /* *Rhswidth */);
518  if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
519  }
520  if ( acount%Rhsperline != linemod ) {
521  fprintf(out_file,"\n");
522  linemod = (acount-1)%Rhsperline;
523  }
524  }
525  }
526  }
527  }
528  int s = fclose(out_file);
529  GMM_ASSERT1(s == 0, "Error closing file in writeHB_mat_double().");
530  return 1;
531  }
532 
533  template <typename T, typename IND_TYPE, int shift> void
534  HarwellBoeing_IO::write(const char *filename,
535  const csc_matrix<T, IND_TYPE, shift>& A) {
536  write(filename, csc_matrix_ref<const T*, const unsigned*,
537  const unsigned *, shift>
538  (&A.pr[0], &A.ir[0], &A.jc[0], A.nr, A.nc));
539  }
540 
541  template <typename T, typename IND_TYPE, int shift> void
542  HarwellBoeing_IO::write(const char *filename,
543  const csc_matrix<T, IND_TYPE, shift>& A,
544  const std::vector<T> &rhs) {
545  write(filename, csc_matrix_ref<const T*, const unsigned*,
546  const unsigned *, shift>
547  (&A.pr[0], &A.ir[0], &A.jc[0], A.nr, A.nc), rhs);
548  }
549 
550  template <typename T, typename INDI, typename INDJ, int shift> void
551  HarwellBoeing_IO::write(const char *filename,
552  const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A) {
553  const char *t = 0;
554  if (is_complex_double__(T()))
555  if (mat_nrows(A) == mat_ncols(A)) t = "CUA"; else t = "CRA";
556  else
557  if (mat_nrows(A) == mat_ncols(A)) t = "RUA"; else t = "RRA";
558  writeHB_mat_double(filename, int(mat_nrows(A)), int(mat_ncols(A)),
559  A.jc[mat_ncols(A)], A.jc, A.ir,
560  (const double *)A.pr,
561  0, 0, 0, 0, "GMM CSC MATRIX", "CSCMAT",
562  t, 0, 0, 0, 0, "F", shift);
563  }
564 
565  template <typename T, typename INDI, typename INDJ, int shift> void
566  HarwellBoeing_IO::write(const char *filename,
567  const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A,
568  const std::vector<T> &rhs) {
569  const char *t = 0;
570  if (is_complex_double__(T()))
571  if (mat_nrows(A) == mat_ncols(A)) t = "CUA"; else t = "CRA";
572  else
573  if (mat_nrows(A) == mat_ncols(A)) t = "RUA"; else t = "RRA";
574  int Nrhs = gmm::vect_size(rhs) / mat_nrows(A);
575  writeHB_mat_double(filename, int(mat_nrows(A)), int(mat_ncols(A)),
576  A.jc[mat_ncols(A)], A.jc, A.ir,
577  (const double *)A.pr,
578  Nrhs, (const double *)(&rhs[0]), 0, 0,
579  "GMM CSC MATRIX", "CSCMAT",
580  t, 0, 0, 0, 0, "F ", shift);
581  }
582 
583 
584  template <typename MAT> void
585  HarwellBoeing_IO::write(const char *filename, const MAT& A) {
586  gmm::csc_matrix<typename gmm::linalg_traits<MAT>::value_type>
587  tmp(gmm::mat_nrows(A), gmm::mat_ncols(A));
588  gmm::copy(A,tmp);
589  HarwellBoeing_IO::write(filename, tmp);
590  }
591 
592  /** save a "double" or "std::complex<double>" csc matrix into a
593  HarwellBoeing file
594  */
595  template <typename T, typename IND_TYPE, int shift> inline void
596  Harwell_Boeing_save(const std::string &filename,
597  const csc_matrix<T, IND_TYPE, shift>& A)
598  { HarwellBoeing_IO::write(filename.c_str(), A); }
599 
600  /** save a reference on "double" or "std::complex<double>" csc matrix
601  into a HarwellBoeing file
602  */
603  template <typename T, typename INDI, typename INDJ, int shift> inline void
604  Harwell_Boeing_save(const std::string &filename,
605  const csc_matrix_ref<T, INDI, INDJ, shift>& A)
606  { HarwellBoeing_IO::write(filename.c_str(), A); }
607 
608  /** save a "double" or "std::complex<double>" generic matrix
609  into a HarwellBoeing file making a copy in a csc matrix
610  */
611  template <typename MAT> inline void
612  Harwell_Boeing_save(const std::string &filename, const MAT& A) {
613  gmm::csc_matrix<typename gmm::linalg_traits<MAT>::value_type>
614  tmp(gmm::mat_nrows(A), gmm::mat_ncols(A));
615  gmm::copy(A, tmp);
616  HarwellBoeing_IO::write(filename.c_str(), tmp);
617  }
618 
619  template <typename MAT, typename VECT> inline void
620  Harwell_Boeing_save(const std::string &filename, const MAT& A,
621  const VECT &RHS) {
622  typedef typename gmm::linalg_traits<MAT>::value_type T;
623  gmm::csc_matrix<T> tmp(gmm::mat_nrows(A), gmm::mat_ncols(A));
624  gmm::copy(A, tmp);
625  std::vector<T> tmprhs(gmm::vect_size(RHS));
626  gmm::copy(RHS, tmprhs);
627  HarwellBoeing_IO::write(filename.c_str(), tmp, tmprhs);
628  }
629 
630  /** load a "double" or "std::complex<double>" csc matrix from a
631  HarwellBoeing file
632  */
633  template <typename T, typename IND_TYPE, int shift> void
634  Harwell_Boeing_load(const std::string &filename, csc_matrix<T, IND_TYPE, shift>& A) {
635  HarwellBoeing_IO h(filename.c_str()); h.read(A);
636  }
637 
638  /** load a "double" or "std::complex<double>" generic matrix from a
639  HarwellBoeing file
640  */
641  template <typename MAT> void
642  Harwell_Boeing_load(const std::string &filename, MAT& A) {
643  csc_matrix<typename gmm::linalg_traits<MAT>::value_type> csc;
644  Harwell_Boeing_load(filename, csc);
645  resize(A, mat_nrows(csc), mat_ncols(csc));
646  copy(csc, A);
647  }
648 
649  /*************************************************************************/
650  /* */
651  /* Functions to read and write MatrixMarket format. */
652  /* */
653  /*************************************************************************/
654 
655  /*
656  * Matrix Market I/O library for ANSI C
657  *
658  * See http://math.nist.gov/MatrixMarket for details.
659  *
660  *
661  */
662 
663 #define MM_MAX_LINE_LENGTH 1025
664 #define MatrixMarketBanner "%%MatrixMarket"
665 #define MM_MAX_TOKEN_LENGTH 64
666 
667  typedef char MM_typecode[4];
668 
669  /******************* MM_typecode query functions *************************/
670 
671 #define mm_is_matrix(typecode) ((typecode)[0]=='M')
672 
673 #define mm_is_sparse(typecode) ((typecode)[1]=='C')
674 #define mm_is_coordinate(typecode) ((typecode)[1]=='C')
675 #define mm_is_dense(typecode) ((typecode)[1]=='A')
676 #define mm_is_array(typecode) ((typecode)[1]=='A')
677 
678 #define mm_is_complex(typecode) ((typecode)[2]=='C')
679 #define mm_is_real(typecode) ((typecode)[2]=='R')
680 #define mm_is_pattern(typecode) ((typecode)[2]=='P')
681 #define mm_is_integer(typecode) ((typecode)[2]=='I')
682 
683 #define mm_is_symmetric(typecode) ((typecode)[3]=='S')
684 #define mm_is_general(typecode) ((typecode)[3]=='G')
685 #define mm_is_skew(typecode) ((typecode)[3]=='K')
686 #define mm_is_hermitian(typecode) ((typecode)[3]=='H')
687 
688  /******************* MM_typecode modify fucntions ************************/
689 
690 #define mm_set_matrix(typecode) ((*typecode)[0]='M')
691 #define mm_set_coordinate(typecode) ((*typecode)[1]='C')
692 #define mm_set_array(typecode) ((*typecode)[1]='A')
693 #define mm_set_dense(typecode) mm_set_array(typecode)
694 #define mm_set_sparse(typecode) mm_set_coordinate(typecode)
695 
696 #define mm_set_complex(typecode) ((*typecode)[2]='C')
697 #define mm_set_real(typecode) ((*typecode)[2]='R')
698 #define mm_set_pattern(typecode) ((*typecode)[2]='P')
699 #define mm_set_integer(typecode) ((*typecode)[2]='I')
700 
701 
702 #define mm_set_symmetric(typecode) ((*typecode)[3]='S')
703 #define mm_set_general(typecode) ((*typecode)[3]='G')
704 #define mm_set_skew(typecode) ((*typecode)[3]='K')
705 #define mm_set_hermitian(typecode) ((*typecode)[3]='H')
706 
707 #define mm_clear_typecode(typecode) ((*typecode)[0]=(*typecode)[1]= \
708  (*typecode)[2]=' ',(*typecode)[3]='G')
709 
710 #define mm_initialize_typecode(typecode) mm_clear_typecode(typecode)
711 
712 
713  /******************* Matrix Market error codes ***************************/
714 
715 
716 #define MM_COULD_NOT_READ_FILE 11
717 #define MM_PREMATURE_EOF 12
718 #define MM_NOT_MTX 13
719 #define MM_NO_HEADER 14
720 #define MM_UNSUPPORTED_TYPE 15
721 #define MM_LINE_TOO_LONG 16
722 #define MM_COULD_NOT_WRITE_FILE 17
723 
724 
725  /******************** Matrix Market internal definitions *****************
726 
727  MM_matrix_typecode: 4-character sequence
728 
729  object sparse/ data storage
730  dense type scheme
731 
732  string position: [0] [1] [2] [3]
733 
734  Matrix typecode: M(atrix) C(oord) R(eal) G(eneral)
735  A(array) C(omplex) H(ermitian)
736  P(attern) S(ymmetric)
737  I(nteger) K(kew)
738 
739  ***********************************************************************/
740 
741 #define MM_MTX_STR "matrix"
742 #define MM_ARRAY_STR "array"
743 #define MM_DENSE_STR "array"
744 #define MM_COORDINATE_STR "coordinate"
745 #define MM_SPARSE_STR "coordinate"
746 #define MM_COMPLEX_STR "complex"
747 #define MM_REAL_STR "real"
748 #define MM_INT_STR "integer"
749 #define MM_GENERAL_STR "general"
750 #define MM_SYMM_STR "symmetric"
751 #define MM_HERM_STR "hermitian"
752 #define MM_SKEW_STR "skew-symmetric"
753 #define MM_PATTERN_STR "pattern"
754 
755  inline char *mm_typecode_to_str(MM_typecode matcode) {
756  char buffer[MM_MAX_LINE_LENGTH];
757  const char *types[4] = {0,0,0,0};
758  /* int error =0; */
759  /* int i; */
760 
761  /* check for MTX type */
762  if (mm_is_matrix(matcode))
763  types[0] = MM_MTX_STR;
764  /*
765  else
766  error=1;
767  */
768  /* check for CRD or ARR matrix */
769  if (mm_is_sparse(matcode))
770  types[1] = MM_SPARSE_STR;
771  else
772  if (mm_is_dense(matcode))
773  types[1] = MM_DENSE_STR;
774  else
775  return NULL;
776 
777  /* check for element data type */
778  if (mm_is_real(matcode))
779  types[2] = MM_REAL_STR;
780  else
781  if (mm_is_complex(matcode))
782  types[2] = MM_COMPLEX_STR;
783  else
784  if (mm_is_pattern(matcode))
785  types[2] = MM_PATTERN_STR;
786  else
787  if (mm_is_integer(matcode))
788  types[2] = MM_INT_STR;
789  else
790  return NULL;
791 
792 
793  /* check for symmetry type */
794  if (mm_is_general(matcode))
795  types[3] = MM_GENERAL_STR;
796  else if (mm_is_symmetric(matcode))
797  types[3] = MM_SYMM_STR;
798  else if (mm_is_hermitian(matcode))
799  types[3] = MM_HERM_STR;
800  else if (mm_is_skew(matcode))
801  types[3] = MM_SKEW_STR;
802  else
803  return NULL;
804 
805  SECURE_SPRINTF4(buffer, sizeof(buffer), "%s %s %s %s", types[0], types[1],
806  types[2], types[3]);
807  return SECURE_STRDUP(buffer);
808 
809  }
810 
811  inline int mm_read_banner(FILE *f, MM_typecode *matcode) {
812  char line[MM_MAX_LINE_LENGTH];
813  char banner[MM_MAX_TOKEN_LENGTH];
814  char mtx[MM_MAX_TOKEN_LENGTH];
815  char crd[MM_MAX_TOKEN_LENGTH];
816  char data_type[MM_MAX_TOKEN_LENGTH];
817  char storage_scheme[MM_MAX_TOKEN_LENGTH];
818  char *p;
819  gmm::standard_locale sl;
820  /* int ret_code; */
821 
822  mm_clear_typecode(matcode);
823 
824  if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
825  return MM_PREMATURE_EOF;
826 
827 #ifdef GMM_SECURE_CRT
828  if (sscanf_s(line, "%s %s %s %s %s", banner, sizeof(banner),
829  mtx, sizeof(mtx), crd, sizeof(crd), data_type,
830  sizeof(data_type), storage_scheme,
831  sizeof(storage_scheme)) != 5)
832 #else
833  if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd,
834  data_type, storage_scheme) != 5)
835 #endif
836  return MM_PREMATURE_EOF;
837 
838  for (p=mtx; *p!='\0'; *p=char(tolower(*p)),p++) {}; /* convert to lower case */
839  for (p=crd; *p!='\0'; *p=char(tolower(*p)),p++) {};
840  for (p=data_type; *p!='\0'; *p=char(tolower(*p)),p++) {};
841  for (p=storage_scheme; *p!='\0'; *p=char(tolower(*p)),p++) {};
842 
843  /* check for banner */
844  if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
845  return MM_NO_HEADER;
846 
847  /* first field should be "mtx" */
848  if (strcmp(mtx, MM_MTX_STR) != 0)
849  return MM_UNSUPPORTED_TYPE;
850  mm_set_matrix(matcode);
851 
852 
853  /* second field describes whether this is a sparse matrix (in coordinate
854  storgae) or a dense array */
855 
856 
857  if (strcmp(crd, MM_SPARSE_STR) == 0)
858  mm_set_sparse(matcode);
859  else
860  if (strcmp(crd, MM_DENSE_STR) == 0)
861  mm_set_dense(matcode);
862  else
863  return MM_UNSUPPORTED_TYPE;
864 
865 
866  /* third field */
867 
868  if (strcmp(data_type, MM_REAL_STR) == 0)
869  mm_set_real(matcode);
870  else
871  if (strcmp(data_type, MM_COMPLEX_STR) == 0)
872  mm_set_complex(matcode);
873  else
874  if (strcmp(data_type, MM_PATTERN_STR) == 0)
875  mm_set_pattern(matcode);
876  else
877  if (strcmp(data_type, MM_INT_STR) == 0)
878  mm_set_integer(matcode);
879  else
880  return MM_UNSUPPORTED_TYPE;
881 
882 
883  /* fourth field */
884 
885  if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
886  mm_set_general(matcode);
887  else
888  if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
889  mm_set_symmetric(matcode);
890  else
891  if (strcmp(storage_scheme, MM_HERM_STR) == 0)
892  mm_set_hermitian(matcode);
893  else
894  if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
895  mm_set_skew(matcode);
896  else
897  return MM_UNSUPPORTED_TYPE;
898 
899  return 0;
900  }
901 
902  inline int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz ) {
903  char line[MM_MAX_LINE_LENGTH];
904  /* int ret_code;*/
905  int num_items_read;
906 
907  /* set return null parameter values, in case we exit with errors */
908  *M = *N = *nz = 0;
909 
910  /* now continue scanning until you reach the end-of-comments */
911  do {
912  if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
913  return MM_PREMATURE_EOF;
914  } while (line[0] == '%');
915 
916  /* line[] is either blank or has M,N, nz */
917  if (SECURE_NONCHAR_SSCANF(line, "%d %d %d", M, N, nz) == 3) return 0;
918  else
919  do {
920  num_items_read = SECURE_NONCHAR_FSCANF(f, "%d %d %d", M, N, nz);
921  if (num_items_read == EOF) return MM_PREMATURE_EOF;
922  }
923  while (num_items_read != 3);
924 
925  return 0;
926  }
927 
928 
929  inline int mm_read_mtx_crd_data(FILE *f, int, int, int nz, int II[],
930  int J[], double val[], MM_typecode matcode) {
931  int i;
932  if (mm_is_complex(matcode)) {
933  for (i=0; i<nz; i++)
934  if (SECURE_NONCHAR_FSCANF(f, "%d %d %lg %lg", &II[i], &J[i],
935  &val[2*i], &val[2*i+1])
936  != 4) return MM_PREMATURE_EOF;
937  }
938  else if (mm_is_real(matcode)) {
939  for (i=0; i<nz; i++) {
940  if (SECURE_NONCHAR_FSCANF(f, "%d %d %lg\n", &II[i], &J[i], &val[i])
941  != 3) return MM_PREMATURE_EOF;
942 
943  }
944  }
945  else if (mm_is_pattern(matcode)) {
946  for (i=0; i<nz; i++)
947  if (SECURE_NONCHAR_FSCANF(f, "%d %d", &II[i], &J[i])
948  != 2) return MM_PREMATURE_EOF;
949  }
950  else return MM_UNSUPPORTED_TYPE;
951 
952  return 0;
953  }
954 
955  inline int mm_write_mtx_crd(const char *fname, int M, int N, int nz,
956  int II[], int J[], const double val[],
957  MM_typecode matcode) {
958  FILE *f;
959  int i;
960 
961  if (strcmp(fname, "stdout") == 0)
962  f = stdout;
963  else {
964  SECURE_FOPEN(&f, fname, "w");
965  if (f == NULL)
966  return MM_COULD_NOT_WRITE_FILE;
967  }
968 
969  /* print banner followed by typecode */
970  fprintf(f, "%s ", MatrixMarketBanner);
971  char *str = mm_typecode_to_str(matcode);
972  fprintf(f, "%s\n", str);
973  free(str);
974 
975  /* print matrix sizes and nonzeros */
976  fprintf(f, "%d %d %d\n", M, N, nz);
977 
978  /* print values */
979  if (mm_is_pattern(matcode))
980  for (i=0; i<nz; i++)
981  fprintf(f, "%d %d\n", II[i], J[i]);
982  else
983  if (mm_is_real(matcode))
984  for (i=0; i<nz; i++)
985  fprintf(f, "%d %d %20.16g\n", II[i], J[i], val[i]);
986  else
987  if (mm_is_complex(matcode))
988  for (i=0; i<nz; i++)
989  fprintf(f, "%d %d %20.16g %20.16g\n", II[i], J[i], val[2*i],
990  val[2*i+1]);
991  else {
992  if (f != stdout) fclose(f);
993  return MM_UNSUPPORTED_TYPE;
994  }
995 
996  if (f !=stdout) fclose(f);
997  return 0;
998  }
999 
1000 
1001  /** matrix input/output for MatrixMarket storage */
1002  class MatrixMarket_IO {
1003  FILE *f;
1004  bool isComplex, isSymmetric, isHermitian;
1005  int row, col, nz;
1006  MM_typecode matcode;
1007  public:
1008  MatrixMarket_IO() : f(0) {}
1009  MatrixMarket_IO(const char *filename) : f(0) { open(filename); }
1010  ~MatrixMarket_IO() { if (f) fclose(f); f = 0; }
1011 
1012  int nrows() const { return row; }
1013  int ncols() const { return col; }
1014  int nnz() const { return nz; }
1015  int is_complex() const { return isComplex; }
1016  int is_symmetric() const { return isSymmetric; }
1017  int is_hermitian() const { return isHermitian; }
1018 
1019  /* open filename and reads header */
1020  void open(const char *filename);
1021  /* read opened file */
1022  template <typename Matrix> void read(Matrix &A);
1023  /* write a matrix */
1024  template <typename T, typename IND_TYPE, int shift> static void
1025  write(const char *filename, const csc_matrix<T, IND_TYPE, shift>& A);
1026  template <typename T, typename INDI, typename INDJ, int shift> static void
1027  write(const char *filename,
1028  const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A);
1029  template <typename MAT> static void
1030  write(const char *filename, const MAT& A);
1031  };
1032 
1033  /** load a matrix-market file */
1034  template <typename Matrix> inline void
1035  MatrixMarket_load(const char *filename, Matrix& A) {
1036  MatrixMarket_IO mm; mm.open(filename);
1037  mm.read(A);
1038  }
1039  /** write a matrix-market file */
1040  template <typename T, typename IND_TYPE, int shift> void
1041  MatrixMarket_save(const char *filename, const csc_matrix<T, IND_TYPE, shift>& A) {
1042  MatrixMarket_IO mm; mm.write(filename, A);
1043  }
1044 
1045  template <typename T, typename INDI, typename INDJ, int shift> inline void
1046  MatrixMarket_save(const char *filename,
1047  const csc_matrix_ref<T, INDI, INDJ, shift>& A) {
1048  MatrixMarket_IO mm; mm.write(filename, A);
1049  }
1050 
1051 
1052  inline void MatrixMarket_IO::open(const char *filename) {
1053  gmm::standard_locale sl;
1054  if (f) { fclose(f); }
1055  SECURE_FOPEN(&f, filename, "r");
1056  GMM_ASSERT1(f, "Sorry, cannot open file " << filename);
1057  int s1 = mm_read_banner(f, &matcode);
1058  GMM_ASSERT1(s1 == 0, "Sorry, cannnot find the matrix market banner in "
1059  << filename);
1060  int s2 = mm_is_coordinate(matcode), s3 = mm_is_matrix(matcode);
1061  GMM_ASSERT1(s2 > 0 && s3 > 0,
1062  "file is not coordinate storage or is not a matrix");
1063  int s4 = mm_is_pattern(matcode);
1064  GMM_ASSERT1(s4 == 0,
1065  "the file does only contain the pattern of a sparse matrix");
1066  int s5 = mm_is_skew(matcode);
1067  GMM_ASSERT1(s5 == 0, "not currently supporting skew symmetric");
1068  isSymmetric = mm_is_symmetric(matcode) || mm_is_hermitian(matcode);
1069  isHermitian = mm_is_hermitian(matcode);
1070  isComplex = mm_is_complex(matcode);
1071  mm_read_mtx_crd_size(f, &row, &col, &nz);
1072  }
1073 
1074  template <typename Matrix> void MatrixMarket_IO::read(Matrix &A) {
1075  gmm::standard_locale sl;
1076  typedef typename linalg_traits<Matrix>::value_type T;
1077  GMM_ASSERT1(f, "no file opened!");
1078  GMM_ASSERT1(!is_complex_double__(T()) || isComplex,
1079  "Bad MM matrix format (complex matrix expected)");
1080  GMM_ASSERT1(is_complex_double__(T()) || !isComplex,
1081  "Bad MM matrix format (real matrix expected)");
1082  A = Matrix(row, col);
1083  gmm::clear(A);
1084 
1085  std::vector<int> II(nz), J(nz);
1086  std::vector<typename Matrix::value_type> PR(nz);
1087  mm_read_mtx_crd_data(f, row, col, nz, &II[0], &J[0],
1088  (double*)&PR[0], matcode);
1089 
1090  for (size_type i = 0; i < size_type(nz); ++i) {
1091  A(II[i]-1, J[i]-1) = PR[i];
1092 
1093  // FIXED MM Format
1094  if (mm_is_hermitian(matcode) && (II[i] != J[i]) ) {
1095  A(J[i]-1, II[i]-1) = gmm::conj(PR[i]);
1096  }
1097 
1098  if (mm_is_symmetric(matcode) && (II[i] != J[i]) ) {
1099  A(J[i]-1, II[i]-1) = PR[i];
1100  }
1101 
1102  if (mm_is_skew(matcode) && (II[i] != J[i]) ) {
1103  A(J[i]-1, II[i]-1) = -PR[i];
1104  }
1105  }
1106  }
1107 
1108  template <typename T, typename IND_TYPE, int shift> void
1109  MatrixMarket_IO::write(const char *filename, const csc_matrix<T, IND_TYPE, shift>& A) {
1110  write(filename, csc_matrix_ref<const T*, const unsigned*,
1111  const unsigned*,shift>
1112  (&A.pr[0], &A.ir[0], &A.jc[0], A.nr, A.nc));
1113  }
1114 
1115  template <typename T, typename INDI, typename INDJ, int shift> void
1116  MatrixMarket_IO::write(const char *filename,
1117  const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A) {
1118  gmm::standard_locale sl;
1119  static MM_typecode t1 = {'M', 'C', 'R', 'G'};
1120  static MM_typecode t2 = {'M', 'C', 'C', 'G'};
1121  MM_typecode t;
1122 
1123  if (is_complex_double__(T())) std::copy(&(t2[0]), &(t2[0])+4, &(t[0]));
1124  else std::copy(&(t1[0]), &(t1[0])+4, &(t[0]));
1125  size_type nz = A.jc[mat_ncols(A)];
1126  std::vector<int> II(nz), J(nz);
1127  for (size_type j=0; j < mat_ncols(A); ++j) {
1128  for (size_type i = A.jc[j]; i < A.jc[j+1]; ++i) {
1129  II[i] = A.ir[i] + 1 - shift;
1130  J[i] = int(j + 1);
1131  }
1132  }
1133  mm_write_mtx_crd(filename, int(mat_nrows(A)), int(mat_ncols(A)),
1134  int(nz), &II[0], &J[0], (const double *)A.pr, t);
1135  }
1136 
1137 
1138  template <typename MAT> void
1139  MatrixMarket_IO::write(const char *filename, const MAT& A) {
1140  gmm::csc_matrix<typename gmm::linalg_traits<MAT>::value_type>
1141  tmp(gmm::mat_nrows(A), gmm::mat_ncols(A));
1142  gmm::copy(A,tmp);
1143  MatrixMarket_IO::write(filename, tmp);
1144  }
1145 
1146  template<typename VEC> static void vecsave(std::string fname, const VEC& V,
1147  bool binary=false, std::string Vformat="") {
1148  if (binary) {
1149  std::ofstream f(fname.c_str(), std::ofstream::binary);
1150  for (size_type i=0; i < gmm::vect_size(V); ++i)
1151  f.write(reinterpret_cast<const char*>(&V[i]), sizeof(V[i]));
1152  }
1153  else {
1154  if (Vformat.empty()){
1155  std::ofstream f(fname.c_str()); f.imbue(std::locale("C"));
1156  f.precision(16);
1157  for (size_type i=0; i < gmm::vect_size(V); ++i) f << V[i] << "\n";
1158  }
1159  else {
1160  FILE* f = fopen(fname.c_str(), "w");
1161  for (size_type i=0; i < gmm::vect_size(V); ++i) fprintf(f, Vformat.c_str(), V[i]);
1162  fclose(f);
1163  }
1164  }
1165  }
1166 
1167  template<typename VEC> static void vecload(std::string fname, const VEC& V_,
1168  bool binary=false) {
1169  VEC &V(const_cast<VEC&>(V_));
1170  if (binary) {
1171  std::ifstream f(fname.c_str(), std::ifstream::binary);
1172  for (size_type i=0; i < gmm::vect_size(V); ++i)
1173  f.read(reinterpret_cast<char*>(&V[i]), sizeof(V[i]));
1174  }
1175  else {
1176  std::ifstream f(fname.c_str()); f.imbue(std::locale("C"));
1177  for (size_type i=0; i < gmm::vect_size(V); ++i) f >> V[i];
1178  }
1179  }
1180 }
1181 
1182 
1183 #endif // GMM_INOUTPUT_H
gmm::MatrixMarket_IO
matrix input/output for MatrixMarket storage
Definition: gmm_inoutput.h:1002
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
gmm::Harwell_Boeing_save
void Harwell_Boeing_save(const std::string &filename, const csc_matrix< T, IND_TYPE, shift > &A)
save a "double" or "std::complex<double>" csc matrix into a HarwellBoeing file
Definition: gmm_inoutput.h:596
gmm::MatrixMarket_save
void MatrixMarket_save(const char *filename, const csc_matrix< T, IND_TYPE, shift > &A)
write a matrix-market file
Definition: gmm_inoutput.h:1041
gmm::Harwell_Boeing_load
void Harwell_Boeing_load(const std::string &filename, csc_matrix< T, IND_TYPE, shift > &A)
load a "double" or "std::complex<double>" csc matrix from a HarwellBoeing file
Definition: gmm_inoutput.h:634
gmm::MatrixMarket_load
void MatrixMarket_load(const char *filename, Matrix &A)
load a matrix-market file
Definition: gmm_inoutput.h:1035
gmm_kernel.h
Include the base gmm files.
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
gmm::HarwellBoeing_IO
matrix input/output for Harwell-Boeing format
Definition: gmm_inoutput.h:113
gmm::HarwellBoeing_IO::open
void open(const char *filename)
open filename and reads header
Definition: gmm_inoutput.h:271
gmm::HarwellBoeing_IO::read
void read(csc_matrix< T, IND_TYPE, shift > &A)
read the opened file
Definition: gmm_inoutput.h:327