Skip to content
Snippets Groups Projects
Select Git revision
  • 89bbb8639eead7ad8f342c5f5ae5bb5f79a2b909
  • master default protected
  • hierarchical-basis
  • alphashapes
  • bl
  • relaying
  • new_export_boris
  • oras_vs_osm
  • reassign_partitions
  • distributed_fwi
  • rename-classes
  • fix/fortran-api-example-t4
  • robust_partitions
  • reducing_files
  • fix_overlaps
  • 3115-issue-fix
  • 3023-Fillet2D-Update
  • convert_fdivs
  • tmp_jcjc24
  • fixedMeshIF
  • save_edges
  • gmsh_4_14_0
  • gmsh_4_13_1
  • gmsh_4_13_0
  • gmsh_4_12_2
  • gmsh_4_12_1
  • gmsh_4_12_0
  • gmsh_4_11_1
  • gmsh_4_11_0
  • gmsh_4_10_5
  • gmsh_4_10_4
  • gmsh_4_10_3
  • gmsh_4_10_2
  • gmsh_4_10_1
  • gmsh_4_10_0
  • gmsh_4_9_5
  • gmsh_4_9_4
  • gmsh_4_9_3
  • gmsh_4_9_2
  • gmsh_4_9_1
  • gmsh_4_9_0
41 results

GModelIO_KEY.cpp

Blame
  • gmm_inoutput.h 39.33 KiB
    // -*- c++ -*- (enables emacs c++ mode)
    //===========================================================================
    //
    // Copyright (C) 2003-2008 Yves Renard
    //
    // 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 2.1 of the License,  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 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 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 General Public License.  This exception does not however
    // invalidate any other reasons why the executable file might be covered by
    // the GNU 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) { 
          if(fgets(buf, BUFSIZ, f) == NULL) return buf; ++lcount;
          int s = SECURE_NONCHAR_SSCANF(buf,"%*s");
          GMM_ASSERT1(s >= 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 = strlen(s);
    	for (int j=last+1;j>=0;j--) {
    	  s[j] = s[j-1];
    	  if ( s[j] == '+' || s[j] == '-' ) {
    	    s[j-1] = 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;           /* Indicates 'E','D', or 'F' float format */
          char line[BUFSIZ];
    
          /*  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);
          }
    
          cout << "Valwidth = " << Valwidth << endl; getchar();
        
          /*  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] = 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)");
        if (A.pr) { delete[] A.pr; delete[] A.ir; delete[] A.jc; }
        A.nc = ncols(); A.nr = nrows();
        A.pr = 0;
        A.jc = new IND_TYPE[ncols()+1];
        A.ir = new IND_TYPE[nnz()];
        A.pr = new T[nnz()];
        readHB_data(A.jc, A.ir, (double*)A.pr);
        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];
        
        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') *strchr(Valfmt,'D') = '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 (Rhsflag == 'D') *strchr(Rhsfmt,'D') = '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<T*, unsigned*, unsigned *, shift>
    	  (A.pr, A.ir, A.jc, 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<T*, unsigned*, unsigned *, shift>
    	  (A.pr, A.ir, A.jc, 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, mat_nrows(A), 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, mat_nrows(A), 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;
        /*    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=tolower(*p),p++);  /* convert to lower case */
        for (p=crd; *p!='\0'; *p=tolower(*p),p++);  
        for (p=data_type; *p!='\0'; *p=tolower(*p),p++);
        for (p=storage_scheme; *p!='\0'; *p=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 I[],
    				  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", &I[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", &I[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", &I[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 I[], 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", I[i], J[i]);
        else
          if (mm_is_real(matcode))
            for (i=0; i<nz; i++)
    	  fprintf(f, "%d %d %20.16g\n", I[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", I[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) {
        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) {
        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> I(nz), J(nz);
        std::vector<typename Matrix::value_type> PR(nz);
        mm_read_mtx_crd_data(f, row, col, nz, &I[0], &J[0],
    			 (double*)&PR[0], matcode);
        
        for (size_type i = 0; i < size_type(nz); ++i) A(I[i]-1, J[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<T*,unsigned*,unsigned*,shift>
    	  (A.pr, A.ir, A.jc, 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) {
        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> I(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) {
    	I[i] = A.ir[i] + 1 - shift;
    	J[i] = j + 1;
          }
        }
        mm_write_mtx_crd(filename, mat_nrows(A), mat_ncols(A),
    		     nz, &I[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) {
        std::ofstream f(fname.c_str()); f.precision(16);
        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_) {
        VEC &V(const_cast<VEC&>(V_));
        std::ifstream f(fname.c_str());
        for (size_type i=0; i < gmm::vect_size(V); ++i) f >> V[i]; 
      }
    }
    
    
    #endif //  GMM_INOUTPUT_H