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