EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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
41 #include "EpetraExt_OperatorOut.h"
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"
48 #include "Epetra_SerialDenseVector.h"
49 #include "Epetra_IntSerialDenseVector.h"
50 #include "Epetra_Import.h"
51 #include "Epetra_CrsMatrix.h"
52 
53 using namespace EpetraExt;
54 namespace EpetraExt {
55 
56 int 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 
63 int 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 
111 int 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 }
174 int 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 }
198 int 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
#define mm_set_coordinate(typecode)
static Epetra_Map Create_Root_Map(const Epetra_Map &usermap, int root=0)
#define mm_set_real(typecode)
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...
#define mm_set_matrix(typecode)
bool UniqueGIDs() const
int mm_write_banner(FILE *f, MM_typecode matcode)
virtual const Epetra_Map & OperatorDomainMap() const =0
virtual int MyPID() const =0
virtual int SumAll(double *PartialSums, double *GlobalSums, int Count) const =0
virtual int Apply(const Epetra_MultiVector &X, Epetra_MultiVector &Y) const =0
virtual const Epetra_Map & OperatorRangeMap() const =0
int mm_write_mtx_crd_size(FILE *f, long long M, long long N, long long nz)
virtual const Epetra_Comm & Comm() const =0
int get_nz(const Epetra_Operator &A, long long &nz)
char MM_typecode[4]
int OperatorToHandle(FILE *handle, const Epetra_Operator &A)
int LID(int GID) const
bool MyGID(int GID_in) const
const Epetra_Comm & Comm() const
int OperatorToMatlabFile(const char *filename, const Epetra_Operator &A)
Writes an Epetra_Operator object to a file that is compatible with Matlab.
#define mm_initialize_typecode(typecode)
int writeOperatorStrip(FILE *handle, const Epetra_MultiVector &y, const Epetra_Map &rootDomainMap, const Epetra_Map &rootRangeMap, long long startColumn)