EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
EpetraExt_MultiVectorOut.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_BlockMap.h"
45#include "Epetra_Map.h"
46#include "Epetra_Vector.h"
47#include "Epetra_IntVector.h"
54#include "Epetra_Import.h"
55#include "Epetra_CrsMatrix.h"
56
57using namespace EpetraExt;
58namespace EpetraExt {
59
60int MultiVectorToMatlabFile( const char *filename, const Epetra_MultiVector & A) {
61
62 std::FILE * handle = 0;
63 if (A.Map().Comm().MyPID()==0) { // Only PE 0 does this section
64 handle = fopen(filename,"w");
65 if (!handle) return(-1);
66 }
67 if (MultiVectorToMatlabHandle(handle, A)) return(-1); // Everybody calls this routine
68
69 if (A.Map().Comm().MyPID()==0) // Only PE 0 opened a file
70 if (fclose(handle)) return(-1);
71 return(0);
72}
73
74int MultiVectorToMatrixMarketFile( const char *filename, const Epetra_MultiVector & A,
75 const char * matrixName,
76 const char *matrixDescription,
77 bool writeHeader) {
78 long long M = A.GlobalLength64();
79 int N = A.NumVectors();
80
81 FILE * handle = 0;
82
83 if (A.Map().Comm().MyPID()==0) { // Only PE 0 does this section
84
85 handle = fopen(filename,"w");
86 if (!handle) return(-1);
87 MM_typecode matcode;
88 mm_initialize_typecode(&matcode);
89 mm_set_matrix(&matcode);
90 mm_set_array(&matcode);
91 mm_set_real(&matcode);
92
93 if (writeHeader==true) { // Only write header if requested (true by default)
94
95 if (mm_write_banner(handle, matcode)) return(-1);
96
97 if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
98 if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
99
100 if (mm_write_mtx_array_size(handle, M, N)) return(-1);
101 }
102 }
103
104 if (MultiVectorToMatrixMarketHandle(handle, A)) return(-1); // Everybody calls this routine
105
106 if (A.Map().Comm().MyPID()==0) // Only PE 0 opened a file
107 if (fclose(handle)) return(-1);
108 return(0);
109}
110
111int MultiVectorToMatlabHandle(FILE * handle, const Epetra_MultiVector & A) {
112 return(MultiVectorToHandle(handle, A, false));
113}
115 return(MultiVectorToHandle(handle, A, true));
116}
117
118template<typename int_type>
119int TMultiVectorToHandle(FILE * handle, const Epetra_MultiVector & A, bool mmFormat) {
120
121 Epetra_BlockMap bmap = A.Map();
122 const Epetra_Comm & comm = bmap.Comm();
123 int numProc = comm.NumProc();
124
125 if (numProc==1)
126 writeMultiVector(handle, A, mmFormat);
127 else {
128
129 // In the case of more than one column in the multivector, and writing to MatrixMarket
130 // format, we call this function recursively, passing each vector of the multivector
131 // individually so that we can get all of it written to file before going on to the next
132 // multivector
133 if (A.NumVectors()>1 && mmFormat) {
134 for (int i=0; i<A.NumVectors(); i++)
135 if (MultiVectorToHandle(handle, *(A(i)), mmFormat)) return(-1);
136 return(0);
137 }
138
139 Epetra_Map map((int_type) -1, bmap.NumMyPoints(), (int_type) 0, comm);
140 // Create a veiw of this multivector using a map (instead of block map)
141 Epetra_MultiVector A1(View, map, A.Pointers(), A.NumVectors());
142 int numRows = map.NumMyElements();
143
144 Epetra_Map allGidsMap((int_type) -1, numRows, (int_type) 0,comm);
145
146 typename Epetra_GIDTypeVector<int_type>::impl allGids(allGidsMap);
147 for (int i=0; i<numRows; i++) allGids[i] = (int_type) map.GID64(i);
148
149 // Now construct a MultiVector on PE 0 by strip-mining the rows of the input matrix A.
150 int numChunks = numProc;
151 int stripSize = allGids.GlobalLength64()/numChunks;
152 int remainder = allGids.GlobalLength64()%numChunks;
153 int curStart = 0;
154 int curStripSize = 0;
156 if (comm.MyPID()==0)
157 importGidList.Size(stripSize+1); // Set size of vector to max needed
158 for (int i=0; i<numChunks; i++) {
159 if (comm.MyPID()==0) { // Only PE 0 does this part
160 curStripSize = stripSize;
161 if (i<remainder) curStripSize++; // handle leftovers
162 for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
163 curStart += curStripSize;
164 }
165 // The following import map will be non-trivial only on PE 0.
166 Epetra_Map importGidMap((int_type) -1, curStripSize, importGidList.Values(), (int_type) 0, comm);
167 Epetra_Import gidImporter(importGidMap, allGidsMap);
168 typename Epetra_GIDTypeVector<int_type>::impl importGids(importGidMap);
169 if (importGids.Import(allGids, gidImporter, Insert)) return(-1);
170
171 // importGids now has a list of GIDs for the current strip of matrix rows.
172 // Use these values to build another importer that will get rows of the matrix.
173
174 // The following import map will be non-trivial only on PE 0.
175 Epetra_Map importMap((int_type) -1, importGids.MyLength(), importGids.Values(), (int_type) 0, comm);
176 Epetra_Import importer(importMap, map);
177 Epetra_MultiVector importA(importMap, A1.NumVectors());
178 if (importA.Import(A1, importer, Insert)) return(-1);
179
180 // Finally we are ready to write this strip of the matrix to ostream
181 if (writeMultiVector(handle, importA, mmFormat)) return(-1);
182 }
183 }
184 return(0);
185}
186
187int MultiVectorToHandle(FILE * handle, const Epetra_MultiVector & A, bool mmFormat) {
188#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
189 if(A.Map().GlobalIndicesInt()) {
190 return TMultiVectorToHandle<int>(handle, A, mmFormat);
191 }
192 else
193#endif
194#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
195 if(A.Map().GlobalIndicesLongLong()) {
196 return TMultiVectorToHandle<long long>(handle, A, mmFormat);
197 }
198 else
199#endif
200 throw "EpetraExt::MultiVectorToHandle: GlobalIndices type unknown";
201}
202
203int writeMultiVector(FILE * handle, const Epetra_MultiVector & A, bool mmFormat) {
204
205 int ierr = 0;
206 long long length = A.GlobalLength64();
207 int numVectors = A.NumVectors();
208 const Epetra_Comm & comm = A.Map().Comm();
209 if (comm.MyPID()!=0) {
210 if (A.MyLength()!=0) ierr = -1;
211 }
212 else {
213 if (length!=A.MyLength()) ierr = -1;
214 for (int j=0; j<numVectors; j++) {
215 for (long long i=0; i<length; i++) {
216 double val = A[j][i];
217 if (mmFormat)
218 fprintf(handle, "%22.16e\n", val);
219 else
220 fprintf(handle, "%22.16e ", val);
221 }
222 if (!mmFormat) fprintf(handle, "%s", "\n");
223 }
224 }
225 int ierrGlobal;
226 comm.MinAll(&ierr, &ierrGlobal, 1); // If any processor has -1, all return -1
227 return(ierrGlobal);
228}
229} // namespace EpetraExt
230
#define mm_set_array(typecode)
char MM_typecode[4]
#define mm_initialize_typecode(typecode)
#define mm_set_matrix(typecode)
#define mm_set_real(typecode)
Insert
View
int NumMyPoints() const
long long GID64(int LID) const
bool GlobalIndicesInt() const
const Epetra_Comm & Comm() const
int NumMyElements() const
bool GlobalIndicesLongLong() const
virtual int NumProc() const=0
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const=0
virtual int MyPID() const=0
const Epetra_BlockMap & Map() const
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int NumVectors() const
int MyLength() const
double ** Pointers() const
long long GlobalLength64() const
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int MultiVectorToMatrixMarketHandle(FILE *handle, const Epetra_MultiVector &A)
int writeMultiVector(FILE *handle, const Epetra_MultiVector &A, bool mmFormat)
int mm_write_mtx_array_size(FILE *f, long long M, long long N)
int mm_write_banner(FILE *f, MM_typecode matcode)
int MultiVectorToMatlabFile(const char *filename, const Epetra_MultiVector &A)
Writes an Epetra_MultiVector object to a file that is compatible with Matlab.
int TMultiVectorToHandle(FILE *handle, const Epetra_MultiVector &A, bool mmFormat)
int MultiVectorToMatlabHandle(FILE *handle, const Epetra_MultiVector &A)
int MultiVectorToHandle(FILE *handle, const Epetra_MultiVector &A, bool mmFormat)
int MultiVectorToMatrixMarketFile(const char *filename, const Epetra_MultiVector &A, const char *matrixName, const char *matrixDescription, bool writeHeader)
Writes an Epetra_MultiVector object to a Matrix Market format file.
const int N