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