Ifpack Package Browser (Single Doxygen Collection) Development
read_matrix.cpp
Go to the documentation of this file.
00001 #include <fstream>
00002 
00003 #include "Teuchos_CommHelpers.hpp"
00004 #include "Epetra_Comm.h"
00005 #include "Epetra_Map.h"
00006 #include "Epetra_CrsMatrix.h"
00007 #include "Epetra_Vector.h"
00008 #include "Trilinos_Util.h"
00009 
00010 Epetra_CrsMatrix*
00011 read_matrix_mm(const std::string& mm_file,
00012                const Epetra_Comm& comm)
00013 {
00014   int my_proc = comm.MyPID();
00015 
00016   int num_global_rows = 0;
00017   int nnz_per_row = 0;
00018 
00019   std::ifstream* infile = NULL;
00020   infile = new std::ifstream(mm_file.c_str());
00021   if (infile == NULL || !*infile) {
00022     throw std::runtime_error("Failed to open file "+mm_file);
00023   }
00024 
00025   std::ifstream& in = *infile;
00026 
00027   //first skip over the file header, which has
00028   //lines beginning with '%'.
00029   std::string line;
00030   do {
00031     getline(in, line);
00032   } while(line[0] == '%');
00033 
00034   //now get the matrix dimensions.
00035 
00036   int numrows, numcols, nnz;
00037   std::istringstream isstr(line);
00038   isstr >> numrows >> numcols >> nnz;
00039 
00040   //make sure we successfully read the three ints from that line.
00041   if (isstr.fail()) {
00042     throw std::runtime_error("Failed to parse matrix-market header.");
00043   }
00044 
00045   if (my_proc == 0) {
00046     num_global_rows = numrows;
00047     nnz_per_row = nnz/numrows;
00048   }
00049 
00050   comm.Broadcast(&num_global_rows, 1, 0);
00051   comm.Broadcast(&nnz_per_row, 1, 0);
00052 
00053   const int indexBase = 0;
00054   Epetra_Map rowmap(num_global_rows, indexBase, comm);
00055 
00056   Epetra_CrsMatrix* A = new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);
00057 
00058   Teuchos::Array<int> col;
00059   Teuchos::Array<double> coef;
00060 
00061   int irow=0, icol=0;
00062   int g_row=-1, last_row=-1;
00063   double val=0;
00064 
00065   while(!in.eof()) {
00066     getline(in, line);
00067     std::istringstream isstr(line);
00068     isstr >> irow >> icol >> val;
00069   
00070     if (isstr.fail()) continue;
00071     if (!rowmap.MyGID(irow-1)) continue;
00072 
00073     g_row = irow-1;
00074     if (g_row != last_row) {
00075       if (col.size() > 0) {
00076         A->InsertGlobalValues(last_row, col.size(), &coef[0], &col[0] );
00077         col.clear();
00078         coef.clear();
00079       }
00080       last_row = g_row;
00081     }
00082     col.push_back(icol-1);
00083     coef.push_back(val);
00084   }
00085 
00086   if (col.size() > 0) {
00087     A->InsertGlobalValues(g_row, col.size(), &coef[0], &col[0]);
00088   }
00089 
00090   A->FillComplete();
00091 
00092   return A;
00093 }
00094 
00095 Epetra_Vector*
00096 read_vector_mm(const std::string& mm_file,
00097                const Epetra_Comm& comm)
00098 {
00099   int my_proc = comm.MyPID();
00100 
00101   int num_global_rows = 0;
00102 
00103   std::ifstream* infile = NULL;
00104   if (my_proc == 0) {
00105     infile = new std::ifstream(mm_file.c_str());
00106     if (infile == NULL || !*infile) {
00107       throw std::runtime_error("Failed to open file "+mm_file);
00108     }
00109 
00110     std::ifstream& in = *infile;
00111 
00112     //first skip over the file header, which has
00113     //lines beginning with '%'.
00114     std::string line;
00115     do {
00116       getline(in, line);
00117     } while(line[0] == '%');
00118 
00119     //now get the matrix dimensions.
00120 
00121     int numrows, numcols;
00122     std::istringstream isstr(line);
00123     isstr >> numrows >> numcols;
00124 
00125     //make sure we successfully read the ints from that line.
00126     if (isstr.fail()) {
00127       throw std::runtime_error("Failed to parse matrix-market header.");
00128     }
00129 
00130     num_global_rows = numrows;
00131   }
00132 
00133   comm.Broadcast(&num_global_rows, 1, 0);
00134 
00135   const int indexBase = 0;
00136   Epetra_Map rowmap(num_global_rows, indexBase, comm);
00137 
00138   Epetra_Vector* b = new Epetra_Vector(rowmap, 1);
00139 
00140   if (my_proc == 0) {
00141     int irow=0, icol=0;
00142     double val=0;
00143 
00144     std::string line;
00145     std::ifstream& in = *infile;
00146     while(!in.eof()) {
00147       getline(in, line);
00148       std::istringstream isstr(line);
00149       isstr >> val;
00150     
00151       if (isstr.fail()) continue;
00152 
00153       b->ReplaceGlobalValue(irow++, icol, val);
00154     }
00155   }
00156 
00157   return b;
00158 }
00159 
00160 void read_matrix_hb(const std::string& hb_file,
00161                     const Epetra_Comm& Comm,
00162                     Epetra_CrsMatrix*& A,
00163                     Epetra_Vector*& b)
00164 {
00165   Epetra_Map* Map = NULL;
00166   Epetra_Vector* x = NULL;
00167   Epetra_Vector* xexact = NULL;
00168   Trilinos_Util_ReadHb2Epetra(const_cast<char*>(hb_file.c_str()), Comm, Map,
00169                              A, x, b, xexact);
00170   delete x;
00171   delete xexact;
00172 }
00173 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines