EpetraExt  Development
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EpetraExt_PutBlockMap.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_PutBlockMap.h"
43 #include "Epetra_Comm.h"
44 #include "Epetra_BlockMap.h"
45 #include "Epetra_Map.h"
46 #include "Epetra_IntVector.h"
47 #include "Epetra_IntSerialDenseVector.h"
48 #include "Epetra_Import.h"
49 
50 using namespace Matlab;
51 namespace Matlab {
52 
53 int CopyBlockMap(mxArray* matlabA, const Epetra_BlockMap& map) {
54 
55  int valueCount = 0;
56  const Epetra_Comm & comm = map.Comm();
57  int numProc = comm.NumProc();
58  bool doSizes = !map.ConstantElementSize();
59 
60  if (numProc==1) {
61  int * myElements = map.MyGlobalElements();
62  int * elementSizeList = 0;
63  if (doSizes) elementSizeList = map.ElementSizeList();
64  return(DoCopyBlockMap(matlabA, valueCount, map.NumGlobalElements(), myElements, elementSizeList, doSizes));
65  }
66 
67  int numRows = map.NumMyElements();
68 
69  Epetra_Map allGidsMap(-1, numRows, 0,comm);
70 
71  Epetra_IntVector allGids(allGidsMap);
72  for (int i=0; i<numRows; i++) allGids[i] = map.GID(i);
73 
74  Epetra_IntVector allSizes(allGidsMap);
75  for (int i=0; i<numRows; i++) allSizes[i] = map.ElementSize(i);
76 
77  // Now construct a Map on PE 0 by strip-mining the rows of the input matrix map.
78  int numChunks = numProc;
79  int stripSize = allGids.GlobalLength()/numChunks;
80  int remainder = allGids.GlobalLength()%numChunks;
81  int curStart = 0;
82  int curStripSize = 0;
83  Epetra_IntSerialDenseVector importGidList;
84  Epetra_IntSerialDenseVector importSizeList;
85  int numImportGids = 0;
86  if (comm.MyPID()==0) {
87  importGidList.Size(stripSize+1); // Set size of vector to max needed
88  if (doSizes) importSizeList.Size(stripSize+1); // Set size of vector to max needed
89  }
90  for (int i=0; i<numChunks; i++) {
91  if (comm.MyPID()==0) { // Only PE 0 does this part
92  curStripSize = stripSize;
93  if (i<remainder) curStripSize++; // handle leftovers
94  for (int j=0; j<curStripSize; j++) importGidList[j] = j + curStart;
95  curStart += curStripSize;
96  }
97  // The following import map will be non-trivial only on PE 0.
98  Epetra_Map importGidMap(-1, curStripSize, importGidList.Values(), 0, comm);
99  Epetra_Import gidImporter(importGidMap, allGidsMap);
100 
101  Epetra_IntVector importGids(importGidMap);
102  if (importGids.Import(allGids, gidImporter, Insert)) return(-1);
103  Epetra_IntVector importSizes(importGidMap);
104  if (doSizes) if (importSizes.Import(allSizes, gidImporter, Insert)) return(-1);
105 
106  // importGids (and importSizes, if non-trivial block map)
107  // now have a list of GIDs (and sizes, respectively) for the current strip of map.
108 
109  int * myElements = importGids.Values();
110  int * elementSizeList = 0;
111  if (doSizes) elementSizeList = importSizes.Values();
112  // Finally we are ready to write this strip of the map to file
113  if (comm.MyPID()==0) {
114  DoCopyBlockMap(matlabA, valueCount, importGids.MyLength(), myElements, elementSizeList, doSizes);
115  }
116  }
117  return(0);
118 }
119 
120 int DoCopyBlockMap(mxArray* matlabA, int& valueCount, int length, const int * v1, const int * v2, bool doSizes) {
121  // declare and get initial values of all matlabA pointers
122  double* matlabAvaluesPtr = mxGetPr(matlabA);
123  int* matlabAcolumnIndicesPtr = mxGetJc(matlabA);
124  int* matlabArowIndicesPtr = mxGetIr(matlabA);
125 
126  // set all matlabA pointers to the proper offset
127  matlabAvaluesPtr += valueCount;
128  matlabArowIndicesPtr += valueCount;
129  int numGidsDone = valueCount;
130  if (doSizes) {
131  numGidsDone /= 2;
132  }
133 
134  matlabAcolumnIndicesPtr += numGidsDone;
135 
136  for (int i=0; i<length; i++) {
137  *matlabAcolumnIndicesPtr++ = valueCount;
138  *matlabArowIndicesPtr++ = 0;
139  *matlabAvaluesPtr++ = v1[i]; // row GID of block map
140 
141  if (doSizes) {
142  *matlabAvaluesPtr++ = v2[i]; // size of block
143  valueCount += 2;
144  *matlabArowIndicesPtr++ = 1;
145  }
146  else {
147  valueCount++;
148  }
149  }
150 
151  return(0);
152 }
153 } // namespace Matlab
int NumGlobalElements() const
int ElementSize() const
int CopyBlockMap(mxArray *matlabA, const Epetra_BlockMap &map)
int MyGlobalElements(int *MyGlobalElementList) const
int Size(int Length_in)
bool ConstantElementSize() const
int * Values() const
int GlobalLength() const
int DoCopyBlockMap(mxArray *matlabA, int &valueCount, int length, const int *v1, const int *v2, bool doSizes)
virtual int MyPID() const =0
int NumMyElements() const
int GID(int LID) const
int * ElementSizeList() const
const Epetra_Comm & Comm() const
virtual int NumProc() const =0
int MyLength() const