Ifpack Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LL/read_matrix.cpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // Ifpack: Object-Oriented Algebraic Preconditioner Package
6 // Copyright (2002) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 #include <fstream>
45 
46 #include "Teuchos_CommHelpers.hpp"
47 #include "Epetra_Comm.h"
48 #include "Epetra_Map.h"
49 #include "Epetra_CrsMatrix.h"
50 #include "Epetra_Vector.h"
51 #include "Trilinos_Util.h"
52 
54 read_matrix_mm(const std::string& mm_file,
55  const Epetra_Comm& comm)
56 {
57  int my_proc = comm.MyPID();
58 
59  long long num_global_rows = 0;
60  int nnz_per_row = 0;
61 
62  std::ifstream* infile = NULL;
63  infile = new std::ifstream(mm_file.c_str());
64  if (infile == NULL || !*infile) {
65  throw std::runtime_error("Failed to open file "+mm_file);
66  }
67 
68  std::ifstream& in = *infile;
69 
70  //first skip over the file header, which has
71  //lines beginning with '%'.
72  std::string line;
73  do {
74  getline(in, line);
75  } while(line[0] == '%');
76 
77  //now get the matrix dimensions.
78 
79  int numrows, numcols, nnz;
80  std::istringstream isstr(line);
81  isstr >> numrows >> numcols >> nnz;
82 
83  //make sure we successfully read the three ints from that line.
84  if (isstr.fail()) {
85  throw std::runtime_error("Failed to parse matrix-market header.");
86  }
87 
88  if (my_proc == 0) {
89  num_global_rows = numrows;
90  nnz_per_row = nnz/numrows;
91  }
92 
93  comm.Broadcast(&num_global_rows, 1, 0);
94  comm.Broadcast(&nnz_per_row, 1, 0);
95 
96  const int indexBase = 0;
97  Epetra_Map rowmap(num_global_rows, indexBase, comm);
98 
99  Epetra_CrsMatrix* A = new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);
100 
103 
104  int irow=0, icol=0;
105  int g_row=-1, last_row=-1;
106  double val=0;
107 
108  while(!in.eof()) {
109  getline(in, line);
110  std::istringstream isstr(line);
111  isstr >> irow >> icol >> val;
112 
113  if (isstr.fail()) continue;
114  if (!rowmap.MyGID(irow-1)) continue;
115 
116  g_row = irow-1;
117  if (g_row != last_row) {
118  if (col.size() > 0) {
119  A->InsertGlobalValues(last_row, col.size(), &coef[0], &col[0] );
120  col.clear();
121  coef.clear();
122  }
123  last_row = g_row;
124  }
125  col.push_back(icol-1);
126  coef.push_back(val);
127  }
128 
129  if (col.size() > 0) {
130  A->InsertGlobalValues(g_row, col.size(), &coef[0], &col[0]);
131  }
132 
133  A->FillComplete();
134 
135  return A;
136 }
137 
139 read_vector_mm(const std::string& mm_file,
140  const Epetra_Comm& comm)
141 {
142  int my_proc = comm.MyPID();
143 
144  long long num_global_rows = 0;
145 
146  std::ifstream* infile = NULL;
147  if (my_proc == 0) {
148  infile = new std::ifstream(mm_file.c_str());
149  if (infile == NULL || !*infile) {
150  throw std::runtime_error("Failed to open file "+mm_file);
151  }
152 
153  std::ifstream& in = *infile;
154 
155  //first skip over the file header, which has
156  //lines beginning with '%'.
157  std::string line;
158  do {
159  getline(in, line);
160  } while(line[0] == '%');
161 
162  //now get the matrix dimensions.
163 
164  int numrows, numcols;
165  std::istringstream isstr(line);
166  isstr >> numrows >> numcols;
167 
168  //make sure we successfully read the ints from that line.
169  if (isstr.fail()) {
170  throw std::runtime_error("Failed to parse matrix-market header.");
171  }
172 
173  num_global_rows = numrows;
174  }
175 
176  comm.Broadcast(&num_global_rows, 1, 0);
177 
178  const int indexBase = 0;
179  Epetra_Map rowmap(num_global_rows, indexBase, comm);
180 
181  Epetra_Vector* b = new Epetra_Vector(rowmap, 1);
182 
183  if (my_proc == 0) {
184  int irow=0, icol=0;
185  double val=0;
186 
187  std::string line;
188  std::ifstream& in = *infile;
189  while(!in.eof()) {
190  getline(in, line);
191  std::istringstream isstr(line);
192  isstr >> val;
193 
194  if (isstr.fail()) continue;
195 
196  b->ReplaceGlobalValue(irow++, icol, val);
197  }
198  }
199 
200  return b;
201 }
202 
203 void read_matrix_hb(const std::string& hb_file,
204  const Epetra_Comm& Comm,
206  Epetra_Vector*& b)
207 {
208  Epetra_Map* Map = NULL;
209  Epetra_Vector* x = NULL;
210  Epetra_Vector* xexact = NULL;
211  Trilinos_Util_ReadHb2Epetra64(const_cast<char*>(hb_file.c_str()), Comm, Map,
212  A, x, b, xexact);
213  delete x;
214  delete xexact;
215 }
216 
void read_matrix_hb(const std::string &hb_file, const Epetra_Comm &Comm, Epetra_CrsMatrix *&A, Epetra_Vector *&b)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
bool MyGID(int GID_in) const
virtual int Broadcast(double *MyVals, int Count, int Root) const =0
void push_back(const value_type &x)
size_type size() const
Epetra_Vector * read_vector_mm(const std::string &mm_file, const Epetra_Comm &comm)
Epetra_CrsMatrix * read_matrix_mm(const std::string &mm_file, const Epetra_Comm &comm)
Definition: read_matrix.cpp:54