librheolef-dev 5.93-2 / usr / include / rheolef / inhb.h

# ifndef _SKIT_INHB_H
# define _SKIT_INHB_H
///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef 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 General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
//
// IN HB: Matrix & Vector Input in the Harwell-Boeing format
//
// read_harwell_boeing(..) read on the specified stream and load
//	the matrix and the optional vectors(s)
//
// author: Pierre.Saramito@imag.fr
//
// date: 21 january 1997
//
// note: for a description of the Harwell-Boeing standard, see:
//
//            Duff, et al.,  ACM TOMS Vol.15, No.1, March 1989
// 
// TODO: types not supported:
//       ".S." need to be symmetrized (may be)
//       ".Z." need to be anti-symmetrized (may be)
//       "R.E" Elemented Real (could be)
//       "C.." Complex values (difficult)
// and for rhs:
//       "M.." sparse rhs
//       
//
// We summarize the format here for convenience (from SPARSKIT2/inout.f).
// 
//   a) all lines in inout are assumed to be 80 character long.
//
//   b) the file consists of a header followed by the block of the 
//      start pointers followed by the block of the
//      indices, followed by the block of the floats values and
//      finally the numerical values of the right-hand-side if a 
//      right hand side is supplied. 
//
//   c) the file starts by a header which contains four lines if no
//      right hand side is supplied and five lines otherwise.
//      * first line contains the title (72 characters long) followed by
//        the 8-character identifier (name of the matrix, called key)
//       [ A72,A8 ]
//      * second line contains the number of lines for each
//        of the following data blocks (4 of them) and the total number 
//        of lines excluding the header.
//       [5i14]
//      * the third line contains a three character string identifying
//        the type of matrices, and the number of pointers (ncol CSR, nrow CSC),
//        and indices (nrow CRS, ncol CSC) nonzero entries (nnz).
//        The type is defined by:
//	  First Character: R Real matrix C Complex matrix P Pattern only (no numerical values supplied) 
//	  Second Character: S Symmetric U Unsymmetric H Hermitian Z Skew symmetric R Rectangular 
// 	  Third Character: A Assembled E Elemental matrices (unassembled) 
//       [A3,11X,4I14]
//      * The fourth line contains the variable fortran format
//        for the following data blocks.
//       [2A16,2A20] 
//      * The fifth line is only present if right-hand-sides are 
//        supplied. It consists of three one character-strings containing
//        the storage format for the right-hand-sides 
//        ('F'= full,'M'=sparse=same as matrix), an initial guess 
//        indicator ('G' for yes), an exact solution indicator 
//        ('X' for yes), followed by the number of right-hand-sides
//        and then the number of indices. 
//       [A3,11X,2I14] 
//
//    d) The three following blocks follow the header as described above.
//
//    e) In case the right hand-side are in sparse formats then 
//       the fourth block uses the same storage format as for the matrix
//       to describe the NRHS right hand sides provided, with a column
//       being replaced by a right hand side.
//

//
// read header: return true if success
//
# include "rheolef/inhb-aux.h"
namespace rheolef { 
template <class Size>
bool
read_harwell_boeing (
	std::istream& s, 
	char title[73], 
	char key[9], 
	char type[4], 
	Size& nptr, 
	Size& nidx, 
	Size& nnz, 
	Size& neltvl, 
        char ptrfmt[17], 
	char idxfmt[17], 
	char valfmt[21], 
	char rhsfmt[21],
        Size& totcrd, 
	Size& ptrcrd, 
	Size& idxcrd, 
	Size& valcrd, 
	Size& rhscrd, 
        char rhstype[4], 
	Size& nrhs, 
	Size& nrhsix,
	const char *file_name, 
	Size& line_no)
{
	// First line: [ A72,A8 ]
	if (! (s.get (title, 73) && s.get (key, 9)) ) {
            return hb_error (line_no, file_name, "(A72, A8) format required (ERR1)");
	}
	if (!endofline (s, 10)) {
            return hb_error (line_no, file_name, "(A72, A8) format required (ERR2)");
        }
	line_no++;

	// Second line: [5i14]
	if (! (getI(s,&totcrd) && getI(s,&ptrcrd) && getI(s,&idxcrd) && getI (s,&valcrd)) ) {
	    return hb_error (line_no, file_name, "(5I14) format required");
	}
	warning_macro ("totcrd " << totcrd);
	warning_macro ("ptrcrd " << ptrcrd);
	warning_macro ("idxcrd " << idxcrd);
	warning_macro ("valcrd " << valcrd);
	line_no++;

	// rhscrd is optional, and eat line
	if (getIoptional(s,&rhscrd)) {
	    if (!endofline (s, 10)) 
	        return hb_error (2, file_name, "(5I14) format required; line too long");
        } else {
	    if (!endofline (s, 24))
	        return hb_error (2, file_name, "(5I14) format required; line too long");
	}
	line_no++;
	
	// Third line: [A3,11X,4I14]
	char X11 [12];
	if (! (s.get (type, 4) && s.get(X11,12) && getI(s,&nidx) && getI(s,&nptr) && getI(s,&nnz)) )
	    return hb_error (line_no, file_name, "(A3,11X,4I14) format required");

	warning_macro ("X11 "  << X11);
	warning_macro ("nidx " << nidx);
	warning_macro ("nptr " << nptr);
	warning_macro ("nnz "  << nnz);
	// be quiet we type...
	upcase (type);

	// neltvl is optional; if rhscrd != 0 or full rhs, takes default value else 0
	if (getIoptional(s,&neltvl)) {
	    if (!endofline (s, 10))
	        return hb_error (line_no, file_name, "(A3,11X,4I14) format required; line too long");
	} else {
	    if (!endofline (s, 24))
	        return hb_error (line_no, file_name, "(A3,11X,4I14) format required; line too long");
	}
	line_no++;
	
	// Fourth line: [2A16,2A20]
	if (! (s.get(ptrfmt, 17) && s.get(idxfmt, 17) && s.get(valfmt, 21) && s.get(rhsfmt, 21)) )
	    return hb_error (line_no, file_name, "(2A16,2A20) format required");

	if (!endofline (s, 8))
	    return hb_error (line_no, file_name, "(2A16,2A20) format required; line too long");
	line_no++;

	// Fifth line, if right-hand side: [A3,11X,2I14] 
	if (rhscrd != 0 ) { 

	    // has rhs
	    if (! (s.get(rhstype, 4) && s.get(X11,12) && getI(s,&nrhs)) )
		return hb_error (line_no, file_name, "(A3,11X,2I14) format required");
	    
	    // nrhsix is optional
	    if (getIoptional(s,&nrhsix)) {
	        if (!endofline (s, 38))
		    return hb_error (line_no, file_name, "(A3,11X,2I14) format required: line too long");
	    } else {
	        if (!endofline (s, 52))
		    return hb_error (line_no, file_name, "(A3,11X,2I14) format required: line too long");
	    }
	    line_no++;
	} else {
	    // no rhs: set default values
	    rhstype [0] = 0;
	    nrhs = 0;
	    nrhsix = 0;
	}
	return true;
}
//
// for debug purpose
//
template <class Size>
void
trace_harwell_boeing (

	std::ostream& s,
	const char title[73], 
	const char key[9], 
	const char type[4],
	Size nptr, 
	Size nidx, 
	Size nnz, 
	Size neltvl,
       	const char ptrfmt[17], 
	const char idxfmt[17], 
	const char valfmt[21], 
	const char rhsfmt [21],
       	Size totcrd, 
	Size ptrcrd, 
	Size idxcrd, 
	Size valcrd, 
	Size rhscrd,
       	const char rhstype [4], 
	Size nrhs, 
	Size nrhsix)
{
	s << "             =======================================================================\n";
	s <<  std::endl;

	s << "title("  << strlen(title) << ") = |" << title << "|\n";
	s << "key(" << strlen(key)   << ")    = |" << key << "|\n";
	s << "type(" << strlen(type)  << ")   = |" << type << "|\n";
	s <<  std::endl;

	s << "nptr=" << nptr << "\n";
	s << "nidx=" << nidx << "\n";
	s << "nnz=" << nnz << "\n";
	s << "neltvl=" << neltvl << "\n";
	s <<  std::endl;
	
	s << "ptrfmt(" << strlen(ptrfmt) << ") = |" << ptrfmt << "|\n";
	s << "idxfmt(" << strlen(idxfmt) << ") = |" << idxfmt << "|\n";
	s << "valfmt(" << strlen(valfmt) << ") = |" << valfmt << "|\n";
	s << "rhsfmt(" << strlen(rhsfmt) << ") = |" << rhsfmt << "|\n";
	s <<  std::endl;
	
	s << "totcrd=" << totcrd << "\n";
	s << "ptrcrd=" << ptrcrd << "\n";
	s << "idxcrd=" << idxcrd << "\n";
	s << "valcrd=" << valcrd << "\n";
	s << "rhscrd=" << rhscrd << "\n";
	s <<  std::endl;
	
	s << "rhstype(" << strlen(rhsfmt) << ") = |" << rhstype << "|\n";
	s << "nrhs=" << nrhs << "\n";
	s << "nrhsix=" << nrhsix << "\n";
	s <<  std::endl;
}
//
// read matrix data
//
template <class Size, class IteratorValue, class IteratorSize>
bool
read_harwell_boeing (

	std::istream& s, 

    // input
	const char title[73],
	const char key[9],
	const char type[4], 
	Size nptr,
	Size nidx,
	Size nnz,
	Size neltvl, 
        const char ptrfmt[17],
	const char idxfmt[17],
	const char valfmt[21],
	const char rhsfmt[21],
        Size totcrd,
	Size ptrcrd,
	Size idxcrd,
	Size valcrd,
	Size rhscrd, 
        const char rhstype[4],
	Size nrhs,
	Size nrhsix,
	const char *file_name,
	Size& line_no,

    // output values already allocated [nptr+1,nnz,nnz]
	IteratorValue val,
	IteratorSize  ptr,
	IteratorSize  idx)
{
	// matrix type
	if (type [0] == 'C')
	    return hb_error (3, file_name, "complex valued matrix "<<type<<" not yet supported");
	if (type [0] == 'P' && type [2] == 'A')
	    return hb_error (3, file_name, "valued elemental matrix "<<type<<" not yet supported");
	bool has_value = (type [0] != 'P');

	// in some case, lines have are more than 81 char, be carreful
	// e.g. valfmt = "(4D22.16)" => 88 char per line
	const int line_size_max = 256;
	char      read_buf [line_size_max];

	// ptrfmt = "(16I5)" => nptrperline=16, ptrwidth=5
	// idxfmt idem
    	int nptrperline, ptrwidth;
	int nidxperline, idxwidth;

	// valfmt = "(4D22.16)" => nvalperline=4, valwidth=22, valflag='D'
    	int  nvalperline, valwidth;
    	int  nrhsperline, rhswidth;
    	bool val_has_D,   rhs_has_D;

        // parse the formats
    	if (! parse_ifmt (ptrfmt, line_no, file_name, &nptrperline, &ptrwidth))
	        return hb_error (line_no, file_name, "unexpected pointer format specification: " << ptrfmt);
    	
	if (! parse_ifmt (idxfmt, line_no, file_name, &nidxperline, &idxwidth))
	        return hb_error (line_no, file_name, "unexpected index format specification: " << idxfmt);
	
	if (valcrd && ! parse_rfmt (valfmt, line_no, file_name, &nvalperline, &valwidth, &val_has_D))
	        return hb_error (line_no, file_name, "unexpected value format specification: " << valfmt);
	
	if (rhscrd && ! parse_rfmt (rhsfmt, line_no, file_name, &nrhsperline, &rhswidth, &rhs_has_D))
	        return hb_error (line_no, file_name, "unexpected right-hand-side format specification" << rhsfmt);

	// check line pointers: try to correct it, if there is some mistakes..
	Size ptrcrd1 = (nptr+1)/nptrperline;
	if ((nptr+1) % nptrperline != 0) ptrcrd1++;
	Size idxcrd1 = nnz/nidxperline;
	if (nnz % nidxperline != 0) idxcrd1++;
	Size valcrd1;
	if (has_value) {
	    valcrd1 = nnz/nvalperline;
	    if (nnz % nvalperline != 0) valcrd1++;
	} else {
	    valcrd1 = 0;
	}
# ifdef TODO
	// TODO: vecteurs creux: le calcul est incorrect
	Size rhscrd1 = nrhs*nidx/4; 
	if ( nrhs*nidx%4 != 0) rhscrd1++;
	Size totcrd1 = ptrcrd1+idxcrd1+valcrd1+rhscrd1;
# endif /* TODO */

	if (ptrcrd1 != ptrcrd) {
	  hb_warning(2,file_name,"pointer size #2 = "<<ptrcrd<<", but "<<ptrcrd1<<" was expected (fixed)");
	  ptrcrd = ptrcrd1;
	}
	if (idxcrd1 != idxcrd) {
	  hb_warning(2,file_name,"index size #3 = "<<idxcrd<<", but "<<idxcrd1<<" was expected (fixed)");
	  idxcrd = idxcrd1;
	}
	if (valcrd1 != valcrd) {
	  hb_warning(2,file_name,"value size #4 = "<<valcrd<<", but "<<valcrd1<<" was expected (fixed)");
	  valcrd = valcrd1;
	}
# ifdef TODO
	// TODO: vecteurs creux: la verif est incorrecte
	if (rhscrd1 != rhscrd) {
	  hb_warning(2,file_name,"line-pointer #5 = "<<rhscrd<<", but "<<rhscrd1<<" was expected (fixed)");
	  rhscrd = rhscrd1;
	}
	if (totcrd1 != totcrd) {
	  hb_warning(2,file_name,"line-pointer #1 = "<<totcrd<<", but "<<totcrd1<<" was expected (fixed)");
	  totcrd = totcrd1;
	}
# endif /* TODO */
	// setw(n) does not work on input flot: use C format
	char cvalfmt[20];
	sprintf (cvalfmt, "%%lf");

	// read pointer array
    	Size iptr = 0;
        for (Size iline = 0; iline < ptrcrd; iline++) {

	    for (int iptrperline = 0; iptrperline < nptrperline; iptrperline++) {

		if (iptr == nptr+1) break;
	        s.get (read_buf, ptrwidth+1);
		int tmp;
		sscanf (read_buf, "%d", &tmp);
          	ptr[iptr] = tmp;
		ptr[iptr]--; // c_offset style
		iptr++;
       	    }
	    if (!skip_eoln (s, file_name, line_no)) return false;
	    line_no++;
    	}
	// read index array
	Size iidx = 0;
    	for (Size iline = 0; iline < idxcrd; iline++) {

	    for (int iidxperline = 0; iidxperline < nidxperline; iidxperline++) {

		if (iidx == nnz) break;
	        s.get (read_buf, idxwidth+1);
		int tmp;
		sscanf (read_buf, "%d", &tmp);
          	idx[iidx] = tmp;
		idx[iidx]--; // c_offset style
		iidx++;
            }
	    if (!skip_eoln (s, file_name, line_no)) return false;
	    line_no++;
    	}
	// read values array
	if (!has_value) return true;
	Size ival = 0;
    	if (val_has_D) {
	    
	    for (Size iline = 0; iline < valcrd; iline++) {

	        for (int ivalperline = 0; ivalperline < nvalperline; ivalperline++) {

          	    if (ival == nnz) break;
	            s.get (read_buf, valwidth+1);
	    	    convertDtoE(read_buf);
          	    double tmp;
		    sscanf (read_buf, cvalfmt, &tmp);
		    val[ival] = tmp;
          	    ival++;
                }
	        if (!skip_eoln (s, file_name, line_no)) return false;
	        line_no++;
    	    }
	} else {
	    
	    for (Size iline = 0; iline < valcrd; iline++) {

	        for (int ivalperline = 0; ivalperline < nvalperline; ivalperline++) {

          	    if (ival == nnz) break;
	            s.get (read_buf, valwidth+1);
          	    double tmp;
          	    sscanf (read_buf, cvalfmt, &tmp);
		    val[ival] = tmp;
          	    ival++;
                }
	        if (!skip_eoln (s, file_name, line_no)) return false;
	        line_no++;
    	    }
	}
	return true;
}
//
// read the (ivec)-th guess, exact solution or right-hand-side 
// in the HB file,
// assuming that 0..(ivec-1) vectors has been read using this function
// and that output array has correct size `size'
//
template <class Size, class IteratorValue>
bool
read_harwell_boeing (
	std::istream& s, 

    // input
	Size ivec,
	Size size,
	const char rhsfmt[21],
	Size nrhs,
	const char *file_name,
	Size& line_no,

    // output value already allocated [nidx]
	IteratorValue u)
{
	int nrhsperline, rhswidth;
	bool rhs_has_D;
	if (! parse_rfmt (rhsfmt, line_no, file_name, &nrhsperline, &rhswidth, &rhs_has_D)) {
	    return hb_error (line_no, file_name, "invalid right-hand-side format specification " << rhsfmt);
	}
	char crhsfmt[20];
	sprintf (crhsfmt, "%%lf");

	// need to know the index of the rhs, to skip end-of-lines
    	int istart = ivec*size;

	const int line_size_max = 256;
	char      read_buf [line_size_max];

	for (Size i = 0; i < size; i++) {

    	    s.get (read_buf, rhswidth+1);
	    if (rhs_has_D) convertDtoE(read_buf);
            double tmp;
	    sscanf (read_buf, crhsfmt, &tmp);
	    u[i] = tmp;

	    if ((istart + i + 1) % nrhsperline == 0) {
		if (!skip_eoln (s, file_name, line_no)) return false;
		line_no++;
	    }
    	}
	if (ivec + 1 == nrhs && (size + istart) % nrhsperline != 0) {

	    // last vector of the vector block: eat eoln or reaches eof
	    if (endofline(s,256)) {
		line_no++;
	    }
	    line_no++;
	}
	return true;
}
}// namespace rheolef
# endif // _SKIT_INHB_H