Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
23 bool 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 
34 bool 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 
68 void 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 
125 void 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 
143 void 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 
297 namespace Tpetra {
298 namespace Utils {
299 
300 TPETRA_ETI_MANGLING_TYPEDEFS()
301 
302 TPETRA_INSTANTIATE_SLGN(TPETRA_MATRIXIO_INSTANT)
303 
304 } // namespace Utils
305 } // namespace Tpetra
306 
307 #endif // HAVE_TPETRA_EXPLICIT_INSTANTIATION
Generated on Thu Oct 2 2025 07:10:52 for Tpetra parallel linear algebra by   doxygen 1.8.5