EpetraExt Package Browser (Single Doxygen Collection)  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_PutRowMatrix.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 
42 #include "EpetraExt_PutRowMatrix.h"
43 #include "Epetra_Comm.h"
44 #include "Epetra_Map.h"
45 #include "Epetra_Vector.h"
46 #include "Epetra_IntVector.h"
49 #include "Epetra_Import.h"
50 #include "Epetra_RowMatrix.h"
51 #include "Epetra_CrsMatrix.h"
52 
53 using namespace Matlab;
54 namespace Matlab {
55 
56 int CopyRowMatrix(mxArray* matlabA, const Epetra_RowMatrix& A) {
57  int valueCount = 0;
58  //int* valueCount = &temp;
59 
60  Epetra_Map map = A.RowMatrixRowMap();
61  const Epetra_Comm & comm = map.Comm();
62  int numProc = comm.NumProc();
63 
64  if (numProc==1)
65  DoCopyRowMatrix(matlabA, valueCount, A);
66  else {
67  int numRows = map.NumMyElements();
68 
69  //cout << "creating allGidsMap\n";
70  Epetra_Map allGidsMap(-1, numRows, 0,comm);
71  //cout << "done creating allGidsMap\n";
72 
73  Epetra_IntVector allGids(allGidsMap);
74  for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
75 
76  // Now construct a RowMatrix on PE 0 by strip-mining the rows of the input matrix A.
77  int numChunks = numProc;
78  int stripSize = allGids.GlobalLength()/numChunks;
79  int remainder = allGids.GlobalLength()%numChunks;
80  int curStart = 0;
81  int curStripSize = 0;
82  Epetra_IntSerialDenseVector importGidList;
83  int numImportGids = 0;
84  if (comm.MyPID()==0)
85  importGidList.Size(stripSize+1); // Set size of vector to max needed
86  for (int i=0; i<numChunks; i++) {
87  if (comm.MyPID()==0) { // Only PE 0 does this part
88  curStripSize = stripSize;
89  if (i<remainder) curStripSize++; // handle leftovers
90  for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
91  curStart += curStripSize;
92  }
93  // The following import map will be non-trivial only on PE 0.
94  //cout << "creating importGidMap\n";
95  Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
96  //cout << "done creating importGidMap\n";
97  Epetra_Import gidImporter(importGidMap, allGidsMap);
98  Epetra_IntVector importGids(importGidMap);
99  if (importGids.Import(allGids, gidImporter, Insert)) return(-1);
100 
101  // importGids now has a list of GIDs for the current strip of matrix rows.
102  // Use these values to build another importer that will get rows of the matrix.
103 
104  // The following import map will be non-trivial only on PE 0.
105  //cout << "creating importMap\n";
106  //cout << "A.RowMatrixRowMap().MinAllGID: " << A.RowMatrixRowMap().MinAllGID() << "\n";
107  Epetra_Map importMap(-1, importGids.MyLength(), importGids.Values(), A.RowMatrixRowMap().MinAllGID(), comm);
108  //cout << "done creating importMap\n";
109  Epetra_Import importer(importMap, map);
110  Epetra_CrsMatrix importA(Copy, importMap, 0);
111  if (importA.Import(A, importer, Insert)) return(-1);
112  if (importA.FillComplete()) return(-1);
113 
114  // Finally we are ready to write this strip of the matrix to ostream
115  if (DoCopyRowMatrix(matlabA, valueCount, importA)) return(-1);
116  }
117  }
118 
119  if (A.RowMatrixRowMap().Comm().MyPID() == 0) {
120  // set max cap
121  int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
122  matlabAcolumnIndicesPtr[A.NumGlobalRows()] = valueCount;
123  }
124 
125  return(0);
126 }
127 
128 int DoCopyRowMatrix(mxArray* matlabA, int& valueCount, const Epetra_RowMatrix& A) {
129  //cout << "doing DoCopyRowMatrix\n";
130  int ierr = 0;
131  int numRows = A.NumGlobalRows();
132  //cout << "numRows: " << numRows << "\n";
133  Epetra_Map rowMap = A.RowMatrixRowMap();
134  Epetra_Map colMap = A.RowMatrixColMap();
135  int minAllGID = rowMap.MinAllGID();
136 
137  const Epetra_Comm & comm = rowMap.Comm();
138  //cout << "did global setup\n";
139  if (comm.MyPID()!=0) {
140  if (A.NumMyRows()!=0) ierr = -1;
141  if (A.NumMyCols()!=0) ierr = -1;
142  }
143  else {
144  // declare and get initial values of all matlabA pointers
145  double* matlabAvaluesPtr = mxGetPr(matlabA);
146  int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
147  int* matlabArowIndicesPtr = mxGetIr(matlabA);
148 
149  // set all matlabA pointers to the proper offset
150  matlabAvaluesPtr += valueCount;
151  matlabArowIndicesPtr += valueCount;
152 
153  if (numRows!=A.NumMyRows()) ierr = -1;
156  //cout << "did proc0 setup\n";
157  for (int i=0; i<numRows; i++) {
158  //cout << "extracting a row\n";
159  int I = rowMap.GID(i);
160  int numEntries = 0;
161  if (A.ExtractMyRowCopy(i, values.Length(), numEntries,
162  values.Values(), indices.Values())) return(-1);
163  matlabAcolumnIndicesPtr[I - minAllGID] = valueCount; // set the starting index of column I
164  double* serialValuesPtr = values.Values();
165  for (int j=0; j<numEntries; j++) {
166  int J = colMap.GID(indices[j]);
167  *matlabAvaluesPtr = *serialValuesPtr++;
168  *matlabArowIndicesPtr = J;
169  // increment matlabA pointers
170  matlabAvaluesPtr++;
171  matlabArowIndicesPtr++;
172  valueCount++;
173  }
174  }
175  //cout << "proc0 row extraction for this chunck is done\n";
176  }
177 
178 /*
179  if (comm.MyPID() == 0) {
180  cout << "printing matlabA pointers\n";
181  double* matlabAvaluesPtr = mxGetPr(matlabA);
182  int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
183  int* matlabArowIndicesPtr = mxGetIr(matlabA);
184  for(int i=0; i < numRows; i++) {
185  for(int j=0; j < A.MaxNumEntries(); j++) {
186  cout << "*matlabAvaluesPtr: " << *matlabAvaluesPtr++ << " *matlabAcolumnIndicesPtr: " << *matlabAcolumnIndicesPtr++ << " *matlabArowIndicesPtr" << *matlabArowIndicesPtr++ << "\n";
187  }
188  }
189 
190  cout << "done printing matlabA pointers\n";
191  }
192  */
193 
194  int ierrGlobal;
195  comm.MinAll(&ierr, &ierrGlobal, 1); // If any processor has -1, all return -1
196  return(ierrGlobal);
197 }
198 
199 } // namespace Matlab
virtual const Epetra_Map & RowMatrixRowMap() const =0
int Size(int Length_in)
virtual int MinAll(double *PartialMins, double *GlobalMins, int Count) const =0
int * Values() const
int GlobalLength() 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 GID(int LID) const
int CopyRowMatrix(mxArray *matlabA, const Epetra_RowMatrix &A)
virtual int NumMyRows() const =0
int DoCopyRowMatrix(mxArray *matlabA, int &valueCount, const Epetra_RowMatrix &A)
int MinAllGID() const
const Epetra_Comm & Comm() const
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
virtual int NumGlobalRows() const =0
int MyLength() const