dune-istl  2.2.0
matrixmarket.hh
Go to the documentation of this file.
00001 #ifndef DUNE_MATRIXMARKET_HH
00002 #define DUNE_MATRIXMARKET_HH
00003 
00004 #include<ostream>
00005 #include<istream>
00006 #include<fstream>
00007 #include<sstream>
00008 #include<limits>
00009 #include<ios>
00010 #include"matrixutils.hh"
00011 #include "bcrsmatrix.hh"
00012 #include"owneroverlapcopy.hh"
00013 #include<dune/common/fmatrix.hh>
00014 #include<dune/common/tuples.hh>
00015 #include<dune/common/misc.hh>
00016 
00017 namespace Dune
00018 {
00019 
00045   namespace
00046   {
00056     template<class T>
00057     struct mm_numeric_type{
00058       enum{
00062         is_numeric=false
00063       };
00064     };
00065 
00066     template<>
00067     struct mm_numeric_type<int>
00068     {
00069       enum{
00073         is_numeric=true
00074       };
00075       
00076       static std::string str()
00077       {
00078         return "integer";
00079       }
00080     };
00081     
00082     template<>
00083     struct mm_numeric_type<double>
00084     {
00085       enum{
00089         is_numeric=true
00090       };
00091       
00092       static std::string str()
00093       {
00094         return "real";
00095       }
00096     };
00097     
00098     template<>
00099     struct mm_numeric_type<float>
00100     {
00101       enum{
00105         is_numeric=true
00106       };
00107       
00108       static std::string str()
00109       {
00110         return "real";
00111       }
00112     };
00113 
00114     template<>
00115     struct mm_numeric_type<std::complex<double> >
00116     {
00117       enum{
00121         is_numeric=true
00122       };
00123       
00124       static std::string str()
00125       {
00126         return "complex";
00127       }
00128     };
00129     
00130     template<>
00131     struct mm_numeric_type<std::complex<float> >
00132     {
00133       enum{
00137         is_numeric=true
00138       };
00139       
00140       static std::string str()
00141       {
00142         return "complex";
00143       }
00144     };
00145 
00154     template<class M>
00155     struct mm_header_printer;
00156     
00157     template<typename T, typename A, int i, int j>
00158     struct mm_header_printer<BCRSMatrix<FieldMatrix<T,i,j>,A> >
00159     {
00160       static void print(std::ostream& os)
00161       {
00162         os<<"%%MatrixMarket matrix coordinate ";
00163         os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
00164       }
00165     };
00166 
00167     template<typename B, typename A>
00168     struct mm_header_printer<BlockVector<B,A> >
00169     {
00170       static void print(std::ostream& os)
00171       {
00172         os<<"%%MatrixMarket matrix array ";
00173         os<<mm_numeric_type<typename B::field_type>::str()<<" general"<<std::endl;
00174       }
00175     };
00176 
00177     template<typename T, int j>
00178     struct mm_header_printer<FieldVector<T,j> >
00179     {
00180       static void print(std::ostream& os)
00181       {
00182         os<<"%%MatrixMarket matrix array ";
00183         os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
00184       }
00185     };
00186 
00187     template<typename T, int i, int j>
00188     struct mm_header_printer<FieldMatrix<T,i,j> >
00189     {
00190       static void print(std::ostream& os)
00191       {
00192         os<<"%%MatrixMarket matrix array ";
00193         os<<mm_numeric_type<T>::str()<<" general"<<std::endl;
00194       }
00195     };
00196 
00205     template<class M>
00206     struct mm_block_structure_header;
00207     
00208     template<typename T, typename A, int i>
00209     struct mm_block_structure_header<BlockVector<FieldVector<T,i>,A> >
00210     {
00211       typedef BlockVector<FieldVector<T,i>,A> M;
00212       
00213       static void print(std::ostream& os, const M& m)
00214       {
00215         os<<"% ISTL_STRUCT blocked ";
00216         os<<i<<" "<<1<<std::endl;
00217       }
00218     };
00219 
00220     template<typename T, typename A, int i, int j>
00221     struct mm_block_structure_header<BCRSMatrix<FieldMatrix<T,i,j>,A> >
00222     {
00223       typedef BCRSMatrix<FieldMatrix<T,i,j>,A> M;
00224       
00225       static void print(std::ostream& os, const M& m)
00226       {
00227         os<<"% ISTL_STRUCT blocked ";
00228         os<<i<<" "<<j<<std::endl;
00229       }
00230     };
00231     
00232 
00233     template<typename T, int i, int j>
00234     struct mm_block_structure_header<FieldMatrix<T,i,j> >
00235     {
00236       typedef FieldMatrix<T,i,j> M;
00237       
00238       static void print(std::ostream& os, const M& m)
00239       {}
00240     };
00241 
00242     template<typename T, int i>
00243     struct mm_block_structure_header<FieldVector<T,i> >
00244     {
00245       typedef FieldVector<T,i> M;
00246       
00247       static void print(std::ostream& os, const M& m)
00248       {}
00249     };
00250 
00251     enum LineType{ MM_HEADER, MM_ISTLSTRUCT, DATA };
00252     enum{ MM_MAX_LINE_LENGTH=1025 };
00253     
00254     enum MM_TYPE { coordinate_type, array_type, unknown_type };
00255     
00256     enum MM_CTYPE { integer_type, double_type, complex_type, pattern, unknown_ctype };
00257     
00258     enum MM_STRUCTURE { general, symmetric, skew_symmetric, hermitian, unknown_structure  };
00259     
00260     struct MMHeader
00261     {
00262       MMHeader()
00263         : type(coordinate_type), ctype(double_type), structure(general)
00264       {}
00265       MM_TYPE type;
00266       MM_CTYPE ctype;
00267       MM_STRUCTURE structure;
00268     };
00269     
00270     bool lineFeed(std::istream& file)
00271     {
00272       char c;
00273       if(!file.eof())
00274         c=file.peek();
00275       else 
00276         return false;
00277       // ignore whitespace
00278       while(c==' ')
00279         {
00280           file.get();
00281           if(file.eof())
00282             return false;
00283           c=file.peek();
00284         }
00285       
00286       if(c=='\n'){
00287       /* eat the line feed */
00288         file.get();
00289         return true;
00290       }
00291       return false;
00292     }
00293     
00294     void skipComments(std::istream& file)
00295     {
00296       lineFeed(file);
00297       char c=file.peek();
00298       // ignore comment lines
00299       while(c=='%')
00300         {
00301           /* disgard the rest of the line */
00302           file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00303           c=file.peek();
00304         }
00305     }
00306     
00307 
00308     bool readMatrixMarketBanner(std::istream& file, MMHeader& mmHeader)
00309     {
00310       std::string buffer;
00311       char c;
00312       file >> buffer;
00313       c=buffer[0];
00314       mmHeader=MMHeader();
00315       if(c!='%')
00316         return false;
00317       std::cout<<buffer<<std::endl;
00318       /* read the banner */
00319       if(buffer!="%%MatrixMarket"){
00320         /* disgard the rest of the line */
00321         file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00322         return false;
00323       }
00324 
00325       if(lineFeed(file))
00326         /* premature end of line */
00327         return false;
00328 
00329       /* read the matrix_type */
00330       file >> buffer;
00331 
00332       if(buffer != "matrix")
00333         {
00334           /* disgard the rest of the line */
00335           file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00336           return false;
00337         }
00338 
00339       if(lineFeed(file))
00340         /* premature end of line */
00341         return false;
00342         
00343       /* The type of the matrix */
00344       file >> buffer;
00345         
00346       if(buffer.empty())
00347         return false;
00348 
00349       std::transform(buffer.begin(), buffer.end(), buffer.begin(), 
00350                      tolower);
00351 
00352       switch(buffer[0])
00353         {
00354         case 'a':
00355           /* sanity check */
00356           if(buffer != "array")
00357             {
00358               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00359               return false;
00360             }
00361           mmHeader.type=array_type;
00362           break;
00363         case 'c':
00364           /* sanity check */
00365           if(buffer != "coordinate")
00366             {
00367               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00368               return false;
00369             }
00370           mmHeader.type=coordinate_type;
00371           break;
00372         default:
00373           file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00374           return false;
00375         }
00376         
00377       if(lineFeed(file))
00378         /* premature end of line */
00379         return false;
00380         
00381       /* The numeric type used. */
00382       file >> buffer;
00383         
00384       if(buffer.empty())
00385         return false;
00386         
00387       std::transform(buffer.begin(), buffer.end(), buffer.begin(), 
00388                      tolower);
00389       switch(buffer[0])
00390         {
00391         case 'i':
00392           /* sanity check */
00393           if(buffer != "integer")
00394             {
00395               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00396               return false;
00397             }
00398           mmHeader.ctype=integer_type;
00399           break;
00400         case 'r':
00401           /* sanity check */
00402           if(buffer != "real")
00403             {
00404               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00405               return false;
00406             }
00407           mmHeader.ctype=double_type;
00408           break;
00409         case 'c':
00410           /* sanity check */
00411           if(buffer != "complex")
00412             {
00413               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00414               return false;
00415             }
00416           mmHeader.ctype=complex_type;
00417           break;
00418         case 'p':
00419           /* sanity check */
00420           if(buffer != "pattern")
00421             {
00422               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00423               return false;
00424             }
00425           mmHeader.ctype=pattern;
00426           break;
00427         default:
00428           file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00429           return false;
00430         }
00431 
00432       if(lineFeed(file))
00433         return false;
00434         
00435       file >> buffer;
00436         
00437       std::transform(buffer.begin(), buffer.end(), buffer.begin(), 
00438                      tolower);
00439       switch(buffer[0])
00440         {
00441         case 'g':
00442           /* sanity check */
00443           if(buffer != "general")
00444             {
00445               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00446               return false;
00447             }
00448           mmHeader.structure=general;
00449           break;
00450         case 'h':
00451           /* sanity check */
00452           if(buffer != "hermitian")
00453             {
00454               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00455               return false;
00456             }
00457           mmHeader.structure=hermitian;
00458           break;
00459         case 's':
00460           if(buffer.size()==1){
00461             file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00462             return false;
00463           }
00464             
00465           switch(buffer[1])
00466             {
00467             case 'y':
00468               /* sanity check */
00469               if(buffer != "symmetric")
00470                 {
00471                   file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00472                   return false;
00473                 }
00474               mmHeader.structure=symmetric;
00475               break;
00476             case 'k':
00477               /* sanity check */
00478               if(buffer != "skew-symmetric")
00479                 {
00480                   file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00481                   return false;
00482                 }
00483               mmHeader.structure=skew_symmetric;
00484               break;
00485             default:
00486               file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00487               return false;
00488             }
00489         default:
00490           file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00491           return false;
00492         }
00493       file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00494       c=file.peek();
00495       return true;
00496     
00497     }
00498 
00499     void readNextLine(std::istream& file, std::ostringstream& line, LineType& type)
00500     {
00501       char c;
00502       std::size_t index=0;
00503       
00504       //empty lines will be disgarded and we will simply read the next line
00505       while(index==0&&!file.eof())
00506         {
00507           // strip spaces
00508           while(!file.eof() && (c=file.get())==' ');
00509           
00510           //read the rest of the line until comment
00511           while(!file.eof() && (c=file.get())=='\n'){
00512             switch(c)
00513               {
00514               case '%':
00515                 /* disgard the rest of the line */
00516                 file.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00517               }
00518           }
00519         }
00520       
00521       //      buffer[index]='\0';
00522     }
00523 
00524     template<std::size_t brows, std::size_t bcols>
00525     Dune::tuple<std::size_t, std::size_t, std::size_t>
00526     calculateNNZ(std::size_t rows, std::size_t cols, std::size_t entries, const MMHeader& header)
00527     {
00528       std::size_t blockrows=rows/brows;
00529       std::size_t blockcols=cols/bcols;
00530       std::size_t blocksize=brows*bcols;
00531       std::size_t blockentries=0;
00532     
00533       switch(header.structure)
00534         {
00535         case general:
00536           blockentries = entries/blocksize; break;
00537         case skew_symmetric:
00538           blockentries = 2*entries/blocksize; break;
00539         case symmetric:
00540           blockentries = (2*entries-rows)/blocksize; break;
00541         case hermitian:
00542           blockentries = (2*entries-rows)/blocksize; break;
00543         default:
00544           throw Dune::NotImplemented();
00545         }
00546       return Dune::make_tuple(blockrows, blockcols, blockentries);
00547     }
00548 
00549     /* 
00550      *   @brief Storage class for the column index and the numeric value.
00551      *
00552      * \tparam T Either a NumericWrapper of the numeric type or PatternDummy
00553      *     for MatrixMarket pattern case.
00554      */
00555     template<typename T>
00556     struct IndexData : public T
00557     {
00558       std::size_t index;
00559     };
00560 
00561   
00572     template<typename T>
00573     struct NumericWrapper
00574     {
00575       T number;
00576       operator T&()
00577       {
00578         return number;
00579       }
00580     };
00581 
00585     struct PatternDummy
00586     {};
00587   
00588     template<>
00589     struct NumericWrapper<PatternDummy>
00590     {};
00591 
00592     template<typename T>
00593     std::istream& operator>>(std::istream& is, NumericWrapper<T>& num)
00594     {
00595       return is>>num.number;
00596     }
00597   
00598     std::istream& operator>>(std::istream& is, NumericWrapper<PatternDummy>& num)
00599     {    
00600       return is;
00601     }
00602 
00608     template<typename T>
00609     bool operator<(const IndexData<T>& i1, const IndexData<T>& i2)
00610     {
00611       return i1.index<i2.index;
00612     };
00613 
00619     template<typename T>
00620     std::istream& operator>>(std::istream& is, IndexData<T>& data)
00621     {
00622       is>>data.index;
00623       /* MatrixMarket indices are one based. Decrement for C++ */
00624       --data.index;
00625       return is>>data.number;
00626     }
00627 
00634     template<typename D, int brows, int bcols>
00635     struct MatrixValuesSetter
00636     {
00642       template<typename M>
00643       void operator()(const std::vector<std::set<IndexData<D> > >& rows,
00644                       M& matrix)
00645       {    
00646         for(typename M::RowIterator iter=matrix.begin();
00647             iter!= matrix.end(); ++iter)
00648           {
00649             for(typename M::size_type brow=iter.index()*brows, 
00650                   browend=iter.index()*brows+brows;
00651                 brow<browend; ++brow)
00652               {
00653                 typedef typename std::set<IndexData<D> >::const_iterator Siter;
00654                 for(Siter siter=rows[brow].begin(), send=rows[brow].end();
00655                     siter != send; ++siter)
00656                   (*iter)[siter->index/bcols][brow%brows][siter->index%bcols]=siter->number;
00657               }
00658           }
00659       }     
00660     };
00661   
00662     template<int brows, int bcols>
00663     struct MatrixValuesSetter<PatternDummy,brows,bcols>
00664     {
00665       template<typename M>
00666       void operator()(const std::vector<std::set<IndexData<PatternDummy> > >& rows,
00667                       M& matrix)
00668       {}
00669     };
00670   
00671     template<typename T, typename A, int brows, int bcols, typename D>
00672     void readSparseEntries(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
00673                            std::istream& file, std::size_t entries,
00674                            const MMHeader& mmHeader, const D&)
00675     {
00676       typedef Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A> Matrix;
00677       // First path
00678       // store entries together with column index in a speparate
00679       // data structure
00680       std::vector<std::set<IndexData<D> > > rows(matrix.N()*brows);
00681     
00682       for(; entries>0;--entries){
00683         std::size_t row;
00684         IndexData<D> data;
00685         skipComments(file);
00686         file>>row;
00687         --row; // Index was 1 based.
00688         assert(row/bcols<matrix.N());
00689         file>>data;
00690         assert(data.index/bcols<matrix.M());
00691         rows[row].insert(data);
00692       }
00693 
00694       // TODO extend to capture the nongeneral cases.
00695       if(mmHeader.structure!= general)
00696         DUNE_THROW(Dune::NotImplemented, "Only general is supported right now!");
00697     
00698       // Setup the matrix sparsity pattern
00699       int nnz=0;
00700       for(typename Matrix::CreateIterator iter=matrix.createbegin();
00701           iter!= matrix.createend(); ++iter)
00702         {
00703           for(std::size_t brow=iter.index()*brows, browend=iter.index()*brows+brows;
00704               brow<browend; ++brow)
00705             {
00706               typedef typename std::set<IndexData<D> >::const_iterator Siter;
00707               for(Siter siter=rows[brow].begin(), send=rows[brow].end();
00708                   siter != send; ++siter, ++nnz)
00709                 iter.insert(siter->index/bcols);
00710             }
00711         }
00712     
00713       //Set the matrix values
00714       matrix=0;
00715 
00716       MatrixValuesSetter<D,brows,bcols> Setter;
00717     
00718       Setter(rows, matrix);
00719     }
00720   } // end anonymous namespace
00721 
00722   class MatrixMarketFormatError : public Dune::Exception
00723   {};
00724 
00725 
00726   void mm_read_header(std::size_t& rows, std::size_t& cols, MMHeader& header, std::istream& istr,
00727                       bool isVector)
00728   {
00729     if(!readMatrixMarketBanner(istr, header)){
00730       std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
00731                 << "%%MatrixMarket matrix coordinate real general"<<std::endl;
00732       // Go to the beginning of the file
00733       istr.clear() ;
00734       istr.seekg(0, std::ios::beg);
00735       if(isVector)
00736         header.type=array_type;
00737     }
00738   
00739     skipComments(istr);
00740 
00741     if(lineFeed(istr))
00742       throw MatrixMarketFormatError();
00743   
00744     istr >> rows;
00745   
00746     if(lineFeed(istr))
00747       throw MatrixMarketFormatError();
00748     istr >> cols;
00749   }
00750 
00751   template<typename T, typename A, int entries>
00752   void mm_read_vector_entries(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
00753                               std::size_t size,
00754                               std::istream& istr)
00755   {
00756     for(int i=0;size>0; ++i, --size){
00757       T val;
00758       istr>>val;
00759       vector[i/entries][i%entries]=val;
00760     }
00761   }
00762   
00763 
00770   template<typename T, typename A, int entries>
00771   void readMatrixMarket(Dune::BlockVector<Dune::FieldVector<T,entries>,A>& vector,
00772                         std::istream& istr)
00773   {
00774     MMHeader header;
00775     std::size_t rows, cols;
00776     mm_read_header(rows,cols,header,istr, true);
00777     if(cols!=1)
00778       DUNE_THROW(MatrixMarketFormatError, "cols!=1, therefore this is no vector!");
00779     
00780     if(header.type!=array_type)
00781       DUNE_THROW(MatrixMarketFormatError, "Vectors have to be stored in array format!");
00782     
00783     std::size_t size=rows/entries;
00784     if(size*entries!=rows)
00785       DUNE_THROW(MatrixMarketFormatError, "Block size of vector is not correct1");
00786     
00787     vector.resize(size);
00788       
00789     istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00790 
00791     mm_read_vector_entries(vector, rows, istr);
00792   }
00793   
00794   
00801   template<typename T, typename A, int brows, int bcols>
00802   void readMatrixMarket(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>& matrix,
00803                         std::istream& istr)
00804   {
00805 
00806     typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,brows,bcols> > Matrix;
00807   
00808     MMHeader header;
00809     if(!readMatrixMarketBanner(istr, header)){
00810       std::cerr << "First line was not a correct Matrix Market banner. Using default:\n"
00811                 << "%%MatrixMarket matrix coordinate real general"<<std::endl;
00812       // Go to the beginning of the file
00813       istr.clear() ;
00814       istr.seekg(0, std::ios::beg);
00815     }
00816     skipComments(istr);
00817   
00818     std::size_t rows, cols, entries;
00819   
00820     if(lineFeed(istr))
00821       throw MatrixMarketFormatError();
00822   
00823     istr >> rows;
00824   
00825     if(lineFeed(istr))
00826       throw MatrixMarketFormatError();
00827     istr >> cols;
00828   
00829     if(lineFeed(istr))
00830       throw MatrixMarketFormatError();
00831   
00832     istr >>entries;
00833   
00834     std::size_t nnz, blockrows, blockcols;
00835 
00836     Dune::tie(blockrows, blockcols, nnz) = calculateNNZ<brows, bcols>(rows, cols, entries, header);
00837   
00838     istr.ignore(std::numeric_limits<std::streamsize>::max(),'\n');
00839 
00840   
00841     matrix.setSize(blockrows, blockcols);
00842     matrix.setBuildMode(Dune::BCRSMatrix<Dune::FieldMatrix<T,brows,bcols>,A>::row_wise);
00843   
00844     if(header.type==array_type)
00845       DUNE_THROW(Dune::NotImplemented, "Array format currently not supported for matrices!");
00846     
00847     readSparseEntries(matrix, istr, entries, header, NumericWrapper<double>());
00848   }     
00849   
00850   template<typename M>
00851   struct mm_multipliers
00852   {};
00853   
00854   template<typename B, int i, int j, typename A>
00855   struct mm_multipliers<BCRSMatrix<FieldMatrix<B,i,j>,A> >
00856   {
00857     enum{
00858       rows = i,
00859       cols = j
00860     };
00861   };
00862     
00863   template<typename B, int i, int j>
00864   void mm_print_entry(const FieldMatrix<B,i,j>& entry, 
00865                       typename FieldMatrix<B,i,j>::size_type rowidx,
00866                       typename FieldMatrix<B,i,j>::size_type colidx,
00867                       std::ostream& ostr)
00868   {
00869     typedef typename FieldMatrix<B,i,j>::const_iterator riterator;
00870     typedef typename FieldMatrix<B,i,j>::ConstColIterator citerator;
00871     for(riterator row=entry.begin(); row != entry.end(); ++row, ++rowidx){
00872       int coli=colidx;
00873       for(citerator col = row->begin(); col != row->end(); ++col, ++coli)
00874         ostr<< rowidx<<" "<<coli<<" "<<*col<<std::endl;
00875     }
00876   }
00877   
00878   // Write a vector entry
00879   template<typename V>
00880   void mm_print_vector_entry(const V& entry, std::ostream& ostr, 
00881                              const integral_constant<int,1>&)
00882   {
00883     ostr<<entry<<std::endl;
00884   }
00885   
00886   // Write a vector
00887   template<typename V>
00888   void mm_print_vector_entry(const V& vector, std::ostream& ostr, 
00889                              const integral_constant<int,0>&)
00890   {
00891     // Is the entry a supported numeric type?
00892     const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
00893     typedef typename V::const_iterator VIter;
00894     
00895     for(VIter i=vector.begin(); i != vector.end(); ++i)
00896         
00897       mm_print_vector_entry(*i, ostr, 
00898                             integral_constant<int,isnumeric>());
00899   }
00900 
00901   template<typename T, typename A, int i>
00902   std::size_t countEntries(const BlockVector<FieldVector<T,i>,A>& vector)
00903   {
00904     return vector.size()*i;
00905   }
00906 
00907   // Version for writing vectors.
00908   template<typename V>
00909   void writeMatrixMarket(const V& vector, std::ostream& ostr,
00910                          const integral_constant<int,0>&)
00911   {
00912     ostr<<countEntries(vector)<<" "<<1<<std::endl;
00913     const int isnumeric = mm_numeric_type<typename V::block_type>::is_numeric;
00914     mm_print_vector_entry(vector,ostr, integral_constant<int,isnumeric>());
00915   }
00916   
00917   // Versions for writing matrices
00918   template<typename M>
00919   void writeMatrixMarket(const M& matrix,
00920                          std::ostream& ostr,
00921                          const integral_constant<int,1>&)
00922   {
00923     ostr<<matrix.M()*mm_multipliers<M>::rows<<" "
00924         <<matrix.N()*mm_multipliers<M>::cols<<" "
00925         <<countNonZeros(matrix)<<std::endl;
00926     
00927     typedef typename M::const_iterator riterator;
00928     typedef typename M::ConstColIterator citerator;
00929     for(riterator row=matrix.begin(); row != matrix.end(); ++row)
00930       for(citerator col = row->begin(); col != row->end(); ++col)
00931         // Matrix Market indexing start with 1!
00932         mm_print_entry(*col, row.index()*mm_multipliers<M>::rows+1, 
00933                        col.index()*mm_multipliers<M>::cols+1, ostr);
00934   }
00935 
00936 
00940   template<typename M>
00941   void writeMatrixMarket(const M& matrix,
00942                         std::ostream& ostr)
00943   {
00944     // Write header information
00945     mm_header_printer<M>::print(ostr);
00946     mm_block_structure_header<M>::print(ostr,matrix);
00947     // Choose the correct function for matrix and vector
00948     writeMatrixMarket(matrix,ostr,integral_constant<int,IsMatrix<M>::value>());
00949   }
00950 
00951   
00962   template<typename M>
00963   void storeMatrixMarket(const M& matrix,
00964                          std::string filename)
00965   {
00966     std::ofstream file(filename.c_str());
00967     file.setf(std::ios::scientific,std::ios::floatfield);
00968     writeMatrixMarket(matrix, file);
00969     file.close();
00970   }
00971   
00972 #if HAVE_MPI
00973 
00986   template<typename M, typename G, typename L>
00987   void storeMatrixMarket(const M& matrix,
00988                          std::string filename,
00989                          const OwnerOverlapCopyCommunication<G,L>& comm,
00990                          bool storeIndices=true)
00991   {
00992     // Get our rank
00993     int rank = comm.communicator().rank();
00994     // Write the local matrix
00995     std::ostringstream rfilename;
00996     rfilename<<filename <<"_"<<rank<<".mm";
00997     std::cout<< rfilename.str()<<std::endl;
00998     std::ofstream file(rfilename.str().c_str());
00999     file.setf(std::ios::scientific,std::ios::floatfield);
01000     writeMatrixMarket(matrix, file);
01001     file.close();
01002     
01003     if(!storeIndices)
01004       return;
01005     
01006     // Write the global to local index mapping
01007     rfilename.str("");
01008     rfilename<<filename<<"_"<<rank<<".idx";
01009     file.open(rfilename.str().c_str());
01010     file.setf(std::ios::scientific,std::ios::floatfield);
01011     typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
01012     typedef typename IndexSet::const_iterator Iterator;
01013     for(Iterator iter = comm.indexSet().begin(); 
01014         iter != comm.indexSet().end(); ++iter){
01015       file << iter->global()<<" "<<(std::size_t)iter->local()<<" "
01016            <<(int)iter->local().attribute()<<" "<<(int)iter->local().isPublic()<<std::endl;
01017     }
01018     // Store neighbour information for efficient remote indices setup.
01019     file<<"neighbours:";
01020     const std::set<int>& neighbours=comm.remoteIndices().getNeighbours();
01021     typedef std::set<int>::const_iterator SIter;
01022     for(SIter neighbour=neighbours.begin(); neighbour != neighbours.end(); ++neighbour){
01023       file<<" "<< *neighbour;
01024     }
01025     file.close();    
01026   }
01027 
01042   template<typename M, typename G, typename L>
01043   void loadMatrixMarket(M& matrix,
01044                         const std::string& filename,
01045                         OwnerOverlapCopyCommunication<G,L>& comm,
01046                         bool readIndices=true)
01047   {
01048     typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet::LocalIndex LocalIndex;
01049     typedef typename LocalIndex::Attribute Attribute;
01050     // Get our rank
01051     int rank = comm.communicator().rank();
01052     // load local matrix
01053     std::ostringstream rfilename;
01054     rfilename<<filename <<"_"<<rank<<".mm";
01055     std::ifstream file;
01056     file.open(rfilename.str().c_str(), std::ios::in);
01057     if(!file)
01058       DUNE_THROW(IOError, "Could not open file" << rfilename.str().c_str());
01059     //if(!file.is_open())
01060     //
01061     readMatrixMarket(matrix,file);
01062     file.close();
01063     
01064     if(!readIndices)
01065       return;
01066     
01067     // read indices
01068     typedef typename OwnerOverlapCopyCommunication<G,L>::ParallelIndexSet IndexSet;
01069     IndexSet& pis=comm.pis;
01070     rfilename.str("");
01071     rfilename<<filename<<"_"<<rank<<".idx";
01072     file.open(rfilename.str().c_str());
01073     if(pis.size()!=0)
01074       DUNE_THROW(InvalidIndexSetState, "Index set is not empty!");
01075     
01076     pis.beginResize();
01077     while(!file.eof() && file.peek()!='n'){
01078       G g;
01079       file >>g;
01080       std::size_t l;
01081       file >>l;
01082       int c;
01083       file >>c;
01084       bool b;
01085       file >> b;
01086       pis.add(g,LocalIndex(l,Attribute(c),b));
01087       lineFeed(file);
01088     }
01089     pis.endResize();
01090     if(!file.eof()){
01091       // read neighbours
01092       std::string s;
01093       file>>s;
01094       if(s!="neighbours:")
01095         DUNE_THROW(MatrixMarketFormatError, "was exspecting the string: \"neighbours:\"");
01096       std::set<int> nb;
01097       while(!file.eof()){
01098         int i;
01099         file >> i;
01100         nb.insert(i);
01101       }
01102       file.close();
01103       comm.ri.setNeighbours(nb);
01104     }
01105     comm.ri.template rebuild<false>();
01106   }
01107 
01108   #endif
01109   
01120   template<typename M>
01121   void loadMatrixMarket(M& matrix,
01122                         const std::string& filename)
01123   {
01124     std::ifstream file;
01125     file.open(filename.c_str(), std::ios::in);
01126     if(!file)
01127       DUNE_THROW(IOError, "Could not open file" << filename);
01128     readMatrixMarket(matrix,file);
01129     file.close();
01130   }
01131   
01133 }
01134 #endif