Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
AztecOO_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
54read_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
101 Teuchos::Array<long long> col;
102 Teuchos::Array<double> coef;
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
139read_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
203void 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
Epetra_CrsMatrix * read_matrix_mm(const std::string &mm_file, const Epetra_Comm &comm)
void read_matrix_hb(const std::string &hb_file, const Epetra_Comm &Comm, Epetra_CrsMatrix *&A, Epetra_Vector *&b)
Epetra_Vector * read_vector_mm(const std::string &mm_file, const Epetra_Comm &comm)
Copy
bool MyGID(int GID_in) const
virtual int Broadcast(double *MyVals, int Count, int Root) const=0
virtual int MyPID() const=0
int FillComplete(bool OptimizeDataStorage=true)
virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double *Values, const int *Indices)
int ReplaceGlobalValue(int GlobalRow, int VectorIndex, double ScalarValue)