Ifpack Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
AztecOO/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 int 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 //first skip over the file header, which has
69 //lines beginning with '%'.
70 std::string line;
71 do {
72 getline(*infile, line);
73 } while(line[0] == '%');
74
75 //now get the matrix dimensions.
76
77 int numrows, numcols, nnz;
78 std::istringstream isstr(line);
79 isstr >> numrows >> numcols >> nnz;
80
81 //make sure we successfully read the three ints from that line.
82 if (isstr.fail()) {
83 throw std::runtime_error("Failed to parse matrix-market header.");
84 }
85
86 if (my_proc == 0) {
87 num_global_rows = numrows;
88 nnz_per_row = nnz/numrows;
89 }
90
91 comm.Broadcast(&num_global_rows, 1, 0);
92 comm.Broadcast(&nnz_per_row, 1, 0);
93
94 const int indexBase = 0;
95 Epetra_Map rowmap(num_global_rows, indexBase, comm);
96
97 Epetra_CrsMatrix* A = new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);
98
99 Teuchos::Array<int> col;
100 Teuchos::Array<double> coef;
101
102 int irow=0, icol=0;
103 int g_row=-1, last_row=-1;
104 double val=0;
105
106 while(!infile->eof()) {
107 getline(*infile, line);
108 std::istringstream isstr(line);
109 isstr >> irow >> icol >> val;
110
111 if (isstr.fail()) continue;
112 if (!rowmap.MyGID(irow-1)) continue;
113
114 g_row = irow-1;
115 if (g_row != last_row) {
116 if (col.size() > 0) {
117 A->InsertGlobalValues(last_row, col.size(), &coef[0], &col[0] );
118 col.clear();
119 coef.clear();
120 }
121 last_row = g_row;
122 }
123 col.push_back(icol-1);
124 coef.push_back(val);
125 }
126
127 if (col.size() > 0) {
128 A->InsertGlobalValues(g_row, col.size(), &coef[0], &col[0]);
129 }
130
131 A->FillComplete();
132 delete infile;
133 return A;
134}
135
137read_vector_mm(const std::string& mm_file,
138 const Epetra_Comm& comm)
139{
140 int my_proc = comm.MyPID();
141
142 int num_global_rows = 0;
143
144 std::ifstream* infile = NULL;
145 if (my_proc == 0) {
146 infile = new std::ifstream(mm_file.c_str());
147 if (infile == NULL || !*infile) {
148 throw std::runtime_error("Failed to open file "+mm_file);
149 }
150
151 //first skip over the file header, which has
152 //lines beginning with '%'.
153 std::string line;
154 do {
155 getline(*infile, line);
156 } while(line[0] == '%');
157
158 //now get the matrix dimensions.
159
160 int numrows, numcols;
161 std::istringstream isstr(line);
162 isstr >> numrows >> numcols;
163
164 //make sure we successfully read the ints from that line.
165 if (isstr.fail()) {
166 throw std::runtime_error("Failed to parse matrix-market header.");
167 }
168
169 num_global_rows = numrows;
170 }
171
172 comm.Broadcast(&num_global_rows, 1, 0);
173
174 const int indexBase = 0;
175 Epetra_Map rowmap(num_global_rows, indexBase, comm);
176
177 Epetra_Vector* b = new Epetra_Vector(rowmap, 1);
178
179 if (my_proc == 0) {
180 int irow=0, icol=0;
181 double val=0;
182
183 std::string line;
184 std::ifstream& in = *infile;
185 while(!infile->eof()) {
186 getline(in, line);
187 std::istringstream isstr(line);
188 isstr >> val;
189
190 if (isstr.fail()) continue;
191
192 b->ReplaceGlobalValue(irow++, icol, val);
193 }
194 }
195 delete infile;
196 return b;
197}
198
199void read_matrix_hb(const std::string& hb_file,
200 const Epetra_Comm& Comm,
202 Epetra_Vector*& b)
203{
204 Epetra_Map* Map = NULL;
205 Epetra_Vector* x = NULL;
206 Epetra_Vector* xexact = NULL;
207 Trilinos_Util_ReadHb2Epetra(const_cast<char*>(hb_file.c_str()), Comm, Map,
208 A, x, b, xexact);
209 delete Map;
210 delete x;
211 delete xexact;
212}
213
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
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)