Epetra Package Browser (Single Doxygen Collection) Development
Epetra_LinearProblemRedistor.cpp
Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright 2011 Sandia Corporation
00007 // 
00008 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00009 // the U.S. Government retains certain rights in this software.
00010 //
00011 // Redistribution and use in source and binary forms, with or without
00012 // modification, are permitted provided that the following conditions are
00013 // met:
00014 //
00015 // 1. Redistributions of source code must retain the above copyright
00016 // notice, this list of conditions and the following disclaimer.
00017 //
00018 // 2. Redistributions in binary form must reproduce the above copyright
00019 // notice, this list of conditions and the following disclaimer in the
00020 // documentation and/or other materials provided with the distribution.
00021 //
00022 // 3. Neither the name of the Corporation nor the names of the
00023 // contributors may be used to endorse or promote products derived from
00024 // this software without specific prior written permission.
00025 //
00026 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00027 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00028 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00029 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00030 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00031 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00032 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00033 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00034 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00035 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00036 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00037 //
00038 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00039 // 
00040 // ************************************************************************
00041 //@HEADER
00042 
00043 #include "Epetra_Comm.h"
00044 #include "Epetra_Map.h"
00045 #include "Epetra_RowMatrixTransposer.h"
00046 #include "Epetra_CrsMatrix.h"
00047 #include "Epetra_LinearProblem.h"
00048 #include "Epetra_LinearProblemRedistor.h"
00049 #include "Epetra_Util.h"
00050 #include "Epetra_MultiVector.h"
00051 #include "Epetra_IntVector.h"
00052 #include "Epetra_Export.h"
00053 #include "Epetra_Import.h"
00054 //=============================================================================
00055 Epetra_LinearProblemRedistor::Epetra_LinearProblemRedistor(Epetra_LinearProblem * OrigProblem, 
00056                                                            const Epetra_Map & RedistMap)
00057   : OrigProblem_(OrigProblem),
00058     NumProc_(0),
00059     RedistProblem_(0),
00060     RedistMap_((Epetra_Map *) &RedistMap),
00061     Transposer_(0),
00062     Replicate_(false),
00063     ConstructTranspose_(false),
00064     MakeDataContiguous_(false),
00065     MapGenerated_(false),
00066     RedistProblemCreated_(false),
00067     ptr_(0)
00068 {
00069 }
00070 //=============================================================================
00071 Epetra_LinearProblemRedistor::Epetra_LinearProblemRedistor(Epetra_LinearProblem * OrigProblem, 
00072                                                            int NumProc,
00073                                                            bool Replicate)
00074   : OrigProblem_(OrigProblem),
00075     NumProc_(NumProc),
00076     RedistProblem_(0),
00077     RedistMap_(0),
00078     Transposer_(0),
00079     Replicate_(Replicate),
00080     ConstructTranspose_(false),
00081     MakeDataContiguous_(false),
00082     MapGenerated_(false),
00083     RedistProblemCreated_(false),
00084     ptr_(0)
00085 {
00086 }
00087 //=============================================================================
00088 Epetra_LinearProblemRedistor::Epetra_LinearProblemRedistor(const Epetra_LinearProblemRedistor& Source)
00089   : OrigProblem_(Source.OrigProblem_),
00090     NumProc_(Source.NumProc_),
00091     RedistProblem_(Source.RedistProblem_),
00092     RedistMap_(Source.RedistMap_),
00093     Transposer_(Source.Transposer_),
00094     Replicate_(Source.Replicate_),
00095     ConstructTranspose_(Source.ConstructTranspose_),
00096     MakeDataContiguous_(Source.MakeDataContiguous_),
00097     RedistProblemCreated_(Source.RedistProblemCreated_),
00098     ptr_(0)
00099 {
00100 
00101   if (RedistProblem_!=0) RedistProblem_ = new Epetra_LinearProblem(*Source.RedistProblem_);
00102   if (Transposer_!=0) Transposer_ = new Epetra_RowMatrixTransposer(*Source.Transposer_);
00103 }
00104 //=========================================================================
00105 Epetra_LinearProblemRedistor::~Epetra_LinearProblemRedistor(){
00106 
00107   if (ptr_!=0) {delete [] ptr_; ptr_=0;}
00108   if (MapGenerated_ && RedistMap_!=0) {delete RedistMap_; RedistMap_=0;}
00109   if (RedistProblem_!=0) {
00110      // If no tranpose, then we must delete matrix (otherwise the transposer must).
00111     if (!ConstructTranspose_ && RedistProblem_->GetMatrix()!=0)
00112       delete RedistProblem_->GetMatrix();
00113     if (RedistProblem_->GetLHS()!=0) delete RedistProblem_->GetLHS();
00114     if (RedistProblem_->GetRHS()!=0) delete RedistProblem_->GetRHS();
00115     delete RedistProblem_; RedistProblem_=0;
00116   }
00117   if (RedistExporter_!=0) {delete RedistExporter_; RedistExporter_=0;}
00118   if (Transposer_!=0) {delete Transposer_; Transposer_ = 0;}
00119 
00120 }
00121 
00122 //=========================================================================
00123 int Epetra_LinearProblemRedistor::GenerateRedistMap() {
00124 
00125 
00126   if (MapGenerated_) return(0);
00127 
00128   const Epetra_Map & SourceMap = OrigProblem_->GetMatrix()->RowMatrixRowMap();
00129   const Epetra_Comm & Comm = SourceMap.Comm();
00130   int IndexBase = SourceMap.IndexBase();
00131 
00132   int NumProc = Comm.NumProc();
00133 
00134   if (NumProc_!=NumProc) return(-1); // Right now we are only supporting redistribution to all processors.
00135 
00136   // Build a list of contiguous GIDs that will be used in either case below.
00137   int NumMyRedistElements = 0;
00138   if ((Comm.MyPID()==0) || Replicate_) NumMyRedistElements = SourceMap.NumGlobalElements64();
00139   
00140   // Now build the GID list to broadcast SourceMapGIDs to all processors that need it (or just to PE 0).
00141   int * ContigIDs = 0;
00142   if (NumMyRedistElements>0) ContigIDs = new int[NumMyRedistElements];
00143   for (int i=0; i<NumMyRedistElements; i++) ContigIDs[i] = IndexBase + i;
00144   
00145   // Case 1: If the map for the input matrix is not a linear contiguous map, then we have to collect
00146   // the map indices in order to construct a compatible map containing all GIDs.
00147   if (!SourceMap.LinearMap()) {
00148 
00149     // First generate a linear map of the same distribution as RowMatrixRowMap
00150     Epetra_Map SourceLinearMap(-1, SourceMap.NumMyElements(), IndexBase, Comm);
00151     // Generate Int vector containing GIDs of SourceMap
00152     Epetra_IntVector SourceMapGIDs(View, SourceLinearMap, SourceMap.MyGlobalElements());
00153     // Now Build target map for SourceMapGIDs and Importer to get IDs
00154     Epetra_Map GIDsTargetMap(-1, NumMyRedistElements, ContigIDs, IndexBase, Comm);
00155     if (NumMyRedistElements>0) delete [] ContigIDs;
00156 
00157     Epetra_Import GIDsImporter(GIDsTargetMap, SourceMap);
00158     Epetra_IntVector TargetMapGIDs(GIDsTargetMap);
00159 
00160     // Now Send SourceMap GIDs to PE 0, and all processors if Replicate is true
00161     EPETRA_CHK_ERR(TargetMapGIDs.Import(SourceMapGIDs, GIDsImporter, Insert));
00162 
00163     // Finally, create RedistMap containing all GIDs of SourceMap on PE 0, or all PEs of Replicate is true
00164     RedistMap_ = new Epetra_Map(-1, NumMyRedistElements, TargetMapGIDs.Values(), IndexBase, Comm);// FIXME long long
00165 
00166   }
00167   // Case 2: If the map has contiguous IDs then we can simply build the map right away using the list
00168   // of contiguous GIDs directly
00169   else
00170     RedistMap_ = new Epetra_Map(-1, NumMyRedistElements, ContigIDs, IndexBase, Comm);// FIXME long long
00171 
00172   MapGenerated_ = true;
00173 
00174   return(0);
00175 }
00176 
00177 //=========================================================================
00178 int Epetra_LinearProblemRedistor::CreateRedistProblem(const bool ConstructTranspose, 
00179                                                       const bool MakeDataContiguous, 
00180                                                       Epetra_LinearProblem *& RedistProblem) {
00181 
00182   if (RedistProblemCreated_) EPETRA_CHK_ERR(-1);  // This method can only be called once
00183 
00184   Epetra_RowMatrix * OrigMatrix = OrigProblem_->GetMatrix();
00185   Epetra_MultiVector * OrigLHS = OrigProblem_->GetLHS();
00186   Epetra_MultiVector * OrigRHS = OrigProblem_->GetRHS();
00187     
00188   if (OrigMatrix==0) EPETRA_CHK_ERR(-2); // There is no matrix associated with this Problem
00189 
00190 
00191   if (RedistMap_==0) {
00192     EPETRA_CHK_ERR(GenerateRedistMap());
00193   }
00194 
00195   RedistExporter_ = new Epetra_Export(OrigProblem_->GetMatrix()->RowMatrixRowMap(), *RedistMap_);
00196 
00197   RedistProblem_ = new Epetra_LinearProblem();
00198   Epetra_CrsMatrix * RedistMatrix;
00199 
00200   // Check if the tranpose should be create or not
00201   if (ConstructTranspose) {
00202     Transposer_ = new Epetra_RowMatrixTransposer(OrigMatrix);
00203     EPETRA_CHK_ERR(Transposer_->CreateTranspose(MakeDataContiguous, RedistMatrix, RedistMap_));
00204   }
00205   else {
00206     // If not, then just do the redistribution based on the the RedistMap
00207     RedistMatrix = new Epetra_CrsMatrix(Copy, *RedistMap_, 0);
00208     // need to do this next step until we generalize the Import/Export ops for CrsMatrix
00209     Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix); 
00210     EPETRA_CHK_ERR(RedistMatrix->Export(*OrigCrsMatrix, *RedistExporter_, Add));
00211     EPETRA_CHK_ERR(RedistMatrix->FillComplete());
00212   }
00213 
00214   RedistProblem_->SetOperator(RedistMatrix);
00215 
00216   // Now redistribute the RHS and LHS if non-zero
00217 
00218   Epetra_MultiVector * RedistLHS = 0;
00219   Epetra_MultiVector * RedistRHS = 0;
00220 
00221   int ierr = 0;
00222 
00223   if (OrigLHS!=0) {
00224     RedistLHS = new Epetra_MultiVector(*RedistMap_, OrigLHS->NumVectors());
00225     EPETRA_CHK_ERR(RedistLHS->Export(*OrigLHS, *RedistExporter_, Add));
00226   }
00227   else ierr = 1;
00228 
00229   if (OrigRHS!=0) {
00230     RedistRHS = new Epetra_MultiVector(*RedistMap_, OrigLHS->NumVectors());
00231     EPETRA_CHK_ERR(RedistRHS->Export(*OrigRHS, *RedistExporter_, Add));
00232   }
00233   else ierr ++;
00234 
00235   RedistProblem_->SetLHS(RedistLHS);
00236   RedistProblem_->SetRHS(RedistRHS);
00237 
00238   RedistProblemCreated_ = true;
00239 
00240   return(ierr);
00241 }
00242 
00243 //=========================================================================
00244 int Epetra_LinearProblemRedistor::UpdateRedistProblemValues(Epetra_LinearProblem * ProblemWithNewValues) {
00245   
00246   if (!RedistProblemCreated_) EPETRA_CHK_ERR(-1);  // This method can only be called after CreateRedistProblem()
00247 
00248   Epetra_RowMatrix * OrigMatrix = ProblemWithNewValues->GetMatrix();
00249   Epetra_MultiVector * OrigLHS = ProblemWithNewValues->GetLHS();
00250   Epetra_MultiVector * OrigRHS = ProblemWithNewValues->GetRHS();
00251     
00252   if (OrigMatrix==0) EPETRA_CHK_ERR(-2); // There is no matrix associated with this Problem
00253 
00254 
00255   Epetra_CrsMatrix * RedistMatrix = dynamic_cast<Epetra_CrsMatrix *>(RedistProblem_->GetMatrix());
00256 
00257   // Check if the tranpose should be create or not
00258   if (ConstructTranspose_) {
00259     EPETRA_CHK_ERR(Transposer_->UpdateTransposeValues(OrigMatrix));
00260   }
00261   else {
00262     // If not, then just do the redistribution based on the the RedistMap
00263     
00264     EPETRA_CHK_ERR(RedistMatrix->PutScalar(0.0));
00265     // need to do this next step until we generalize the Import/Export ops for CrsMatrix
00266     Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix); 
00267 
00268     if (OrigCrsMatrix==0) EPETRA_CHK_ERR(-3); // Broken for a RowMatrix at this point
00269     EPETRA_CHK_ERR(RedistMatrix->Export(*OrigCrsMatrix, *RedistExporter_, Add));
00270   }
00271 
00272   // Now redistribute the RHS and LHS if non-zero
00273 
00274 
00275   if (OrigLHS!=0) {
00276     EPETRA_CHK_ERR(RedistProblem_->GetLHS()->Export(*OrigLHS, *RedistExporter_, Add));
00277   }
00278 
00279   if (OrigRHS!=0) {
00280     EPETRA_CHK_ERR(RedistProblem_->GetRHS()->Export(*OrigRHS, *RedistExporter_, Add));
00281   }
00282 
00283   return(0);
00284 }
00285 
00286 //=========================================================================
00287 // NOTE: This method should be removed and replaced with calls to Epetra_Util_ExtractHbData()
00288 int Epetra_LinearProblemRedistor::ExtractHbData(int & M, int & N, int & nz, int * & ptr, 
00289                                                 int * & ind, double * & val, int & Nrhs, 
00290                                                 double * & rhs, int & ldrhs, 
00291                                                 double * & lhs, int & ldlhs) const {
00292 
00293   Epetra_CrsMatrix * RedistMatrix = dynamic_cast<Epetra_CrsMatrix *>(RedistProblem_->GetMatrix());
00294   
00295   if (RedistMatrix==0) EPETRA_CHK_ERR(-1); // This matrix is zero or not an Epetra_CrsMatrix
00296   if (!RedistMatrix->IndicesAreContiguous()) { // Data must be contiguous for this to work
00297     EPETRA_CHK_ERR(-2);
00298   }
00299 
00300   M = RedistMatrix->NumMyRows();
00301   N = RedistMatrix->NumMyCols();
00302   nz = RedistMatrix->NumMyNonzeros();
00303   val = (*RedistMatrix)[0];        // Dangerous, but cheap and effective way to access first element in 
00304 
00305   const Epetra_CrsGraph & RedistGraph = RedistMatrix->Graph();
00306   ind = RedistGraph[0];  // list of values and indices
00307 
00308   Epetra_MultiVector * LHS = RedistProblem_->GetLHS();
00309   Epetra_MultiVector * RHS = RedistProblem_->GetRHS();
00310   Nrhs = RHS->NumVectors();
00311   if (Nrhs>1) {
00312     if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-3)}; // Must have strided vectors
00313     if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
00314   }
00315   ldrhs = RHS->Stride();
00316   rhs = (*RHS)[0]; // Dangerous but effective (again)
00317   ldlhs = LHS->Stride();
00318   lhs = (*LHS)[0];
00319 
00320   // Finally build ptr vector
00321 
00322   if (ptr_==0) {
00323     ptr_ = new int[M+1];
00324     ptr_[0] = 0;
00325     for (int i=0; i<M; i++) ptr_[i+1] = ptr_[i] + RedistGraph.NumMyIndices(i);
00326   }
00327   ptr = ptr_;
00328   
00329   return(0);
00330 }
00331 
00332 //=========================================================================
00333 int Epetra_LinearProblemRedistor::UpdateOriginalLHS(Epetra_MultiVector * LHS) {
00334 
00335   if (LHS==0) {EPETRA_CHK_ERR(-1);}
00336   if (!RedistProblemCreated_) EPETRA_CHK_ERR(-2);
00337 
00338   EPETRA_CHK_ERR(LHS->Import(*RedistProblem_->GetLHS(), *RedistExporter_, Average));
00339 
00340   return(0);
00341 }
00342 
00343 //=========================================================================
00344 int Epetra_LinearProblemRedistor::UpdateRedistRHS(Epetra_MultiVector * RHS) {
00345 
00346   if (RHS==0) {EPETRA_CHK_ERR(-1);}
00347   if (!RedistProblemCreated_) EPETRA_CHK_ERR(-2);
00348 
00349   EPETRA_CHK_ERR(RedistProblem_->GetRHS()->Export(*RHS, *RedistExporter_, Insert));
00350 
00351   return(0);
00352 }
00353 
00354 
00355 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines