AxioDL
/
metaforce
mirror of https://github.com/AxioDL/metaforce.git
2
0
Fork 0
Code Issues Projects Releases Wiki Activity
metaforce/gmm/gmm_inoutput.h

1177 lines
41 KiB
C
Raw Permalink Normal View History

DCLN cooking and various bug fixes
2017-10-17 05:51:53 +00:00
/* -*- c++ -*- (enables emacs c++ mode) */
/*===========================================================================
Copyright (C) 2003-2017 Yves Renard, Julien Pommier
This file is a part of GetFEM++
GetFEM++ is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation; either version 3 of the License, or
(at your option) any later version along with the GCC Runtime Library
Exception either version 3.1 or (at your option) any later version.
This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
License and GCC Runtime Library Exception for more details.
You should have received a copy of the GNU Lesser General Public License
along with this program; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
As a special exception, you may use this file as it is a part of a free
software library without restriction. Specifically, if other files
instantiate templates or use macros or inline functions from this file,
or you compile this file and link it with other files to produce an
executable, this file does not by itself cause the resulting executable
to be covered by the GNU Lesser General Public License. This exception
does not however invalidate any other reasons why the executable file
might be covered by the GNU Lesser General Public License.
===========================================================================*/
/**@file gmm_inoutput.h
@author Yves Renard <Yves.Renard@insa-lyon.fr>
@author Julien Pommier <Julien.Pommier@insa-toulouse.fr>
@date July 8, 2003.
@brief Input/output on sparse matrices
Support Harwell-Boeing and Matrix-Market formats.
*/
#ifndef GMM_INOUTPUT_H
#define GMM_INOUTPUT_H
#include <stdio.h>
#include "gmm_kernel.h"
namespace gmm {
/*************************************************************************/
/* */
/* Functions to read and write Harwell Boeing format. */
/* */
/*************************************************************************/
// Fri Aug 15 16:29:47 EDT 1997
//
// Harwell-Boeing File I/O in C
// V. 1.0
//
// National Institute of Standards and Technology, MD.
// K.A. Remington
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// NOTICE
//
// Permission to use, copy, modify, and distribute this software and
// its documentation for any purpose and without fee is hereby granted
// provided that the above copyright notice appear in all copies and
// that both the copyright notice and this permission notice appear in
// supporting documentation.
//
// Neither the Author nor the Institution (National Institute of Standards
// and Technology) make any representations about the suitability of this
// software for any purpose. This software is provided "as is" without
// expressed or implied warranty.
// ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
inline void IOHBTerminate(const char *a) { GMM_ASSERT1(false, a);}
inline bool is_complex_double__(std::complex<double>) { return true; }
inline bool is_complex_double__(double) { return false; }
inline int ParseIfmt(const char *fmt, int* perline, int* width) {
if (SECURE_NONCHAR_SSCANF(fmt, " (%dI%d)", perline, width) != 2) {
*perline = 1;
int s = SECURE_NONCHAR_SSCANF(fmt, " (I%d)", width);
GMM_ASSERT1(s == 1, "invalid HB I-format: " << fmt);
}
return *width;
}
inline int ParseRfmt(const char *fmt, int* perline, int* width,
int* prec, int* flag) {
char p;
*perline = *width = *flag = *prec = 0;
#ifdef GMM_SECURE_CRT
if (sscanf_s(fmt, " (%d%c%d.%d)", perline, &p, sizeof(char), width, prec)
< 3 || !strchr("PEDF", p))
#else
if (sscanf(fmt, " (%d%c%d.%d)", perline, &p, width, prec) < 3
|| !strchr("PEDF", p))
#endif
{
*perline = 1;
#ifdef GMM_SECURE_CRT
int s = sscanf_s(fmt, " (%c%d.%d)", &p, sizeof(char), width, prec);
#else
int s = sscanf(fmt, " (%c%d.%d)", &p, width, prec);
#endif
GMM_ASSERT1(s>=2 && strchr("PEDF",p), "invalid HB REAL format: " << fmt);
}
*flag = p;
return *width;
}
/** matrix input/output for Harwell-Boeing format */
struct HarwellBoeing_IO {
int nrows() const { return Nrow; }
int ncols() const { return Ncol; }
int nnz() const { return Nnzero; }
int is_complex() const { return Type[0] == 'C'; }
int is_symmetric() const { return Type[1] == 'S'; }
int is_hermitian() const { return Type[1] == 'H'; }
HarwellBoeing_IO() { clear(); }
HarwellBoeing_IO(const char *filename) { clear(); open(filename); }
~HarwellBoeing_IO() { close(); }
/** open filename and reads header */
void open(const char *filename);
/** read the opened file */
template <typename T, int shift> void read(csc_matrix<T, shift>& A);
template <typename MAT> void read(MAT &M) IS_DEPRECATED;
template <typename T, int shift>
static void write(const char *filename, const csc_matrix<T, shift>& A);
template <typename T, int shift>
static void write(const char *filename, const csc_matrix<T, shift>& A,
const std::vector<T> &rhs);
template <typename T, typename INDI, typename INDJ, int shift>
static void write(const char *filename,
const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A);
template <typename T, typename INDI, typename INDJ, int shift>
static void write(const char *filename,
const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A,
const std::vector<T> &rhs);
/** static method for saving the matrix */
template <typename MAT> static void write(const char *filename,
const MAT& A) IS_DEPRECATED;
private:
FILE *f;
char Title[73], Key[9], Rhstype[4], Type[4];
int Nrow, Ncol, Nnzero, Nrhs;
char Ptrfmt[17], Indfmt[17], Valfmt[21], Rhsfmt[21];
int Ptrcrd, Indcrd, Valcrd, Rhscrd;
int lcount;
void close() { if (f) fclose(f); clear(); }
void clear() {
Nrow = Ncol = Nnzero = Nrhs = 0; f = 0; lcount = 0;
memset(Type, 0, sizeof Type);
memset(Key, 0, sizeof Key);
memset(Title, 0, sizeof Title);
}
char *getline(char *buf) {
char *p = fgets(buf, BUFSIZ, f); ++lcount;
int s = SECURE_NONCHAR_SSCANF(buf,"%*s");
GMM_ASSERT1(s >= 0 && p != 0,
"blank line in HB file at line " << lcount);
return buf;
}
int substrtoi(const char *p, size_type len) {
char s[100]; len = std::min(len, sizeof s - 1);
SECURE_STRNCPY(s, 100, p, len); s[len] = 0; return atoi(s);
}
double substrtod(const char *p, size_type len, int Valflag) {
char s[100]; len = std::min(len, sizeof s - 1);
SECURE_STRNCPY(s, 100, p, len); s[len] = 0;
if ( Valflag != 'F' && !strchr(s,'E')) {
/* insert a char prefix for exp */
int last = int(strlen(s));
for (int j=last+1;j>=0;j--) {
s[j] = s[j-1];
if ( s[j] == '+' || s[j] == '-' ) {
s[j-1] = char(Valflag);
break;
}
}
}
return atof(s);
}
template <typename IND_TYPE>
int readHB_data(IND_TYPE colptr[], IND_TYPE rowind[],
double val[]) {
/***********************************************************************/
/* This function opens and reads the specified file, interpreting its */
/* contents as a sparse matrix stored in the Harwell/Boeing standard */
/* format and creating compressed column storage scheme vectors to */
/* hold the index and nonzero value information. */
/* */
/* ---------- */
/* **CAVEAT** */
/* ---------- */
/* Parsing real formats from Fortran is tricky, and this file reader */
/* does not claim to be foolproof. It has been tested for cases */
/* when the real values are printed consistently and evenly spaced on */
/* each line, with Fixed (F), and Exponential (E or D) formats. */
/* */
/* ** If the input file does not adhere to the H/B format, the ** */
/* ** results will be unpredictable. ** */
/* */
/***********************************************************************/
int i,ind,col,offset,count;
int Ptrperline, Ptrwidth, Indperline, Indwidth;
int Valperline, Valwidth, Valprec, Nentries;
int Valflag = 'D'; /* Indicates 'E','D', or 'F' float format */
char line[BUFSIZ];
gmm::standard_locale sl;
/* Parse the array input formats from Line 3 of HB file */
ParseIfmt(Ptrfmt,&Ptrperline,&Ptrwidth);
ParseIfmt(Indfmt,&Indperline,&Indwidth);
if ( Type[0] != 'P' ) { /* Skip if pattern only */
ParseRfmt(Valfmt,&Valperline,&Valwidth,&Valprec,&Valflag);
}
/* Read column pointer array: */
offset = 0; /* if base 0 storage is declared (via macro def), */
/* then storage entries are offset by 1 */
for (count = 0, i=0;i<Ptrcrd;i++) {
getline(line);
for (col = 0, ind = 0;ind<Ptrperline;ind++) {
if (count > Ncol) break;
colptr[count] = substrtoi(line+col,Ptrwidth)-offset;
count++; col += Ptrwidth;
}
}
/* Read row index array: */
for (count = 0, i=0;i<Indcrd;i++) {
getline(line);
for (col = 0, ind = 0;ind<Indperline;ind++) {
if (count == Nnzero) break;
rowind[count] = substrtoi(line+col,Indwidth)-offset;
count++; col += Indwidth;
}
}
/* Read array of values: */
if ( Type[0] != 'P' ) { /* Skip if pattern only */
if ( Type[0] == 'C' ) Nentries = 2*Nnzero;
else Nentries = Nnzero;
count = 0;
for (i=0;i<Valcrd;i++) {
getline(line);
if (Valflag == 'D') {
// const_cast Due to aCC excentricity
char *p;
while( (p = const_cast<char *>(strchr(line,'D')) )) *p = 'E';
}
for (col = 0, ind = 0;ind<Valperline;ind++) {
if (count == Nentries) break;
val[count] = substrtod(line+col, Valwidth, Valflag);
count++; col += Valwidth;
}
}
}
return 1;
}
};
inline void HarwellBoeing_IO::open(const char *filename) {
int Totcrd,Neltvl,Nrhsix;
char line[BUFSIZ];
close();
SECURE_FOPEN(&f, filename, "r");
GMM_ASSERT1(f, "could not open " << filename);
/* First line: */
#ifdef GMM_SECURE_CRT
sscanf_s(getline(line), "%c%s", Title, 72, Key, 8);
#else
sscanf(getline(line), "%72c%8s", Title, Key);
#endif
Key[8] = Title[72] = 0;
/* Second line: */
Totcrd = Ptrcrd = Indcrd = Valcrd = Rhscrd = 0;
SECURE_NONCHAR_SSCANF(getline(line), "%d%d%d%d%d", &Totcrd, &Ptrcrd,
&Indcrd, &Valcrd, &Rhscrd);
/* Third line: */
Nrow = Ncol = Nnzero = Neltvl = 0;
#ifdef GMM_SECURE_CRT
if (sscanf_s(getline(line), "%c%d%d%d%d", Type, 3, &Nrow, &Ncol, &Nnzero,
&Neltvl) < 1)
#else
if (sscanf(getline(line), "%3c%d%d%d%d", Type, &Nrow, &Ncol, &Nnzero,
&Neltvl) < 1)
#endif
IOHBTerminate("Invalid Type info, line 3 of Harwell-Boeing file.\n");
for (size_type i = 0; i < 3; ++i) Type[i] = char(toupper(Type[i]));
/* Fourth line: */
#ifdef GMM_SECURE_CRT
if ( sscanf_s(getline(line), "%c%c%c%c",Ptrfmt, 16,Indfmt, 16,Valfmt,
20,Rhsfmt, 20) < 3)
#else
if ( sscanf(getline(line), "%16c%16c%20c%20c",Ptrfmt,Indfmt,Valfmt,
Rhsfmt) < 3)
#endif
IOHBTerminate("Invalid format info, line 4 of Harwell-Boeing file.\n");
Ptrfmt[16] = Indfmt[16] = Valfmt[20] = Rhsfmt[20] = 0;
/* (Optional) Fifth line: */
if (Rhscrd != 0 ) {
Nrhs = Nrhsix = 0;
#ifdef GMM_SECURE_CRT
if ( sscanf_s(getline(line), "%c%d%d", Rhstype, 3, &Nrhs, &Nrhsix) != 1)
#else
if ( sscanf(getline(line), "%3c%d%d", Rhstype, &Nrhs, &Nrhsix) != 1)
#endif
IOHBTerminate("Invalid RHS type information, line 5 of"
" Harwell-Boeing file.\n");
}
}
/* only valid for double and complex<double> csc matrices */
template <typename T, int shift> void
HarwellBoeing_IO::read(csc_matrix<T, shift>& A) {
// typedef typename csc_matrix<T, shift>::IND_TYPE IND_TYPE;
GMM_ASSERT1(f, "no file opened!");
GMM_ASSERT1(Type[0] != 'P',
"Bad HB matrix format (pattern matrices not supported)");
GMM_ASSERT1(!is_complex_double__(T()) || Type[0] != 'R',
"Bad HB matrix format (file contains a REAL matrix)");
GMM_ASSERT1(is_complex_double__(T()) || Type[0] != 'C',
"Bad HB matrix format (file contains a COMPLEX matrix)");
A.nc = ncols(); A.nr = nrows();
A.jc.resize(ncols()+1);
A.ir.resize(nnz());
A.pr.resize(nnz());
readHB_data(&A.jc[0], &A.ir[0], (double*)&A.pr[0]);
for (int i = 0; i <= ncols(); ++i) { A.jc[i] += shift; A.jc[i] -= 1; }
for (int i = 0; i < nnz(); ++i) { A.ir[i] += shift; A.ir[i] -= 1; }
}
template <typename MAT> void
HarwellBoeing_IO::read(MAT &M) {
csc_matrix<typename gmm::linalg_traits<MAT>::value_type> csc;
read(csc);
resize(M, mat_nrows(csc), mat_ncols(csc));
copy(csc, M);
}
template <typename IND_TYPE>
inline int writeHB_mat_double(const char* filename, int M, int N, int nz,
const IND_TYPE colptr[],
const IND_TYPE rowind[],
const double val[], int Nrhs,
const double rhs[], const double guess[],
const double exact[], const char* Title,
const char* Key, const char* Type,
const char* Ptrfmt, const char* Indfmt,
const char* Valfmt, const char* Rhsfmt,
const char* Rhstype, int shift) {
/************************************************************************/
/* The writeHB function opens the named file and writes the specified */
/* matrix and optional right-hand-side(s) to that file in */
/* Harwell-Boeing format. */
/* */
/* For a description of the Harwell Boeing standard, see: */
/* Duff, et al., ACM TOMS Vol.15, No.1, March 1989 */
/* */
/************************************************************************/
FILE *out_file;
int i, entry, offset, j, acount, linemod;
int totcrd, ptrcrd, indcrd, valcrd, rhscrd;
int nvalentries, nrhsentries;
int Ptrperline, Ptrwidth, Indperline, Indwidth;
int Rhsperline, Rhswidth, Rhsprec, Rhsflag;
int Valperline, Valwidth, Valprec;
int Valflag; /* Indicates 'E','D', or 'F' float format */
char pformat[16],iformat[16],vformat[19],rformat[19];
// char *pValflag, *pRhsflag;
gmm::standard_locale sl;
if ( Type[0] == 'C' )
{ nvalentries = 2*nz; nrhsentries = 2*M; }
else
{ nvalentries = nz; nrhsentries = M; }
if ( filename != NULL ) {
SECURE_FOPEN(&out_file, filename, "w");
GMM_ASSERT1(out_file != NULL, "Error: Cannot open file: " << filename);
} else out_file = stdout;
if ( Ptrfmt == NULL ) Ptrfmt = "(8I10)";
ParseIfmt(Ptrfmt, &Ptrperline, &Ptrwidth);
SECURE_SPRINTF1(pformat,sizeof(pformat),"%%%dd",Ptrwidth);
ptrcrd = (N+1)/Ptrperline;
if ( (N+1)%Ptrperline != 0) ptrcrd++;
if ( Indfmt == NULL ) Indfmt = Ptrfmt;
ParseIfmt(Indfmt, &Indperline, &Indwidth);
SECURE_SPRINTF1(iformat,sizeof(iformat), "%%%dd",Indwidth);
indcrd = nz/Indperline;
if ( nz%Indperline != 0) indcrd++;
if ( Type[0] != 'P' ) { /* Skip if pattern only */
if ( Valfmt == NULL ) Valfmt = "(4E21.13)";
ParseRfmt(Valfmt, &Valperline, &Valwidth, &Valprec, &Valflag);
// if (Valflag == 'D') {
// pValflag = (char *) strchr(Valfmt,'D');
// *pValflag = 'E';
// }
if (Valflag == 'F')
SECURE_SPRINTF2(vformat, sizeof(vformat), "%% %d.%df", Valwidth,
Valprec);
else
SECURE_SPRINTF2(vformat, sizeof(vformat), "%% %d.%dE", Valwidth,
Valprec);
valcrd = nvalentries/Valperline;
if ( nvalentries%Valperline != 0) valcrd++;
} else valcrd = 0;
if ( Nrhs > 0 ) {
if ( Rhsfmt == NULL ) Rhsfmt = Valfmt;
ParseRfmt(Rhsfmt,&Rhsperline,&Rhswidth,&Rhsprec, &Rhsflag);
if (Rhsflag == 'F')
SECURE_SPRINTF2(rformat,sizeof(rformat), "%% %d.%df",Rhswidth,Rhsprec);
else
SECURE_SPRINTF2(rformat,sizeof(rformat), "%% %d.%dE",Rhswidth,Rhsprec);
// if (Valflag == 'D') {
// pRhsflag = (char *) strchr(Rhsfmt,'D');
// *pRhsflag = 'E';
// }
rhscrd = nrhsentries/Rhsperline;
if ( nrhsentries%Rhsperline != 0) rhscrd++;
if ( Rhstype[1] == 'G' ) rhscrd+=rhscrd;
if ( Rhstype[2] == 'X' ) rhscrd+=rhscrd;
rhscrd*=Nrhs;
} else rhscrd = 0;
totcrd = 4+ptrcrd+indcrd+valcrd+rhscrd;
/* Print header information: */
fprintf(out_file,"%-72s%-8s\n%14d%14d%14d%14d%14d\n",Title, Key, totcrd,
ptrcrd, indcrd, valcrd, rhscrd);
fprintf(out_file,"%3s%11s%14d%14d%14d%14d\n",Type," ", M, N, nz, 0);
fprintf(out_file,"%-16s%-16s%-20s", Ptrfmt, Indfmt, Valfmt);
if ( Nrhs != 0 ) {
/* Print Rhsfmt on fourth line and */
/* optional fifth header line for auxillary vector information:*/
fprintf(out_file,"%-20s\n%-14s%d\n",Rhsfmt,Rhstype,Nrhs);
}
else
fprintf(out_file,"\n");
offset = 1 - shift; /* if base 0 storage is declared (via macro def), */
/* then storage entries are offset by 1 */
/* Print column pointers: */
for (i = 0; i < N+1; i++) {
entry = colptr[i]+offset;
fprintf(out_file,pformat,entry);
if ( (i+1)%Ptrperline == 0 ) fprintf(out_file,"\n");
}
if ( (N+1) % Ptrperline != 0 ) fprintf(out_file,"\n");
/* Print row indices: */
for (i=0;i<nz;i++) {
entry = rowind[i]+offset;
fprintf(out_file,iformat,entry);
if ( (i+1)%Indperline == 0 ) fprintf(out_file,"\n");
}
if ( nz % Indperline != 0 ) fprintf(out_file,"\n");
/* Print values: */
if ( Type[0] != 'P' ) { /* Skip if pattern only */
for (i=0;i<nvalentries;i++) {
fprintf(out_file,vformat,val[i]);
if ( (i+1)%Valperline == 0 ) fprintf(out_file,"\n");
}
if ( nvalentries % Valperline != 0 ) fprintf(out_file,"\n");
/* Print right hand sides: */
acount = 1;
linemod=0;
if ( Nrhs > 0 ) {
for (j=0;j<Nrhs;j++) {
for (i=0;i<nrhsentries;i++) {
fprintf(out_file,rformat,rhs[i] /* *Rhswidth */);
if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
}
if ( acount%Rhsperline != linemod ) {
fprintf(out_file,"\n");
linemod = (acount-1)%Rhsperline;
}
if ( Rhstype[1] == 'G' ) {
for (i=0;i<nrhsentries;i++) {
fprintf(out_file,rformat,guess[i] /* *Rhswidth */);
if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
}
if ( acount%Rhsperline != linemod ) {
fprintf(out_file,"\n");
linemod = (acount-1)%Rhsperline;
}
}
if ( Rhstype[2] == 'X' ) {
for (i=0;i<nrhsentries;i++) {
fprintf(out_file,rformat,exact[i] /* *Rhswidth */);
if ( acount++%Rhsperline == linemod ) fprintf(out_file,"\n");
}
if ( acount%Rhsperline != linemod ) {
fprintf(out_file,"\n");
linemod = (acount-1)%Rhsperline;
}
}
}
}
}
int s = fclose(out_file);
GMM_ASSERT1(s == 0, "Error closing file in writeHB_mat_double().");
return 1;
}
template <typename T, int shift> void
HarwellBoeing_IO::write(const char *filename,
const csc_matrix<T, shift>& A) {
write(filename, csc_matrix_ref<const T*, const unsigned*,
const unsigned *, shift>
(&A.pr[0], &A.ir[0], &A.jc[0], A.nr, A.nc));
}
template <typename T, int shift> void
HarwellBoeing_IO::write(const char *filename,
const csc_matrix<T, shift>& A,
const std::vector<T> &rhs) {
write(filename, csc_matrix_ref<const T*, const unsigned*,
const unsigned *, shift>
(&A.pr[0], &A.ir[0], &A.jc[0], A.nr, A.nc), rhs);
}
template <typename T, typename INDI, typename INDJ, int shift> void
HarwellBoeing_IO::write(const char *filename,
const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A) {
const char *t = 0;
if (is_complex_double__(T()))
if (mat_nrows(A) == mat_ncols(A)) t = "CUA"; else t = "CRA";
else
if (mat_nrows(A) == mat_ncols(A)) t = "RUA"; else t = "RRA";
writeHB_mat_double(filename, int(mat_nrows(A)), int(mat_ncols(A)),
A.jc[mat_ncols(A)], A.jc, A.ir,
(const double *)A.pr,
0, 0, 0, 0, "GETFEM++ CSC MATRIX", "CSCMAT",
t, 0, 0, 0, 0, "F", shift);
}
template <typename T, typename INDI, typename INDJ, int shift> void
HarwellBoeing_IO::write(const char *filename,
const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A,
const std::vector<T> &rhs) {
const char *t = 0;
if (is_complex_double__(T()))
if (mat_nrows(A) == mat_ncols(A)) t = "CUA"; else t = "CRA";
else
if (mat_nrows(A) == mat_ncols(A)) t = "RUA"; else t = "RRA";
int Nrhs = gmm::vect_size(rhs) / mat_nrows(A);
writeHB_mat_double(filename, int(mat_nrows(A)), int(mat_ncols(A)),
A.jc[mat_ncols(A)], A.jc, A.ir,
(const double *)A.pr,
Nrhs, (const double *)(&rhs[0]), 0, 0,
"GETFEM++ CSC MATRIX", "CSCMAT",
t, 0, 0, 0, 0, "F ", shift);
}
template <typename MAT> void
HarwellBoeing_IO::write(const char *filename, const MAT& A) {
gmm::csc_matrix<typename gmm::linalg_traits<MAT>::value_type>
tmp(gmm::mat_nrows(A), gmm::mat_ncols(A));
gmm::copy(A,tmp);
HarwellBoeing_IO::write(filename, tmp);
}
/** save a "double" or "std::complex<double>" csc matrix into a
HarwellBoeing file
*/
template <typename T, int shift> inline void
Harwell_Boeing_save(const std::string &filename,
const csc_matrix<T, shift>& A)
{ HarwellBoeing_IO::write(filename.c_str(), A); }
/** save a reference on "double" or "std::complex<double>" csc matrix
into a HarwellBoeing file
*/
template <typename T, typename INDI, typename INDJ, int shift> inline void
Harwell_Boeing_save(const std::string &filename,
const csc_matrix_ref<T, INDI, INDJ, shift>& A)
{ HarwellBoeing_IO::write(filename.c_str(), A); }
/** save a "double" or "std::complex<double>" generic matrix
into a HarwellBoeing file making a copy in a csc matrix
*/
template <typename MAT> inline void
Harwell_Boeing_save(const std::string &filename, const MAT& A) {
gmm::csc_matrix<typename gmm::linalg_traits<MAT>::value_type>
tmp(gmm::mat_nrows(A), gmm::mat_ncols(A));
gmm::copy(A, tmp);
HarwellBoeing_IO::write(filename.c_str(), tmp);
}
template <typename MAT, typename VECT> inline void
Harwell_Boeing_save(const std::string &filename, const MAT& A,
const VECT &RHS) {
typedef typename gmm::linalg_traits<MAT>::value_type T;
gmm::csc_matrix<T> tmp(gmm::mat_nrows(A), gmm::mat_ncols(A));
gmm::copy(A, tmp);
std::vector<T> tmprhs(gmm::vect_size(RHS));
gmm::copy(RHS, tmprhs);
HarwellBoeing_IO::write(filename.c_str(), tmp, tmprhs);
}
/** load a "double" or "std::complex<double>" csc matrix from a
HarwellBoeing file
*/
template <typename T, int shift> void
Harwell_Boeing_load(const std::string &filename, csc_matrix<T, shift>& A) {
HarwellBoeing_IO h(filename.c_str()); h.read(A);
}
/** load a "double" or "std::complex<double>" generic matrix from a
HarwellBoeing file
*/
template <typename MAT> void
Harwell_Boeing_load(const std::string &filename, MAT& A) {
csc_matrix<typename gmm::linalg_traits<MAT>::value_type> csc;
Harwell_Boeing_load(filename, csc);
resize(A, mat_nrows(csc), mat_ncols(csc));
copy(csc, A);
}
/*************************************************************************/
/* */
/* Functions to read and write MatrixMarket format. */
/* */
/*************************************************************************/
/*
* Matrix Market I/O library for ANSI C
*
* See http://math.nist.gov/MatrixMarket for details.
*
*
*/
#define MM_MAX_LINE_LENGTH 1025
#define MatrixMarketBanner "%%MatrixMarket"
#define MM_MAX_TOKEN_LENGTH 64
typedef char MM_typecode[4];
/******************* MM_typecode query functions *************************/
#define mm_is_matrix(typecode) ((typecode)[0]=='M')
#define mm_is_sparse(typecode) ((typecode)[1]=='C')
#define mm_is_coordinate(typecode) ((typecode)[1]=='C')
#define mm_is_dense(typecode) ((typecode)[1]=='A')
#define mm_is_array(typecode) ((typecode)[1]=='A')
#define mm_is_complex(typecode) ((typecode)[2]=='C')
#define mm_is_real(typecode) ((typecode)[2]=='R')
#define mm_is_pattern(typecode) ((typecode)[2]=='P')
#define mm_is_integer(typecode) ((typecode)[2]=='I')
#define mm_is_symmetric(typecode) ((typecode)[3]=='S')
#define mm_is_general(typecode) ((typecode)[3]=='G')
#define mm_is_skew(typecode) ((typecode)[3]=='K')
#define mm_is_hermitian(typecode) ((typecode)[3]=='H')
/******************* MM_typecode modify fucntions ************************/
#define mm_set_matrix(typecode) ((*typecode)[0]='M')
#define mm_set_coordinate(typecode) ((*typecode)[1]='C')
#define mm_set_array(typecode) ((*typecode)[1]='A')
#define mm_set_dense(typecode) mm_set_array(typecode)
#define mm_set_sparse(typecode) mm_set_coordinate(typecode)
#define mm_set_complex(typecode) ((*typecode)[2]='C')
#define mm_set_real(typecode) ((*typecode)[2]='R')
#define mm_set_pattern(typecode) ((*typecode)[2]='P')
#define mm_set_integer(typecode) ((*typecode)[2]='I')
#define mm_set_symmetric(typecode) ((*typecode)[3]='S')
#define mm_set_general(typecode) ((*typecode)[3]='G')
#define mm_set_skew(typecode) ((*typecode)[3]='K')
#define mm_set_hermitian(typecode) ((*typecode)[3]='H')
#define mm_clear_typecode(typecode) ((*typecode)[0]=(*typecode)[1]= \
(*typecode)[2]=' ',(*typecode)[3]='G')
#define mm_initialize_typecode(typecode) mm_clear_typecode(typecode)
/******************* Matrix Market error codes ***************************/
#define MM_COULD_NOT_READ_FILE 11
#define MM_PREMATURE_EOF 12
#define MM_NOT_MTX 13
#define MM_NO_HEADER 14
#define MM_UNSUPPORTED_TYPE 15
#define MM_LINE_TOO_LONG 16
#define MM_COULD_NOT_WRITE_FILE 17
/******************** Matrix Market internal definitions *****************
MM_matrix_typecode: 4-character sequence
object sparse/ data storage
dense type scheme
string position: [0] [1] [2] [3]
Matrix typecode: M(atrix) C(oord) R(eal) G(eneral)
A(array) C(omplex) H(ermitian)
P(attern) S(ymmetric)
I(nteger) K(kew)
***********************************************************************/
#define MM_MTX_STR "matrix"
#define MM_ARRAY_STR "array"
#define MM_DENSE_STR "array"
#define MM_COORDINATE_STR "coordinate"
#define MM_SPARSE_STR "coordinate"
#define MM_COMPLEX_STR "complex"
#define MM_REAL_STR "real"
#define MM_INT_STR "integer"
#define MM_GENERAL_STR "general"
#define MM_SYMM_STR "symmetric"
#define MM_HERM_STR "hermitian"
#define MM_SKEW_STR "skew-symmetric"
#define MM_PATTERN_STR "pattern"
inline char *mm_typecode_to_str(MM_typecode matcode) {
char buffer[MM_MAX_LINE_LENGTH];
const char *types[4] = {0,0,0,0};
/* int error =0; */
/* int i; */
/* check for MTX type */
if (mm_is_matrix(matcode))
types[0] = MM_MTX_STR;
/*
else
error=1;
*/
/* check for CRD or ARR matrix */
if (mm_is_sparse(matcode))
types[1] = MM_SPARSE_STR;
else
if (mm_is_dense(matcode))
types[1] = MM_DENSE_STR;
else
return NULL;
/* check for element data type */
if (mm_is_real(matcode))
types[2] = MM_REAL_STR;
else
if (mm_is_complex(matcode))
types[2] = MM_COMPLEX_STR;
else
if (mm_is_pattern(matcode))
types[2] = MM_PATTERN_STR;
else
if (mm_is_integer(matcode))
types[2] = MM_INT_STR;
else
return NULL;
/* check for symmetry type */
if (mm_is_general(matcode))
types[3] = MM_GENERAL_STR;
else if (mm_is_symmetric(matcode))
types[3] = MM_SYMM_STR;
else if (mm_is_hermitian(matcode))
types[3] = MM_HERM_STR;
else if (mm_is_skew(matcode))
types[3] = MM_SKEW_STR;
else
return NULL;
SECURE_SPRINTF4(buffer, sizeof(buffer), "%s %s %s %s", types[0], types[1],
types[2], types[3]);
return SECURE_STRDUP(buffer);
}
inline int mm_read_banner(FILE *f, MM_typecode *matcode) {
char line[MM_MAX_LINE_LENGTH];
char banner[MM_MAX_TOKEN_LENGTH];
char mtx[MM_MAX_TOKEN_LENGTH];
char crd[MM_MAX_TOKEN_LENGTH];
char data_type[MM_MAX_TOKEN_LENGTH];
char storage_scheme[MM_MAX_TOKEN_LENGTH];
char *p;
gmm::standard_locale sl;
/* int ret_code; */
mm_clear_typecode(matcode);
if (fgets(line, MM_MAX_LINE_LENGTH, f) == NULL)
return MM_PREMATURE_EOF;
#ifdef GMM_SECURE_CRT
if (sscanf_s(line, "%s %s %s %s %s", banner, sizeof(banner),
mtx, sizeof(mtx), crd, sizeof(crd), data_type,
sizeof(data_type), storage_scheme,
sizeof(storage_scheme)) != 5)
#else
if (sscanf(line, "%s %s %s %s %s", banner, mtx, crd,
data_type, storage_scheme) != 5)
#endif
return MM_PREMATURE_EOF;
for (p=mtx; *p!='\0'; *p=char(tolower(*p)),p++) {}; /* convert to lower case */
for (p=crd; *p!='\0'; *p=char(tolower(*p)),p++) {};
for (p=data_type; *p!='\0'; *p=char(tolower(*p)),p++) {};
for (p=storage_scheme; *p!='\0'; *p=char(tolower(*p)),p++) {};
/* check for banner */
if (strncmp(banner, MatrixMarketBanner, strlen(MatrixMarketBanner)) != 0)
return MM_NO_HEADER;
/* first field should be "mtx" */
if (strcmp(mtx, MM_MTX_STR) != 0)
return MM_UNSUPPORTED_TYPE;
mm_set_matrix(matcode);
/* second field describes whether this is a sparse matrix (in coordinate
storgae) or a dense array */
if (strcmp(crd, MM_SPARSE_STR) == 0)
mm_set_sparse(matcode);
else
if (strcmp(crd, MM_DENSE_STR) == 0)
mm_set_dense(matcode);
else
return MM_UNSUPPORTED_TYPE;
/* third field */
if (strcmp(data_type, MM_REAL_STR) == 0)
mm_set_real(matcode);
else
if (strcmp(data_type, MM_COMPLEX_STR) == 0)
mm_set_complex(matcode);
else
if (strcmp(data_type, MM_PATTERN_STR) == 0)
mm_set_pattern(matcode);
else
if (strcmp(data_type, MM_INT_STR) == 0)
mm_set_integer(matcode);
else
return MM_UNSUPPORTED_TYPE;
/* fourth field */
if (strcmp(storage_scheme, MM_GENERAL_STR) == 0)
mm_set_general(matcode);
else
if (strcmp(storage_scheme, MM_SYMM_STR) == 0)
mm_set_symmetric(matcode);
else
if (strcmp(storage_scheme, MM_HERM_STR) == 0)
mm_set_hermitian(matcode);
else
if (strcmp(storage_scheme, MM_SKEW_STR) == 0)
mm_set_skew(matcode);
else
return MM_UNSUPPORTED_TYPE;
return 0;
}
inline int mm_read_mtx_crd_size(FILE *f, int *M, int *N, int *nz ) {
char line[MM_MAX_LINE_LENGTH];
/* int ret_code;*/
int num_items_read;
/* set return null parameter values, in case we exit with errors */
*M = *N = *nz = 0;
/* now continue scanning until you reach the end-of-comments */
do {
if (fgets(line,MM_MAX_LINE_LENGTH,f) == NULL)
return MM_PREMATURE_EOF;
} while (line[0] == '%');
/* line[] is either blank or has M,N, nz */
if (SECURE_NONCHAR_SSCANF(line, "%d %d %d", M, N, nz) == 3) return 0;
else
do {
num_items_read = SECURE_NONCHAR_FSCANF(f, "%d %d %d", M, N, nz);
if (num_items_read == EOF) return MM_PREMATURE_EOF;
}
while (num_items_read != 3);
return 0;
}
inline int mm_read_mtx_crd_data(FILE *f, int, int, int nz, int II[],
int J[], double val[], MM_typecode matcode) {
int i;
if (mm_is_complex(matcode)) {
for (i=0; i<nz; i++)
if (SECURE_NONCHAR_FSCANF(f, "%d %d %lg %lg", &II[i], &J[i],
&val[2*i], &val[2*i+1])
!= 4) return MM_PREMATURE_EOF;
}
else if (mm_is_real(matcode)) {
for (i=0; i<nz; i++) {
if (SECURE_NONCHAR_FSCANF(f, "%d %d %lg\n", &II[i], &J[i], &val[i])
!= 3) return MM_PREMATURE_EOF;
}
}
else if (mm_is_pattern(matcode)) {
for (i=0; i<nz; i++)
if (SECURE_NONCHAR_FSCANF(f, "%d %d", &II[i], &J[i])
!= 2) return MM_PREMATURE_EOF;
}
else return MM_UNSUPPORTED_TYPE;
return 0;
}
inline int mm_write_mtx_crd(const char *fname, int M, int N, int nz,
int II[], int J[], const double val[],
MM_typecode matcode) {
FILE *f;
int i;
if (strcmp(fname, "stdout") == 0)
f = stdout;
else {
SECURE_FOPEN(&f, fname, "w");
if (f == NULL)
return MM_COULD_NOT_WRITE_FILE;
}
/* print banner followed by typecode */
fprintf(f, "%s ", MatrixMarketBanner);
char *str = mm_typecode_to_str(matcode);
fprintf(f, "%s\n", str);
free(str);
/* print matrix sizes and nonzeros */
fprintf(f, "%d %d %d\n", M, N, nz);
/* print values */
if (mm_is_pattern(matcode))
for (i=0; i<nz; i++)
fprintf(f, "%d %d\n", II[i], J[i]);
else
if (mm_is_real(matcode))
for (i=0; i<nz; i++)
fprintf(f, "%d %d %20.16g\n", II[i], J[i], val[i]);
else
if (mm_is_complex(matcode))
for (i=0; i<nz; i++)
fprintf(f, "%d %d %20.16g %20.16g\n", II[i], J[i], val[2*i],
val[2*i+1]);
else {
if (f != stdout) fclose(f);
return MM_UNSUPPORTED_TYPE;
}
if (f !=stdout) fclose(f);
return 0;
}
/** matrix input/output for MatrixMarket storage */
class MatrixMarket_IO {
FILE *f;
bool isComplex, isSymmetric, isHermitian;
int row, col, nz;
MM_typecode matcode;
public:
MatrixMarket_IO() : f(0) {}
MatrixMarket_IO(const char *filename) : f(0) { open(filename); }
~MatrixMarket_IO() { if (f) fclose(f); f = 0; }
int nrows() const { return row; }
int ncols() const { return col; }
int nnz() const { return nz; }
int is_complex() const { return isComplex; }
int is_symmetric() const { return isSymmetric; }
int is_hermitian() const { return isHermitian; }
/* open filename and reads header */
void open(const char *filename);
/* read opened file */
template <typename Matrix> void read(Matrix &A);
/* write a matrix */
template <typename T, int shift> static void
write(const char *filename, const csc_matrix<T, shift>& A);
template <typename T, typename INDI, typename INDJ, int shift> static void
write(const char *filename,
const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A);
template <typename MAT> static void
write(const char *filename, const MAT& A);
};
/** load a matrix-market file */
template <typename Matrix> inline void
MatrixMarket_load(const char *filename, Matrix& A) {
MatrixMarket_IO mm; mm.open(filename);
mm.read(A);
}
/** write a matrix-market file */
template <typename T, int shift> void
MatrixMarket_save(const char *filename, const csc_matrix<T, shift>& A) {
MatrixMarket_IO mm; mm.write(filename, A);
}
template <typename T, typename INDI, typename INDJ, int shift> inline void
MatrixMarket_save(const char *filename,
const csc_matrix_ref<T, INDI, INDJ, shift>& A) {
MatrixMarket_IO mm; mm.write(filename, A);
}
inline void MatrixMarket_IO::open(const char *filename) {
gmm::standard_locale sl;
if (f) { fclose(f); }
SECURE_FOPEN(&f, filename, "r");
GMM_ASSERT1(f, "Sorry, cannot open file " << filename);
int s1 = mm_read_banner(f, &matcode);
GMM_ASSERT1(s1 == 0, "Sorry, cannnot find the matrix market banner in "
<< filename);
int s2 = mm_is_coordinate(matcode), s3 = mm_is_matrix(matcode);
GMM_ASSERT1(s2 > 0 && s3 > 0,
"file is not coordinate storage or is not a matrix");
int s4 = mm_is_pattern(matcode);
GMM_ASSERT1(s4 == 0,
"the file does only contain the pattern of a sparse matrix");
int s5 = mm_is_skew(matcode);
GMM_ASSERT1(s5 == 0, "not currently supporting skew symmetric");
isSymmetric = mm_is_symmetric(matcode) || mm_is_hermitian(matcode);
isHermitian = mm_is_hermitian(matcode);
isComplex = mm_is_complex(matcode);
mm_read_mtx_crd_size(f, &row, &col, &nz);
}
template <typename Matrix> void MatrixMarket_IO::read(Matrix &A) {
gmm::standard_locale sl;
typedef typename linalg_traits<Matrix>::value_type T;
GMM_ASSERT1(f, "no file opened!");
GMM_ASSERT1(!is_complex_double__(T()) || isComplex,
"Bad MM matrix format (complex matrix expected)");
GMM_ASSERT1(is_complex_double__(T()) || !isComplex,
"Bad MM matrix format (real matrix expected)");
A = Matrix(row, col);
gmm::clear(A);
std::vector<int> II(nz), J(nz);
std::vector<typename Matrix::value_type> PR(nz);
mm_read_mtx_crd_data(f, row, col, nz, &II[0], &J[0],
(double*)&PR[0], matcode);
for (size_type i = 0; i < size_type(nz); ++i) {
A(II[i]-1, J[i]-1) = PR[i];
// FIXED MM Format
if (mm_is_hermitian(matcode) && (II[i] != J[i]) ) {
A(J[i]-1, II[i]-1) = gmm::conj(PR[i]);
}
if (mm_is_symmetric(matcode) && (II[i] != J[i]) ) {
A(J[i]-1, II[i]-1) = PR[i];
}
if (mm_is_skew(matcode) && (II[i] != J[i]) ) {
A(J[i]-1, II[i]-1) = -PR[i];
}
}
}
template <typename T, int shift> void
MatrixMarket_IO::write(const char *filename, const csc_matrix<T, shift>& A) {
write(filename, csc_matrix_ref<const T*, const unsigned*,
const unsigned*,shift>
(&A.pr[0], &A.ir[0], &A.jc[0], A.nr, A.nc));
}
template <typename T, typename INDI, typename INDJ, int shift> void
MatrixMarket_IO::write(const char *filename,
const csc_matrix_ref<T*, INDI*, INDJ*, shift>& A) {
gmm::standard_locale sl;
static MM_typecode t1 = {'M', 'C', 'R', 'G'};
static MM_typecode t2 = {'M', 'C', 'C', 'G'};
MM_typecode t;
if (is_complex_double__(T())) std::copy(&(t2[0]), &(t2[0])+4, &(t[0]));
else std::copy(&(t1[0]), &(t1[0])+4, &(t[0]));
size_type nz = A.jc[mat_ncols(A)];
std::vector<int> II(nz), J(nz);
for (size_type j=0; j < mat_ncols(A); ++j) {
for (size_type i = A.jc[j]; i < A.jc[j+1]; ++i) {
II[i] = A.ir[i] + 1 - shift;
J[i] = int(j + 1);
}
}
mm_write_mtx_crd(filename, int(mat_nrows(A)), int(mat_ncols(A)),
int(nz), &II[0], &J[0], (const double *)A.pr, t);
}
template <typename MAT> void
MatrixMarket_IO::write(const char *filename, const MAT& A) {
gmm::csc_matrix<typename gmm::linalg_traits<MAT>::value_type>
tmp(gmm::mat_nrows(A), gmm::mat_ncols(A));
gmm::copy(A,tmp);
MatrixMarket_IO::write(filename, tmp);
}
template<typename VEC> static void vecsave(std::string fname, const VEC& V,
bool binary=false) {
if (binary) {
std::ofstream f(fname.c_str(), std::ofstream::binary);
for (size_type i=0; i < gmm::vect_size(V); ++i)
f.write(reinterpret_cast<const char*>(&V[i]), sizeof(V[i]));
}
else {
std::ofstream f(fname.c_str()); f.precision(16); f.imbue(std::locale("C"));
for (size_type i=0; i < gmm::vect_size(V); ++i) f << V[i] << "\n";
}
}
template<typename VEC> static void vecload(std::string fname, const VEC& V_,
bool binary=false) {
VEC &V(const_cast<VEC&>(V_));
if (binary) {
std::ifstream f(fname.c_str(), std::ifstream::binary);
for (size_type i=0; i < gmm::vect_size(V); ++i)
f.read(reinterpret_cast<char*>(&V[i]), sizeof(V[i]));
}
else {
std::ifstream f(fname.c_str()); f.imbue(std::locale("C"));
for (size_type i=0; i < gmm::vect_size(V); ++i) f >> V[i];
}
}
}
#endif // GMM_INOUTPUT_H