EpetraExt Development
Loading...
Searching...
No Matches
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
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
50using namespace Matlab;
51namespace Matlab {
52
53int 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
120int 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
Insert
int MyGlobalElements(int *MyGlobalElementList) const
int * ElementSizeList() const
int GID(int LID) const
int ElementSize() const
int NumGlobalElements() const
const Epetra_Comm & Comm() const
bool ConstantElementSize() const
int NumMyElements() const
virtual int NumProc() const=0
virtual int MyPID() const=0
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
int Size(int Length_in)
int MyLength() const
int * Values() const
int GlobalLength() const
int CopyBlockMap(mxArray *matlabA, const Epetra_BlockMap &map)
int DoCopyBlockMap(mxArray *matlabA, int &valueCount, int length, const int *v1, const int *v2, bool doSizes)