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