Tpetra parallel linear algebra Version of the Day
Loading...
Searching...
No Matches
Tpetra_MatrixIO.cpp
1/*
2// @HEADER
3// ***********************************************************************
4//
5// Tpetra: Templated Linear Algebra Services Package
6// Copyright (2008) Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
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 "Tpetra_MatrixIO.hpp"
45
46#include <cstdio>
47#include <cstdlib>
48#include <cstring>
49#include <functional>
50#include <algorithm>
51#include <iterator>
52#include <exception>
53#include <string>
54#include <cctype>
55#include <fstream>
56
57bool Tpetra::Utils::parseIfmt(Teuchos::ArrayRCP<char> fmt, int &perline, int &width) {
58 TEUCHOS_TEST_FOR_EXCEPT(fmt.size() != 0 && fmt[fmt.size()-1] != '\0');
59 // parses integers n and d out of (nId)
60 bool error = true;
61 std::transform(fmt.begin(), fmt.end(), fmt, static_cast < int(*)(int) > (std::toupper));
62 if (std::sscanf(fmt.getRawPtr(),"(%dI%d)",&perline,&width) == 2) {
63 error = false;
64 }
65 return error;
66}
67
68bool Tpetra::Utils::parseRfmt(Teuchos::ArrayRCP<char> fmt, int &perline, int &width, int &prec, char &valformat) {
69 TEUCHOS_TEST_FOR_EXCEPT(fmt.size() != 0 && fmt[fmt.size()-1] != '\0');
70 std::transform(fmt.begin(), fmt.end(), fmt, static_cast < int(*)(int) > (std::toupper));
71 // find the first left paren '(' and the last right paren ')'
72 Teuchos::ArrayRCP<char>::iterator firstLeftParen = std::find( fmt.begin(), fmt.end(), '(');
73 Teuchos::ArrayRCP<char>::iterator lastRightParen = std::find(std::reverse_iterator<Teuchos::ArrayRCP<char>::iterator>(fmt.end()),
74 std::reverse_iterator<Teuchos::ArrayRCP<char>::iterator>(fmt.begin()),
75 ')').base()-1;
76 // select the substring between the parens, including them
77 // if neither was found, set the string to empty
78 if (firstLeftParen == fmt.end() || lastRightParen == fmt.begin()) {
79 fmt.resize(0 + 1);
80 fmt[0] = '\0';
81 }
82 else {
83 fmt += (firstLeftParen - fmt.begin());
84 size_t newLen = lastRightParen - firstLeftParen + 1;
85 fmt.resize(newLen + 1);
86 fmt[newLen] = '\0';
87 }
88 if (std::find(fmt.begin(),fmt.end(),'P') != fmt.end()) {
89 // not supported
90 return true;
91 }
92 bool error = true;
93 if (std::sscanf(fmt.getRawPtr(),"(%d%c%d.%d)",&perline,&valformat,&width,&prec) == 4) {
94 if (valformat == 'E' || valformat == 'D' || valformat == 'F') {
95 error = false;
96 }
97 }
98 return error;
99}
100
101
102void Tpetra::Utils::readHBHeader(std::ifstream &fin, Teuchos::ArrayRCP<char> &Title, Teuchos::ArrayRCP<char> &Key, Teuchos::ArrayRCP<char> &Type,
103 int &Nrow, int &Ncol, int &Nnzero, int &Nrhs,
104 Teuchos::ArrayRCP<char> &Ptrfmt, Teuchos::ArrayRCP<char> &Indfmt, Teuchos::ArrayRCP<char> &Valfmt, Teuchos::ArrayRCP<char> &Rhsfmt,
105 int &Ptrcrd, int &Indcrd, int &Valcrd, int &Rhscrd, Teuchos::ArrayRCP<char> &Rhstype) {
106 int Totcrd, Neltvl, Nrhsix;
107 const int MAXLINE = 81;
108 char line[MAXLINE];
109 //
110 Title.resize(72 + 1); std::fill(Title.begin(), Title.end(), '\0');
111 Key.resize(8 + 1); std::fill(Key.begin(), Key.end(), '\0');
112 Type.resize(3 + 1); std::fill(Type.begin(), Type.end(), '\0');
113 Ptrfmt.resize(16 + 1); std::fill(Ptrfmt.begin(), Ptrfmt.end(), '\0');
114 Indfmt.resize(16 + 1); std::fill(Indfmt.begin(), Indfmt.end(), '\0');
115 Valfmt.resize(20 + 1); std::fill(Valfmt.begin(), Valfmt.end(), '\0');
116 Rhsfmt.resize(20 + 1); std::fill(Rhsfmt.begin(), Rhsfmt.end(), '\0');
117 //
118 const std::string errStr("Tpetra::Utils::readHBHeader(): Improperly formatted H/B file: ");
119 /* First line: (A72,A8) */
120 fin.getline(line,MAXLINE);
121 TEUCHOS_TEST_FOR_EXCEPTION( std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
122 (void)std::sscanf(line, "%72c%8[^\n]", Title.getRawPtr(), Key.getRawPtr());
123 /* Second line: (5I14) or (4I14) */
124 fin.getline(line,MAXLINE);
125 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
126 if ( std::sscanf(line,"%14d%14d%14d%14d%14d",&Totcrd,&Ptrcrd,&Indcrd,&Valcrd,&Rhscrd) != 5 ) {
127 Rhscrd = 0;
128 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%14d%14d%14d%14d",&Totcrd,&Ptrcrd,&Indcrd,&Valcrd) != 4, std::runtime_error, errStr << "error reading pointers (line 2)");
129 }
130 /* Third line: (A3, 11X, 4I14) */
131 fin.getline(line,MAXLINE);
132 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
133 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line, "%3c%14i%14i%14i%14i", Type.getRawPtr(),&Nrow,&Ncol,&Nnzero,&Neltvl) != 5 , std::runtime_error, errStr << "error reading matrix meta-data (line 3)");
134 std::transform(Type.begin(), Type.end(), Type.begin(), static_cast < int(*)(int) > (std::toupper));
135 /* Fourth line: */
136 fin.getline(line,MAXLINE);
137 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
138 if (Rhscrd != 0) {
139 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%16c%16c%20c%20c",Ptrfmt.getRawPtr(),Indfmt.getRawPtr(),Valfmt.getRawPtr(),Rhsfmt.getRawPtr()) != 4, std::runtime_error, errStr << "error reading formats (line 4)");
140 }
141 else {
142 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%16c%16c%20c",Ptrfmt.getRawPtr(),Indfmt.getRawPtr(),Valfmt.getRawPtr()) != 3, std::runtime_error, errStr << "error reading formats (line 4)");
143 }
144 /* (Optional) Fifth line: */
145 if (Rhscrd != 0 ) {
146 Rhstype.resize(3 + 1,'\0');
147 fin.getline(line,MAXLINE);
148 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%*s") < 0, std::runtime_error, errStr << "error buffering line.");
149 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(line,"%3c%14d%14d", Rhstype.getRawPtr(), &Nrhs, &Nrhsix) != 3, std::runtime_error, errStr << "error reading right-hand-side meta-data (line 5)");
150 }
151}
152
153
154void Tpetra::Utils::readHBInfo(const std::string &filename, int &M, int &N, int &nz, Teuchos::ArrayRCP<char> &Type, int &Nrhs) {
155 std::ifstream fin;
156 int Ptrcrd, Indcrd, Valcrd, Rhscrd;
157 Teuchos::ArrayRCP<char> Title, Key, Rhstype, Ptrfmt, Indfmt, Valfmt, Rhsfmt;
158 try {
159 fin.open(filename.c_str(),std::ifstream::in);
160 Tpetra::Utils::readHBHeader(fin, Title, Key, Type, M, N, nz, Nrhs,
161 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
162 Ptrcrd, Indcrd, Valcrd, Rhscrd, Rhstype);
163 fin.close();
164 }
165 catch (std::exception &e) {
166 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
167 "Tpetra::Utils::readHBInfo() of filename \"" << filename << "\" caught exception: " << std::endl
168 << e.what() << std::endl);
169 }
170}
171
172
173void Tpetra::Utils::readHBMatDouble(const std::string &filename, int &numRows, int &numCols, int &numNZ, std::string &type, Teuchos::ArrayRCP<int> &colPtrs, Teuchos::ArrayRCP<int> &rowInds, Teuchos::ArrayRCP<double> &vals) {
174 // NOTE: if complex, allocate enough space for 2*NNZ and interleave real and imaginary parts (real,imag)
175 // if pattern, don't touch parameter vals; do not allocate space space for values
176 try {
177 std::ifstream fin;
178 int ptrCrd, indCrd, valCrd;
179 Teuchos::ArrayRCP<char> Title, Key, Ptrfmt, Indfmt, Valfmt;
180 const int MAXSIZE = 81;
181 char lineBuf[MAXSIZE];
182 // nitty gritty
183 int ptrsPerLine, ptrWidth, indsPerLine, indWidth, valsPerLine, valWidth, valPrec;
184 char valFlag;
185 //
186 fin.open(filename.c_str(),std::ifstream::in);
187 {
188 // we don't care about RHS-related stuff, so declare those vars in an expiring scope
189 int Nrhs, rhsCrd;
190 Teuchos::ArrayRCP<char> Rhstype, Rhsfmt;
191 Teuchos::ArrayRCP<char> TypeArray;
192 Tpetra::Utils::readHBHeader(fin, Title, Key, TypeArray, numRows, numCols, numNZ, Nrhs,
193 Ptrfmt, Indfmt, Valfmt, Rhsfmt,
194 ptrCrd, indCrd, valCrd, rhsCrd, Rhstype);
195 if (TypeArray.size() > 0) {
196 type.resize(TypeArray.size()-1);
197 std::copy(TypeArray.begin(), TypeArray.end(), type.begin());
198 }
199 }
200 const std::string errStr("Tpetra::Utils::readHBHeader(): Improperly formatted H/B file.");
201 const bool readPatternOnly = (type[0] == 'P' || type[0] == 'p');
202 const bool readComplex = (type[0] == 'C' || type[0] == 'c');
203 /* Parse the array input formats from Line 3 of HB file */
204 TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseIfmt(Ptrfmt,ptrsPerLine,ptrWidth) == true, std::runtime_error,
205 "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
206 TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseIfmt(Indfmt,indsPerLine,indWidth) == true, std::runtime_error,
207 "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
208 if (readPatternOnly == false) {
209 TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseRfmt(Valfmt,valsPerLine,valWidth,valPrec,valFlag) == true, std::runtime_error,
210 "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
211 }
212 // special case this: the reason is that the number of colPtrs read is numCols+1, which is non-zero even if numCols == 0
213 // furthermore, if numCols == 0, there is nothing of interest to read
214 if (numCols == 0) return;
215 // allocate space for column pointers, row indices, and values
216 // if the file is empty, do not touch these ARCPs
217 colPtrs = Teuchos::arcp<int>(numCols+1);
218 if (numNZ > 0) {
219 rowInds = Teuchos::arcp<int>(numNZ);
220 if (readPatternOnly == false) {
221 if (readComplex) {
222 vals = Teuchos::arcp<double>(2*numNZ);
223 }
224 else {
225 vals = Teuchos::arcp<double>(numNZ);
226 }
227 }
228 }
229 /* Read column pointer array:
230 Specifically, read ptrCrd number of lines, and on each line, read ptrsPerLine number of integers, each of width ptrWidth
231 Store these in colPtrs */
232 {
233 int colPtrsRead = 0;
234 char NullSub = '\0';
235 for (int lno=0; lno < ptrCrd; ++lno) {
236 fin.getline(lineBuf, MAXSIZE);
237 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
238 char *linePtr = lineBuf;
239 for (int ptr=0; ptr < ptrsPerLine; ++ptr) {
240 if (colPtrsRead == numCols + 1) break;
241 int cptr;
242 // terminate the string at the end of the current ptr block, saving the character in that location
243 std::swap(NullSub,linePtr[ptrWidth]);
244 // read the ptr
245 std::sscanf(linePtr, "%d", &cptr);
246 // put the saved character back, and put the '\0' back into NullSub for use again
247 std::swap(NullSub,linePtr[ptrWidth]);
248 linePtr += ptrWidth;
249 colPtrs[colPtrsRead++] = cptr;
250 }
251 }
252 TEUCHOS_TEST_FOR_EXCEPT(colPtrsRead != numCols + 1);
253 }
254 /* Read row index array:
255 Specifically, read indCrd number of lines, and on each line, read indsPerLine number of integers, each of width indWidth
256 Store these in rowInds */
257 {
258 char NullSub = '\0';
259 int indicesRead = 0;
260 for (int lno=0; lno < indCrd; ++lno) {
261 fin.getline(lineBuf, MAXSIZE);
262 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
263 char *linePtr = lineBuf;
264 for (int indcntr=0; indcntr < indsPerLine; ++indcntr) {
265 if (indicesRead == numNZ) break;
266 int ind;
267 // terminate the string at the end of the current ind block, saving the character in that location
268 std::swap(NullSub,linePtr[indWidth]);
269 // read the ind
270 std::sscanf(linePtr, "%d", &ind);
271 // put the saved character back, and put the '\0' back into NullSub for use again
272 std::swap(NullSub,linePtr[indWidth]);
273 linePtr += indWidth;
274 rowInds[indicesRead++] = ind;
275 }
276 }
277 TEUCHOS_TEST_FOR_EXCEPT(indicesRead != numNZ);
278 }
279 /* Read array of values:
280 Specifically, read valCrd number of lines, and on each line, read valsPerLine number of real values, each of width/precision valWidth/valPrec
281 Store these in vals
282 If readComplex, then read twice as many non-zeros, and interleave real,imag into vals */
283 if (readPatternOnly == false) {
284 int totalNumVals;
285 if (readComplex) {
286 totalNumVals = 2*numNZ;
287 }
288 else {
289 totalNumVals = numNZ;
290 }
291 char NullSub = '\0';
292 int valsRead = 0;
293 for (int lno=0; lno < valCrd; ++lno) {
294 fin.getline(lineBuf, MAXSIZE);
295 TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
296 // if valFlag == 'D', then we need to convert [dD] in fp vals into [eE] that scanf can parse
297 if (valFlag == 'D') std::replace_if(lineBuf, lineBuf+MAXSIZE, [] (const char c) { return c == 'D'; }, 'E');
298 char *linePtr = lineBuf;
299 for (int valcntr=0; valcntr < valsPerLine; ++valcntr) {
300 if (valsRead == totalNumVals) break;
301 double val;
302 // terminate the string at the end of the current val block, saving the character in that location
303 std::swap(NullSub,linePtr[valWidth]);
304 // read the val
305 std::sscanf(linePtr, "%le", &val);
306 // put the saved character back, and put the '\0' back into NullSub for use again
307 std::swap(NullSub,linePtr[valWidth]);
308 linePtr += valWidth;
309 vals[valsRead++] = val;
310 }
311 }
312 TEUCHOS_TEST_FOR_EXCEPT(valsRead != totalNumVals);
313 }
314 fin.close();
315 }
316 catch (std::exception &e) {
317 TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
318 "Tpetra::Utils::readHBInfo() of filename \"" << filename << "\" caught exception: " << std::endl
319 << e.what() << std::endl);
320 }
321}
322
323#ifdef HAVE_TPETRA_EXPLICIT_INSTANTIATION
324
325#include "TpetraCore_ETIHelperMacros.h"
326#include "Tpetra_MatrixIO_def.hpp"
327
328namespace Tpetra {
329 namespace Utils {
330
331 TPETRA_ETI_MANGLING_TYPEDEFS()
332
333 TPETRA_INSTANTIATE_SLGN(TPETRA_MATRIXIO_INSTANT)
334
335 } // namespace Tpetra::Utils
336} // namespace Tpetra
337
338#endif // HAVE_TPETRA_EXPLICIT_INSTANTIATION
Namespace Tpetra contains the class and methods constituting the Tpetra library.