EpetraExt Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
EpetraExt_OperatorOut.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_MultiVector.h"
47#include "Epetra_IntVector.h"
50#include "Epetra_Import.h"
51#include "Epetra_CrsMatrix.h"
52
53using namespace EpetraExt;
54namespace EpetraExt {
55
56int OperatorToMatlabFile( const char *filename, const Epetra_Operator & A) {
57
58 // Simple wrapper to make it clear what can be used to write to Matlab format
59 EPETRA_CHK_ERR(OperatorToMatrixMarketFile(filename, A, 0, 0, false));
60 return(0);
61}
62
63int OperatorToMatrixMarketFile( const char *filename, const Epetra_Operator & A,
64 const char * matrixName,
65 const char *matrixDescription,
66 bool writeHeader) {
67
68 const Epetra_Map & domainMap = A.OperatorDomainMap();
69 const Epetra_Map & rangeMap = A.OperatorRangeMap();
70
71 if (!domainMap.UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
72 if (!rangeMap.UniqueGIDs()) {EPETRA_CHK_ERR(-2);}
73
74 long long M = rangeMap.NumGlobalElements64();
75 long long N = domainMap.NumGlobalElements64();
76 long long nz = 0;
77
78 FILE * handle = 0;
79
80 // To get count of nonzero terms we do multiplies ...
81 if (get_nz(A, nz)) {EPETRA_CHK_ERR(-1);}
82
83 if (domainMap.Comm().MyPID()==0) { // Only PE 0 does this section
84
85 handle = fopen(filename,"w");
86 if (!handle) {EPETRA_CHK_ERR(-1);}
87 MM_typecode matcode;
88 mm_initialize_typecode(&matcode);
89 mm_set_matrix(&matcode);
90 mm_set_coordinate(&matcode);
91 mm_set_real(&matcode);
92
93
94 if (writeHeader==true) { // Only write header if requested (true by default)
95
96 if (mm_write_banner(handle, matcode)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
97
98 if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
99 if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
100
101 if (mm_write_mtx_crd_size(handle, M, N, nz)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}
102 }
103 }
104
105 if (OperatorToHandle(handle, A)!=0) {if (handle!=0) fclose(handle); EPETRA_CHK_ERR(-1);}// Everybody calls this routine
106
107 if (handle!=0) fclose(handle);
108 return(0);
109}
110
111int OperatorToHandle(FILE * handle, const Epetra_Operator & A) {
112
113 const Epetra_Map & domainMap = A.OperatorDomainMap();
114 const Epetra_Map & rangeMap = A.OperatorRangeMap();
115 long long N = domainMap.NumGlobalElements64();
116
117 //cout << "rangeMap = " << rangeMap << endl;
118 Epetra_Map rootRangeMap = Epetra_Util::Create_Root_Map(rangeMap);
119 //cout << "rootRangeMap = " << rootRangeMap << endl;
120 Epetra_Map rootDomainMap = Epetra_Util::Create_Root_Map(domainMap, -1); // Replicate on all processors
121 Epetra_Import importer(rootRangeMap, rangeMap);
122
123 int chunksize = 5; // Let's do multiple RHS at a time
124 long long numchunks = N/chunksize;
125 int rem = N%chunksize;
126
127 if (rem>0) {
128 Epetra_MultiVector xrem(domainMap, rem);
129 Epetra_MultiVector yrem(rangeMap, rem);
130 Epetra_MultiVector yrem1(rootRangeMap, rem);
131 // Put 1's in slots
132 for (int j=0; j<rem; j++) {
133 long long curGlobalCol = rootDomainMap.GID64(j); // Should return same value on all processors
134 if (domainMap.MyGID(curGlobalCol)) {
135 int curCol = domainMap.LID(curGlobalCol);
136 xrem[j][curCol] = 1.0;
137 }
138 }
139 EPETRA_CHK_ERR(A.Apply(xrem, yrem));
140 EPETRA_CHK_ERR(yrem1.Import(yrem, importer, Insert));
141 EPETRA_CHK_ERR(writeOperatorStrip(handle, yrem1, rootDomainMap, rootRangeMap, 0));
142 }
143
144 if (numchunks>0) {
145 Epetra_MultiVector x(domainMap, chunksize);
146 Epetra_MultiVector y(rangeMap, chunksize);
147 Epetra_MultiVector y1(rootRangeMap, chunksize);
148 for (long long ichunk = 0; ichunk<numchunks; ichunk++) {
149 long long startCol = ichunk*chunksize+rem;
150 // Put 1's in slots
151 for (int j=0; j<chunksize; j++) {
152 long long curGlobalCol = rootDomainMap.GID64(startCol+j); // Should return same value on all processors
153 if (domainMap.MyGID(curGlobalCol)){
154 int curCol = domainMap.LID(curGlobalCol);
155 x[j][curCol] = 1.0;
156 }
157 }
158 EPETRA_CHK_ERR(A.Apply(x, y));
159 EPETRA_CHK_ERR(y1.Import(y, importer, Insert));
160 EPETRA_CHK_ERR(writeOperatorStrip(handle, y1, rootDomainMap, rootRangeMap, startCol));
161 // Put 0's in slots
162 for (int j=0; j<chunksize; j++) {
163 long long curGlobalCol = rootDomainMap.GID64(startCol+j); // Should return same value on all processors
164 if (domainMap.MyGID(curGlobalCol)){
165 int curCol = domainMap.LID(curGlobalCol);
166 x[j][curCol] = 0.0;
167 }
168 }
169 }
170 }
171
172 return(0);
173}
174int writeOperatorStrip(FILE * handle, const Epetra_MultiVector & y, const Epetra_Map & rootDomainMap, const Epetra_Map & rootRangeMap, long long startColumn) {
175
176 long long numRows = y.GlobalLength64();
177 int numCols = y.NumVectors();
178 long long ioffset = 1 - rootRangeMap.IndexBase64(); // Matlab indices start at 1
179 long long joffset = 1 - rootDomainMap.IndexBase64(); // Matlab indices start at 1
180 if (y.Comm().MyPID()!=0) {
181 if (y.MyLength()!=0) {EPETRA_CHK_ERR(-1);}
182 }
183 else {
184 if (numRows!=y.MyLength()) {EPETRA_CHK_ERR(-1);}
185 for (int j=0; j<numCols; j++) {
186 long long J = rootDomainMap.GID64(j + startColumn) + joffset;
187 for (long long i=0; i<numRows; i++) {
188 double val = y[j][i];
189 if (val!=0.0) {
190 long long I = rootRangeMap.GID64(i) + ioffset;
191 fprintf(handle, "%lld %lld %22.16e\n", I, J, val);
192 }
193 }
194 }
195 }
196 return(0);
197}
198int get_nz(const Epetra_Operator & A, long long & nz) {
199
200 const Epetra_Map & domainMap = A.OperatorDomainMap();
201 const Epetra_Map & rangeMap = A.OperatorRangeMap();
202
203 long long N = domainMap.NumGlobalElements64();
204
205 Epetra_Map rootDomainMap = Epetra_Util::Create_Root_Map(domainMap, -1); // Replicate on all processors
206
207
208 int chunksize = 5; // Let's do multiple RHS at a time
209 long long numchunks = N/chunksize;
210 int rem = N%chunksize;
211
212 long long lnz = 0;
213 if (rem>0) {
214 Epetra_MultiVector xrem(domainMap, rem);
215 Epetra_MultiVector yrem(rangeMap, rem);
216 // Put 1's in slots
217 for (int j=0; j<rem; j++) {
218 long long curGlobalCol = rootDomainMap.GID64(j);
219 if (domainMap.MyGID(curGlobalCol)) xrem[j][domainMap.LID(curGlobalCol)] = 1.0;
220 }
221 EPETRA_CHK_ERR(A.Apply(xrem, yrem));
222 for (int j=0; j<rem; j++) {
223 int mylength = yrem.MyLength();
224 for (int i=0; i<mylength; i++)
225 if (yrem[j][i]!=0.0) lnz++;
226 }
227 }
228
229 if (numchunks>0) {
230 Epetra_MultiVector x(domainMap, chunksize);
231 Epetra_MultiVector y(rangeMap, chunksize);
232 for (long long ichunk = 0; ichunk<numchunks; ichunk++) {
233 long long startCol = ichunk*chunksize+rem;
234 // Put 1's in slots
235 for (int j=0; j<chunksize; j++) {
236 long long curGlobalCol = rootDomainMap.GID64(startCol+j);
237 if (domainMap.MyGID(curGlobalCol)) x[j][domainMap.LID(curGlobalCol)] = 1.0;
238 }
239 EPETRA_CHK_ERR(A.Apply(x, y));
240 for (int j=0; j<chunksize; j++) {
241 int mylength = y.MyLength();
242 for (int i=0; i<mylength; i++)
243 if (y[j][i]!=0.0) lnz++;
244 }
245 // Put 0's in slots
246 for (int j=0; j<chunksize; j++) {
247 long long curGlobalCol = rootDomainMap.GID64(startCol+j);
248 if (domainMap.MyGID(curGlobalCol)) x[j][domainMap.LID(curGlobalCol)] = 0.0;
249 }
250 }
251 }
252
253 // Sum up nonzero counts
254 EPETRA_CHK_ERR(A.Comm().SumAll(&lnz, &nz, 1));
255 return(0);
256}
257} // 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
#define EPETRA_CHK_ERR(a)
long long IndexBase64() const
int LID(int GID) const
long long GID64(int LID) const
bool UniqueGIDs() const
long long NumGlobalElements64() const
const Epetra_Comm & Comm() const
bool MyGID(int GID_in) const
virtual int MyPID() const=0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const=0
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
const Epetra_Comm & Comm() const
int NumVectors() const
int MyLength() const
long long GlobalLength64() const
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const=0
virtual const Epetra_Comm & Comm() const=0
virtual const Epetra_Map & OperatorDomainMap() const=0
virtual const Epetra_Map & OperatorRangeMap() const=0
static Epetra_Map Create_Root_Map(const Epetra_Map &usermap, int root=0)
EpetraExt::BlockCrsMatrix: A class for constructing a distributed block matrix.
int writeOperatorStrip(FILE *handle, const Epetra_MultiVector &y, const Epetra_Map &rootDomainMap, const Epetra_Map &rootRangeMap, long long startColumn)
int OperatorToMatlabFile(const char *filename, const Epetra_Operator &A)
Writes an Epetra_Operator object to a file that is compatible with Matlab.
int OperatorToHandle(FILE *handle, const Epetra_Operator &A)
int mm_write_banner(FILE *f, MM_typecode matcode)
int get_nz(const Epetra_Operator &A, long long &nz)
int mm_write_mtx_crd_size(FILE *f, long long M, long long N, long long nz)
const int N
int OperatorToMatrixMarketFile(const char *filename, const Epetra_Operator &A, const char *matrixName, const char *matrixDescription, bool writeHeader)
Writes an Epetra_Operator object to a Matrix Market format file, forming the coefficients by applying...