EpetraExt Development
Loading...
Searching...
No Matches
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
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
57using namespace EpetraExt;
58namespace EpetraExt {
59
60int 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
66int 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
104template<typename int_type>
105int 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
161int 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
177int 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
char MM_typecode[4]
#define mm_set_coordinate(typecode)
#define mm_initialize_typecode(typecode)
#define mm_set_matrix(typecode)
#define mm_set_real(typecode)
Insert
Copy
bool DistributedGlobal() const
bool GlobalIndicesInt() const
const Epetra_Comm & Comm() const
int NumMyElements() const
bool GlobalIndicesLongLong() const
virtual int NumProc() const=0
virtual int MyPID() const=0
int FillComplete(bool OptimizeDataStorage=true)
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual int NumMyRows() const=0
virtual int NumMyCols() const=0
virtual const Epetra_Map & RowMatrixColMap() const=0
virtual const Epetra_Map & RowMatrixRowMap() const=0
virtual int MaxNumEntries() const=0
virtual int ExtractMyRowCopy(int MyRow, int Length, int &NumEntries, double *Values, int *Indices) const=0
double * Values() const
virtual const Epetra_BlockMap & Map() const=0
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int RowMatrixToMatlabFile(const char *filename, const Epetra_RowMatrix &A)
Writes an Epetra_RowMatrix object to a file that is compatible with Matlab.
int mm_write_banner(FILE *f, MM_typecode matcode)
int writeRowMatrix(FILE *handle, const Epetra_RowMatrix &A)
int RowMatrixToHandle(FILE *handle, const Epetra_RowMatrix &A)
int mm_write_mtx_crd_size(FILE *f, long long M, long long N, long long nz)
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.