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 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 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 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 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 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  Tpetra::Utils::readHBHeader(fin, Title, Key, Type, M, N, nz, Nrhs,
127  Ptrfmt, Indfmt, Valfmt, Rhsfmt,
128  Ptrcrd, Indcrd, Valcrd, Rhscrd, Rhstype);
129  fin.close();
130  }
131  catch (std::exception &e) {
132  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
133  "Tpetra::Utils::readHBInfo() of filename \"" << filename << "\" caught exception: " << std::endl
134  << e.what() << std::endl);
135  }
136 }
137 
138 
139 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) {
140  // NOTE: if complex, allocate enough space for 2*NNZ and interleave real and imaginary parts (real,imag)
141  // if pattern, don't touch parameter vals; do not allocate space space for values
142  try {
143  std::ifstream fin;
144  int ptrCrd, indCrd, valCrd;
145  Teuchos::ArrayRCP<char> Title, Key, Ptrfmt, Indfmt, Valfmt;
146  const int MAXSIZE = 81;
147  char lineBuf[MAXSIZE];
148  // nitty gritty
149  int ptrsPerLine, ptrWidth, indsPerLine, indWidth, valsPerLine, valWidth, valPrec;
150  char valFlag;
151  //
152  fin.open(filename.c_str(),std::ifstream::in);
153  {
154  // we don't care about RHS-related stuff, so declare those vars in an expiring scope
155  int Nrhs, rhsCrd;
156  Teuchos::ArrayRCP<char> Rhstype, Rhsfmt;
157  Teuchos::ArrayRCP<char> TypeArray;
158  Tpetra::Utils::readHBHeader(fin, Title, Key, TypeArray, numRows, numCols, numNZ, Nrhs,
159  Ptrfmt, Indfmt, Valfmt, Rhsfmt,
160  ptrCrd, indCrd, valCrd, rhsCrd, Rhstype);
161  if (TypeArray.size() > 0) {
162  type.resize(TypeArray.size()-1);
163  std::copy(TypeArray.begin(), TypeArray.end(), type.begin());
164  }
165  }
166  const std::string errStr("Tpetra::Utils::readHBHeader(): Improperly formatted H/B file.");
167  const bool readPatternOnly = (type[0] == 'P' || type[0] == 'p');
168  const bool readComplex = (type[0] == 'C' || type[0] == 'c');
169  /* Parse the array input formats from Line 3 of HB file */
170  TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseIfmt(Ptrfmt,ptrsPerLine,ptrWidth) == true, std::runtime_error,
171  "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
172  TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseIfmt(Indfmt,indsPerLine,indWidth) == true, std::runtime_error,
173  "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
174  if (readPatternOnly == false) {
175  TEUCHOS_TEST_FOR_EXCEPTION( Tpetra::Utils::parseRfmt(Valfmt,valsPerLine,valWidth,valPrec,valFlag) == true, std::runtime_error,
176  "Tpetra::Utils::readHBMatDouble(): error parsing. Invalid/unsupported file format.");
177  }
178  // special case this: the reason is that the number of colPtrs read is numCols+1, which is non-zero even if numCols == 0
179  // furthermore, if numCols == 0, there is nothing of interest to read
180  if (numCols == 0) return;
181  // allocate space for column pointers, row indices, and values
182  // if the file is empty, do not touch these ARCPs
183  colPtrs = Teuchos::arcp<int>(numCols+1);
184  if (numNZ > 0) {
185  rowInds = Teuchos::arcp<int>(numNZ);
186  if (readPatternOnly == false) {
187  if (readComplex) {
188  vals = Teuchos::arcp<double>(2*numNZ);
189  }
190  else {
191  vals = Teuchos::arcp<double>(numNZ);
192  }
193  }
194  }
195  /* Read column pointer array:
196  Specifically, read ptrCrd number of lines, and on each line, read ptrsPerLine number of integers, each of width ptrWidth
197  Store these in colPtrs */
198  {
199  int colPtrsRead = 0;
200  char NullSub = '\0';
201  for (int lno=0; lno < ptrCrd; ++lno) {
202  fin.getline(lineBuf, MAXSIZE);
203  TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
204  char *linePtr = lineBuf;
205  for (int ptr=0; ptr < ptrsPerLine; ++ptr) {
206  if (colPtrsRead == numCols + 1) break;
207  int cptr;
208  // terminate the string at the end of the current ptr block, saving the character in that location
209  std::swap(NullSub,linePtr[ptrWidth]);
210  // read the ptr
211  std::sscanf(linePtr, "%d", &cptr);
212  // put the saved character back, and put the '\0' back into NullSub for use again
213  std::swap(NullSub,linePtr[ptrWidth]);
214  linePtr += ptrWidth;
215  colPtrs[colPtrsRead++] = cptr;
216  }
217  }
218  TEUCHOS_TEST_FOR_EXCEPT(colPtrsRead != numCols + 1);
219  }
220  /* Read row index array:
221  Specifically, read indCrd number of lines, and on each line, read indsPerLine number of integers, each of width indWidth
222  Store these in rowInds */
223  {
224  char NullSub = '\0';
225  int indicesRead = 0;
226  for (int lno=0; lno < indCrd; ++lno) {
227  fin.getline(lineBuf, MAXSIZE);
228  TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
229  char *linePtr = lineBuf;
230  for (int indcntr=0; indcntr < indsPerLine; ++indcntr) {
231  if (indicesRead == numNZ) break;
232  int ind;
233  // terminate the string at the end of the current ind block, saving the character in that location
234  std::swap(NullSub,linePtr[indWidth]);
235  // read the ind
236  std::sscanf(linePtr, "%d", &ind);
237  // put the saved character back, and put the '\0' back into NullSub for use again
238  std::swap(NullSub,linePtr[indWidth]);
239  linePtr += indWidth;
240  rowInds[indicesRead++] = ind;
241  }
242  }
243  TEUCHOS_TEST_FOR_EXCEPT(indicesRead != numNZ);
244  }
245  /* Read array of values:
246  Specifically, read valCrd number of lines, and on each line, read valsPerLine number of real values, each of width/precision valWidth/valPrec
247  Store these in vals
248  If readComplex, then read twice as many non-zeros, and interleave real,imag into vals */
249  if (readPatternOnly == false) {
250  int totalNumVals;
251  if (readComplex) {
252  totalNumVals = 2*numNZ;
253  }
254  else {
255  totalNumVals = numNZ;
256  }
257  char NullSub = '\0';
258  int valsRead = 0;
259  for (int lno=0; lno < valCrd; ++lno) {
260  fin.getline(lineBuf, MAXSIZE);
261  TEUCHOS_TEST_FOR_EXCEPTION(std::sscanf(lineBuf,"%*s") < 0, std::runtime_error, errStr);
262  // if valFlag == 'D', then we need to convert [dD] in fp vals into [eE] that scanf can parse
263  if (valFlag == 'D') std::replace_if(lineBuf, lineBuf+MAXSIZE, [] (const char c) { return c == 'D'; }, 'E');
264  char *linePtr = lineBuf;
265  for (int valcntr=0; valcntr < valsPerLine; ++valcntr) {
266  if (valsRead == totalNumVals) break;
267  double val;
268  // terminate the string at the end of the current val block, saving the character in that location
269  std::swap(NullSub,linePtr[valWidth]);
270  // read the val
271  std::sscanf(linePtr, "%le", &val);
272  // put the saved character back, and put the '\0' back into NullSub for use again
273  std::swap(NullSub,linePtr[valWidth]);
274  linePtr += valWidth;
275  vals[valsRead++] = val;
276  }
277  }
278  TEUCHOS_TEST_FOR_EXCEPT(valsRead != totalNumVals);
279  }
280  fin.close();
281  }
282  catch (std::exception &e) {
283  TEUCHOS_TEST_FOR_EXCEPTION(true, std::runtime_error,
284  "Tpetra::Utils::readHBInfo() of filename \"" << filename << "\" caught exception: " << std::endl
285  << e.what() << std::endl);
286  }
287 }
288 
289 #ifdef HAVE_TPETRA_EXPLICIT_INSTANTIATION
290 
291 #include "TpetraCore_ETIHelperMacros.h"
292 #include "Tpetra_MatrixIO_def.hpp"
293 
294 namespace Tpetra {
295  namespace Utils {
296 
297  TPETRA_ETI_MANGLING_TYPEDEFS()
298 
299  TPETRA_INSTANTIATE_SLGN(TPETRA_MATRIXIO_INSTANT)
300 
301  } // namespace Tpetra::Utils
302 } // namespace Tpetra
303 
304 #endif // HAVE_TPETRA_EXPLICIT_INSTANTIATION
Generated on Thu Nov 21 2024 09:19:26 for Tpetra parallel linear algebra by   doxygen 1.8.5