Epetra Package Browser (Single Doxygen Collection) Development
Loading...
Searching...
No Matches
Epetra_LinearProblemRedistor.cpp
Go to the documentation of this file.
1
2//@HEADER
3// ************************************************************************
4//
5// Epetra: Linear Algebra Services Package
6// Copyright 2011 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39//
40// ************************************************************************
41//@HEADER
42
43#include "Epetra_Comm.h"
44#include "Epetra_Map.h"
46#include "Epetra_CrsMatrix.h"
49#include "Epetra_Util.h"
50#include "Epetra_MultiVector.h"
51#include "Epetra_IntVector.h"
52#include "Epetra_Export.h"
53#include "Epetra_Import.h"
54//=============================================================================
56 const Epetra_Map & RedistMap)
57 : OrigProblem_(OrigProblem),
58 NumProc_(0),
59 RedistProblem_(0),
60 RedistMap_((Epetra_Map *) &RedistMap),
61 Transposer_(0),
62 Replicate_(false),
63 ConstructTranspose_(false),
64 MakeDataContiguous_(false),
65 MapGenerated_(false),
66 RedistProblemCreated_(false),
67 ptr_(0)
68{
69}
70//=============================================================================
72 int NumProc,
73 bool Replicate)
74 : OrigProblem_(OrigProblem),
75 NumProc_(NumProc),
76 RedistProblem_(0),
77 RedistMap_(0),
78 Transposer_(0),
79 Replicate_(Replicate),
80 ConstructTranspose_(false),
81 MakeDataContiguous_(false),
82 MapGenerated_(false),
83 RedistProblemCreated_(false),
84 ptr_(0)
85{
86}
87//=============================================================================
89 : OrigProblem_(Source.OrigProblem_),
90 NumProc_(Source.NumProc_),
91 RedistProblem_(Source.RedistProblem_),
92 RedistMap_(Source.RedistMap_),
93 Transposer_(Source.Transposer_),
94 Replicate_(Source.Replicate_),
95 ConstructTranspose_(Source.ConstructTranspose_),
96 MakeDataContiguous_(Source.MakeDataContiguous_),
97 RedistProblemCreated_(Source.RedistProblemCreated_),
98 ptr_(0)
99{
100
103}
104//=========================================================================
106
107 if (ptr_!=0) {delete [] ptr_; ptr_=0;}
108 if (MapGenerated_ && RedistMap_!=0) {delete RedistMap_; RedistMap_=0;}
109 if (RedistProblem_!=0) {
110 // If no tranpose, then we must delete matrix (otherwise the transposer must).
112 delete RedistProblem_->GetMatrix();
113 if (RedistProblem_->GetLHS()!=0) delete RedistProblem_->GetLHS();
114 if (RedistProblem_->GetRHS()!=0) delete RedistProblem_->GetRHS();
116 }
118 if (Transposer_!=0) {delete Transposer_; Transposer_ = 0;}
119
120}
121
122//=========================================================================
124
125
126 if (MapGenerated_) return(0);
127
128 const Epetra_Map & SourceMap = OrigProblem_->GetMatrix()->RowMatrixRowMap();
129 const Epetra_Comm & Comm = SourceMap.Comm();
130 int IndexBase = SourceMap.IndexBase();
131
132 int NumProc = Comm.NumProc();
133
134 if (NumProc_!=NumProc) return(-1); // Right now we are only supporting redistribution to all processors.
135
136 // Build a list of contiguous GIDs that will be used in either case below.
137 int NumMyRedistElements = 0;
138 if ((Comm.MyPID()==0) || Replicate_) NumMyRedistElements = SourceMap.NumGlobalElements64();
139
140 // Now build the GID list to broadcast SourceMapGIDs to all processors that need it (or just to PE 0).
141 int * ContigIDs = 0;
142 if (NumMyRedistElements>0) ContigIDs = new int[NumMyRedistElements];
143 for (int i=0; i<NumMyRedistElements; i++) ContigIDs[i] = IndexBase + i;
144
145 // Case 1: If the map for the input matrix is not a linear contiguous map, then we have to collect
146 // the map indices in order to construct a compatible map containing all GIDs.
147 if (!SourceMap.LinearMap()) {
148
149 // First generate a linear map of the same distribution as RowMatrixRowMap
150 Epetra_Map SourceLinearMap(-1, SourceMap.NumMyElements(), IndexBase, Comm);
151 // Generate Int vector containing GIDs of SourceMap
152 Epetra_IntVector SourceMapGIDs(View, SourceLinearMap, SourceMap.MyGlobalElements());
153 // Now Build target map for SourceMapGIDs and Importer to get IDs
154 Epetra_Map GIDsTargetMap(-1, NumMyRedistElements, ContigIDs, IndexBase, Comm);
155 if (NumMyRedistElements>0) delete [] ContigIDs;
156
157 Epetra_Import GIDsImporter(GIDsTargetMap, SourceMap);
158 Epetra_IntVector TargetMapGIDs(GIDsTargetMap);
159
160 // Now Send SourceMap GIDs to PE 0, and all processors if Replicate is true
161 EPETRA_CHK_ERR(TargetMapGIDs.Import(SourceMapGIDs, GIDsImporter, Insert));
162
163 // Finally, create RedistMap containing all GIDs of SourceMap on PE 0, or all PEs of Replicate is true
164 RedistMap_ = new Epetra_Map(-1, NumMyRedistElements, TargetMapGIDs.Values(), IndexBase, Comm);// FIXME long long
165
166 }
167 // Case 2: If the map has contiguous IDs then we can simply build the map right away using the list
168 // of contiguous GIDs directly
169 else
170 RedistMap_ = new Epetra_Map(-1, NumMyRedistElements, ContigIDs, IndexBase, Comm);// FIXME long long
171
172 MapGenerated_ = true;
173
174 return(0);
175}
176
177//=========================================================================
179 const bool MakeDataContiguous,
180 Epetra_LinearProblem *& RedistProblem) {
181
182 if (RedistProblemCreated_) EPETRA_CHK_ERR(-1); // This method can only be called once
183
184 Epetra_RowMatrix * OrigMatrix = OrigProblem_->GetMatrix();
187
188 if (OrigMatrix==0) EPETRA_CHK_ERR(-2); // There is no matrix associated with this Problem
189
190
191 if (RedistMap_==0) {
193 }
194
196
198 Epetra_CrsMatrix * RedistMatrix;
199
200 // Check if the tranpose should be create or not
201 if (ConstructTranspose) {
202 Transposer_ = new Epetra_RowMatrixTransposer(OrigMatrix);
203 EPETRA_CHK_ERR(Transposer_->CreateTranspose(MakeDataContiguous, RedistMatrix, RedistMap_));
204 }
205 else {
206 // If not, then just do the redistribution based on the the RedistMap
207 RedistMatrix = new Epetra_CrsMatrix(Copy, *RedistMap_, 0);
208 // need to do this next step until we generalize the Import/Export ops for CrsMatrix
209 Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix);
210 EPETRA_CHK_ERR(RedistMatrix->Export(*OrigCrsMatrix, *RedistExporter_, Add));
211 EPETRA_CHK_ERR(RedistMatrix->FillComplete());
212 }
213
214 RedistProblem_->SetOperator(RedistMatrix);
215
216 // Now redistribute the RHS and LHS if non-zero
217
218 Epetra_MultiVector * RedistLHS = 0;
219 Epetra_MultiVector * RedistRHS = 0;
220
221 int ierr = 0;
222
223 if (OrigLHS!=0) {
224 RedistLHS = new Epetra_MultiVector(*RedistMap_, OrigLHS->NumVectors());
225 EPETRA_CHK_ERR(RedistLHS->Export(*OrigLHS, *RedistExporter_, Add));
226 }
227 else ierr = 1;
228
229 if (OrigRHS!=0) {
230 RedistRHS = new Epetra_MultiVector(*RedistMap_, OrigLHS->NumVectors());
231 EPETRA_CHK_ERR(RedistRHS->Export(*OrigRHS, *RedistExporter_, Add));
232 }
233 else ierr ++;
234
235 RedistProblem_->SetLHS(RedistLHS);
236 RedistProblem_->SetRHS(RedistRHS);
237
239
240 return(ierr);
241}
242
243//=========================================================================
245
246 if (!RedistProblemCreated_) EPETRA_CHK_ERR(-1); // This method can only be called after CreateRedistProblem()
247
248 Epetra_RowMatrix * OrigMatrix = ProblemWithNewValues->GetMatrix();
249 Epetra_MultiVector * OrigLHS = ProblemWithNewValues->GetLHS();
250 Epetra_MultiVector * OrigRHS = ProblemWithNewValues->GetRHS();
251
252 if (OrigMatrix==0) EPETRA_CHK_ERR(-2); // There is no matrix associated with this Problem
253
254
255 Epetra_CrsMatrix * RedistMatrix = dynamic_cast<Epetra_CrsMatrix *>(RedistProblem_->GetMatrix());
256
257 // Check if the tranpose should be create or not
260 }
261 else {
262 // If not, then just do the redistribution based on the the RedistMap
263
264 EPETRA_CHK_ERR(RedistMatrix->PutScalar(0.0));
265 // need to do this next step until we generalize the Import/Export ops for CrsMatrix
266 Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix);
267
268 if (OrigCrsMatrix==0) EPETRA_CHK_ERR(-3); // Broken for a RowMatrix at this point
269 EPETRA_CHK_ERR(RedistMatrix->Export(*OrigCrsMatrix, *RedistExporter_, Add));
270 }
271
272 // Now redistribute the RHS and LHS if non-zero
273
274
275 if (OrigLHS!=0) {
277 }
278
279 if (OrigRHS!=0) {
281 }
282
283 return(0);
284}
285
286//=========================================================================
287// NOTE: This method should be removed and replaced with calls to Epetra_Util_ExtractHbData()
288int Epetra_LinearProblemRedistor::ExtractHbData(int & M, int & N, int & nz, int * & ptr,
289 int * & ind, double * & val, int & Nrhs,
290 double * & rhs, int & ldrhs,
291 double * & lhs, int & ldlhs) const {
292
293 Epetra_CrsMatrix * RedistMatrix = dynamic_cast<Epetra_CrsMatrix *>(RedistProblem_->GetMatrix());
294
295 if (RedistMatrix==0) EPETRA_CHK_ERR(-1); // This matrix is zero or not an Epetra_CrsMatrix
296 if (!RedistMatrix->IndicesAreContiguous()) { // Data must be contiguous for this to work
297 EPETRA_CHK_ERR(-2);
298 }
299
300 M = RedistMatrix->NumMyRows();
301 N = RedistMatrix->NumMyCols();
302 nz = RedistMatrix->NumMyNonzeros();
303 val = (*RedistMatrix)[0]; // Dangerous, but cheap and effective way to access first element in
304
305 const Epetra_CrsGraph & RedistGraph = RedistMatrix->Graph();
306 ind = RedistGraph[0]; // list of values and indices
307
310 Nrhs = RHS->NumVectors();
311 if (Nrhs>1) {
312 if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-3)}; // Must have strided vectors
313 if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
314 }
315 ldrhs = RHS->Stride();
316 rhs = (*RHS)[0]; // Dangerous but effective (again)
317 ldlhs = LHS->Stride();
318 lhs = (*LHS)[0];
319
320 // Finally build ptr vector
321
322 if (ptr_==0) {
323 ptr_ = new int[M+1];
324 ptr_[0] = 0;
325 for (int i=0; i<M; i++) ptr_[i+1] = ptr_[i] + RedistGraph.NumMyIndices(i);
326 }
327 ptr = ptr_;
328
329 return(0);
330}
331
332//=========================================================================
342
343//=========================================================================
353
354
355
@ Average
#define EPETRA_CHK_ERR(a)
int MyGlobalElements(int *MyGlobalElementList) const
Puts list of global elements on this processor into the user-provided array.
bool LinearMap() const
Returns true if the global ID space is contiguously divided (but not necessarily uniformly) across al...
int IndexBase() const
Index base for this map.
long long NumGlobalElements64() const
const Epetra_Comm & Comm() const
Access function for Epetra_Comm communicator.
int NumMyElements() const
Number of elements on the calling processor.
Epetra_Comm: The Epetra Communication Abstract Base Class.
Definition Epetra_Comm.h:73
virtual int NumProc() const =0
Returns total number of processes.
virtual int MyPID() const =0
Return my process ID.
Epetra_CrsGraph: A class for constructing and using sparse compressed row graphs.
int NumMyIndices(int Row) const
Returns the current number of nonzero entries in specified local row on this processor.
Epetra_CrsMatrix: A class for constructing and using real-valued double-precision sparse compressed r...
int NumMyCols() const
Returns the number of entries in the set of column-indices that appear on this processor.
int FillComplete(bool OptimizeDataStorage=true)
Signal that data entry is complete. Perform transformations to local index space.
int PutScalar(double ScalarConstant)
Initialize all values in the matrix with constant value.
int NumMyNonzeros() const
Returns the number of nonzero entries in the calling processor's portion of the matrix.
const Epetra_CrsGraph & Graph() const
Returns a reference to the Epetra_CrsGraph object associated with this matrix.
int NumMyRows() const
Returns the number of matrix rows owned by the calling processor.
bool IndicesAreContiguous() const
If matrix indices are packed into single array (done in OptimizeStorage()) return true,...
int Import(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Imports an Epetra_DistObject using the Epetra_Import object.
int Export(const Epetra_SrcDistObject &A, const Epetra_Import &Importer, Epetra_CombineMode CombineMode, const Epetra_OffsetIndex *Indexor=0)
Exports an Epetra_DistObject using the Epetra_Import object.
Epetra_Export: This class builds an export object for efficient exporting of off-processor elements.
Epetra_Import: This class builds an import object for efficient importing of off-processor elements.
Epetra_IntVector: A class for constructing and using dense integer vectors on a parallel computer.
int * Values() const
Returns a pointer to an array containing the values of this vector.
Epetra_LinearProblemRedistor: A class for redistributing an Epetra_LinearProblem object.
int UpdateRedistRHS(Epetra_MultiVector *RHSWithNewValues)
Update the values of an already-redistributed RHS.
int CreateRedistProblem(const bool ConstructTranspose, const bool MakeDataContiguous, Epetra_LinearProblem *&RedistProblem)
Generate a new Epetra_LinearProblem as a redistribution of the one passed into the constructor.
virtual ~Epetra_LinearProblemRedistor()
Epetra_LinearProblemRedistor destructor.
Epetra_LinearProblemRedistor(Epetra_LinearProblem *OrigProblem, const Epetra_Map &RedistMap)
Epetra_LinearProblemRedistor constructor using pre-defined layout.
int UpdateRedistProblemValues(Epetra_LinearProblem *ProblemWithNewValues)
Update the values of an already-redistributed problem.
Epetra_RowMatrixTransposer * Transposer_
int UpdateOriginalLHS(Epetra_MultiVector *LHS)
Update LHS of original Linear Problem object.
int ExtractHbData(int &M, int &N, int &nz, int *&ptr, int *&ind, double *&val, int &Nrhs, double *&rhs, int &ldrhs, double *&lhs, int &ldlhs) const
Extract the redistributed problem data in a form usable for other codes that require Harwell-Boeing f...
Epetra_LinearProblem: The Epetra Linear Problem Class.
Epetra_MultiVector * GetRHS() const
Get a pointer to the right-hand-side B.
void SetOperator(Epetra_RowMatrix *A)
Set Operator A of linear problem AX = B using an Epetra_RowMatrix.
void SetRHS(Epetra_MultiVector *B)
Set right-hand-side B of linear problem AX = B.
Epetra_RowMatrix * GetMatrix() const
Get a pointer to the matrix A.
void SetLHS(Epetra_MultiVector *X)
Set left-hand-side X of linear problem AX = B.
Epetra_MultiVector * GetLHS() const
Get a pointer to the left-hand-side X.
Epetra_Map: A class for partitioning vectors and matrices.
Definition Epetra_Map.h:119
Epetra_MultiVector: A class for constructing and using dense multi-vectors, vectors and matrices in p...
int NumVectors() const
Returns the number of vectors in the multi-vector.
bool ConstantStride() const
Returns true if this multi-vector has constant stride between vectors.
int Stride() const
Returns the stride between vectors in the multi-vector (only meaningful if ConstantStride() is true).
Epetra_RowMatrixTransposer: A class for transposing an Epetra_RowMatrix object.
int UpdateTransposeValues(Epetra_RowMatrix *MatrixWithNewValues)
Update the values of an already-redistributed problem.
int CreateTranspose(const bool MakeDataContiguous, Epetra_CrsMatrix *&TransposeMatrix, Epetra_Map *TransposeRowMap=0)
Generate a new Epetra_CrsMatrix as the transpose of an Epetra_RowMatrix passed into the constructor.
Epetra_RowMatrix: A pure virtual class for using real-valued double-precision row matrices.
virtual const Epetra_Map & RowMatrixRowMap() const =0
Returns the Epetra_Map object associated with the rows of this matrix.