EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_RowMatrixOut.cpp
Go to the documentation of this file.
1 //@HEADER
2 // ***********************************************************************
3 //
4 // EpetraExt: Epetra Extended - Linear Algebra Services Package
5 // Copyright (2011) Sandia Corporation
6 //
7 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
8 // the U.S. Government retains certain rights in this software.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 #include "EpetraExt_RowMatrixOut.h"
42 #include "EpetraExt_mmio.h"
43 #include "Epetra_Comm.h"
44 #include "Epetra_Map.h"
45 #include "Epetra_Vector.h"
46 #include "Epetra_IntVector.h"
47 #include "Epetra_LongLongVector.h"
48 #include "Epetra_GIDTypeVector.h"
49 #include "Epetra_SerialDenseVector.h"
50 #include "Epetra_IntSerialDenseVector.h"
51 #include "Epetra_LongLongSerialDenseVector.h"
52 #include "Epetra_GIDTypeSerialDenseVector.h"
53 #include "Epetra_Import.h"
54 #include "Epetra_CrsMatrix.h"
55 #include <limits>
56 
57 using namespace EpetraExt;
58 namespace EpetraExt {
59 
60 int RowMatrixToMatlabFile( const char *filename, const Epetra_RowMatrix & A) {
61 
62  // Simple wrapper to make it clear what can be used to write to Matlab format
63  return(RowMatrixToMatrixMarketFile(filename, A, 0, 0, false));
64 }
65 
66 int RowMatrixToMatrixMarketFile( const char *filename, const Epetra_RowMatrix & A,
67  const char * matrixName,
68  const char *matrixDescription,
69  bool writeHeader) {
70  long long M = A.NumGlobalRows64();
71  long long N = A.NumGlobalCols64();
72  long long nz = A.NumGlobalNonzeros64();
73 
74  FILE * handle = 0;
75 
76  if (A.RowMatrixRowMap().Comm().MyPID()==0) { // Only PE 0 does this section
77 
78  handle = fopen(filename,"w");
79  if (!handle) {EPETRA_CHK_ERR(-1);}
80  MM_typecode matcode;
81  mm_initialize_typecode(&matcode);
82  mm_set_matrix(&matcode);
83  mm_set_coordinate(&matcode);
84  mm_set_real(&matcode);
85 
86  if (writeHeader==true) { // Only write header if requested (true by default)
87 
88  if (mm_write_banner(handle, matcode)!=0) {EPETRA_CHK_ERR(-1);}
89 
90  if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
91  if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
92 
93  if (mm_write_mtx_crd_size(handle, M, N, nz)!=0) {EPETRA_CHK_ERR(-1);}
94  }
95  }
96 
97  if (RowMatrixToHandle(handle, A)!=0) {EPETRA_CHK_ERR(-1);}// Everybody calls this routine
98 
99  if (A.RowMatrixRowMap().Comm().MyPID()==0) // Only PE 0 opened a file
100  if (fclose(handle)!=0) {EPETRA_CHK_ERR(-1);}
101  return(0);
102 }
103 
104 template<typename int_type>
105 int RowMatrixToHandle(FILE * handle, const Epetra_RowMatrix & A) {
106 
107  Epetra_Map map = A.RowMatrixRowMap();
108  const Epetra_Comm & comm = map.Comm();
109  int numProc = comm.NumProc();
110 
111  if (numProc==1 || !A.Map().DistributedGlobal())
112  writeRowMatrix(handle, A);
113  else {
114  int numRows = map.NumMyElements();
115 
116  Epetra_Map allGidsMap((int_type) -1, numRows, (int_type) 0,comm);
117 
118  typename Epetra_GIDTypeVector<int_type>::impl allGids(allGidsMap);
119  for (int i=0; i<numRows; i++) allGids[i] = (int_type) map.GID64(i);
120 
121  // Now construct a RowMatrix on PE 0 by strip-mining the rows of the input matrix A.
122  int numChunks = numProc;
123  int stripSize = allGids.GlobalLength64()/numChunks;
124  int remainder = allGids.GlobalLength64()%numChunks;
125  int curStart = 0;
126  int curStripSize = 0;
128  if (comm.MyPID()==0)
129  importGidList.Size(stripSize+1); // Set size of vector to max needed
130  for (int i=0; i<numChunks; i++) {
131  if (comm.MyPID()==0) { // Only PE 0 does this part
132  curStripSize = stripSize;
133  if (i<remainder) curStripSize++; // handle leftovers
134  for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
135  curStart += curStripSize;
136  }
137  // The following import map will be non-trivial only on PE 0.
138  if (comm.MyPID()>0) assert(curStripSize==0);
139  Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
140  Epetra_Import gidImporter(importGidMap, allGidsMap);
141  typename Epetra_GIDTypeVector<int_type>::impl importGids(importGidMap);
142  if (importGids.Import(allGids, gidImporter, Insert)!=0) {EPETRA_CHK_ERR(-1); }
143 
144  // importGids now has a list of GIDs for the current strip of matrix rows.
145  // Use these values to build another importer that will get rows of the matrix.
146 
147  // The following import map will be non-trivial only on PE 0.
148  Epetra_Map importMap(-1, importGids.MyLength(), importGids.Values(), map.IndexBase64(), comm);
149  Epetra_Import importer(importMap, map);
150  Epetra_CrsMatrix importA(Copy, importMap, 0);
151  if (importA.Import(A, importer, Insert)!=0) {EPETRA_CHK_ERR(-1); }
152  if (importA.FillComplete(A.OperatorDomainMap(), importMap)!=0) {EPETRA_CHK_ERR(-1);}
153 
154  // Finally we are ready to write this strip of the matrix to ostream
155  if (writeRowMatrix(handle, importA)!=0) {EPETRA_CHK_ERR(-1);}
156  }
157  }
158  return(0);
159 }
160 
161 int RowMatrixToHandle(FILE * handle, const Epetra_RowMatrix & A) {
162 #ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
164  return RowMatrixToHandle<int>(handle, A);
165  }
166  else
167 #endif
168 #ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
170  return RowMatrixToHandle<long long>(handle, A);
171  }
172  else
173 #endif
174  throw "EpetraExt::RowMatrixToHandle: GlobalIndices type unknown";
175 }
176 
177 int writeRowMatrix(FILE * handle, const Epetra_RowMatrix & A) {
178 
179  long long numRows_LL = A.NumGlobalRows64();
180  if(numRows_LL > std::numeric_limits<int>::max())
181  throw "EpetraExt::writeRowMatrix: numRows_LL > std::numeric_limits<int>::max()";
182 
183  int numRows = static_cast<int>(numRows_LL);
184  Epetra_Map rowMap = A.RowMatrixRowMap();
185  Epetra_Map colMap = A.RowMatrixColMap();
186  const Epetra_Comm & comm = rowMap.Comm();
187  long long ioffset = 1 - rowMap.IndexBase64(); // Matlab indices start at 1
188  long long joffset = 1 - colMap.IndexBase64(); // Matlab indices start at 1
189  if (comm.MyPID()!=0) {
190  if (A.NumMyRows()!=0) {EPETRA_CHK_ERR(-1);}
191  if (A.NumMyCols()!=0) {EPETRA_CHK_ERR(-1);}
192  }
193  else {
194  if (numRows!=A.NumMyRows()) {EPETRA_CHK_ERR(-1);}
197  for (int i=0; i<numRows; i++) {
198  long long I = rowMap.GID64(i) + ioffset;
199  int numEntries;
200  if (A.ExtractMyRowCopy(i, values.Length(), numEntries,
201  values.Values(), indices.Values())!=0) {EPETRA_CHK_ERR(-1);}
202  for (int j=0; j<numEntries; j++) {
203  long long J = colMap.GID64(indices[j]) + joffset;
204  double val = values[j];
205  fprintf(handle, "%lld %lld %22.16e\n", I, J, val);
206  }
207  }
208  }
209  return(0);
210 }
211 } // namespace EpetraExt
int RowMatrixToHandle(FILE *handle, const Epetra_RowMatrix &A)
bool DistributedGlobal() const
#define mm_set_coordinate(typecode)
virtual const Epetra_Map & RowMatrixRowMap() const =0
#define mm_set_real(typecode)
bool GlobalIndicesLongLong() const
#define mm_set_matrix(typecode)
int mm_write_banner(FILE *f, MM_typecode matcode)
virtual const Epetra_Map & OperatorDomainMap() const =0
bool GlobalIndicesInt() const
virtual int MyPID() const =0
int FillComplete(bool OptimizeDataStorage=true)
virtual int NumMyCols() const =0
int NumMyElements() const
virtual int MaxNumEntries() const =0
int mm_write_mtx_crd_size(FILE *f, long long M, long long N, long long nz)
int RowMatrixToMatlabFile(const char *filename, const Epetra_RowMatrix &A)
Writes an Epetra_RowMatrix object to a file that is compatible with Matlab.
virtual const Epetra_BlockMap & Map() const =0
virtual int NumMyRows() const =0
char MM_typecode[4]
const Epetra_Comm & Comm() const
int writeRowMatrix(FILE *handle, const Epetra_RowMatrix &A)
int RowMatrixToMatrixMarketFile(const char *filename, const Epetra_RowMatrix &A, const char *matrixName, const char *matrixDescription, bool writeHeader)
Writes an Epetra_RowMatrix object to a Matrix Market format file.
virtual int NumProc() const =0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const =0
virtual const Epetra_Map & RowMatrixColMap() const =0
#define mm_initialize_typecode(typecode)