Tpetra parallel linear algebra  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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 
57 bool 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 
68 bool 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 
102 void 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 
154 void 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 
173 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) {
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 
328 namespace 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