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