Epetra_LinearProblemRedistor.cpp

Go to the documentation of this file.
00001 
00002 //@HEADER
00003 // ************************************************************************
00004 // 
00005 //               Epetra: Linear Algebra Services Package 
00006 //                 Copyright (2001) Sandia Corporation
00007 // 
00008 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00009 // license for use of this work by or on behalf of the U.S. Government.
00010 // 
00011 // This library is free software; you can redistribute it and/or modify
00012 // it under the terms of the GNU Lesser General Public License as
00013 // published by the Free Software Foundation; either version 2.1 of the
00014 // License, or (at your option) any later version.
00015 //  
00016 // This library is distributed in the hope that it will be useful, but
00017 // WITHOUT ANY WARRANTY; without even the implied warranty of
00018 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00019 // Lesser General Public License for more details.
00020 //  
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License along with this library; if not, write to the Free Software
00023 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00024 // USA
00025 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00026 // 
00027 // ************************************************************************
00028 //@HEADER
00029 
00030 #include "Epetra_Comm.h"
00031 #include "Epetra_Map.h"
00032 #include "Epetra_RowMatrixTransposer.h"
00033 #include "Epetra_CrsMatrix.h"
00034 #include "Epetra_LinearProblem.h"
00035 #include "Epetra_LinearProblemRedistor.h"
00036 #include "Epetra_Util.h"
00037 #include "Epetra_MultiVector.h"
00038 #include "Epetra_IntVector.h"
00039 #include "Epetra_Export.h"
00040 #include "Epetra_Import.h"
00041 //=============================================================================
00042 Epetra_LinearProblemRedistor::Epetra_LinearProblemRedistor(Epetra_LinearProblem * OrigProblem, 
00043                                                            const Epetra_Map & RedistMap)
00044   : OrigProblem_(OrigProblem),
00045     NumProc_(0),
00046     RedistProblem_(0),
00047     RedistMap_((Epetra_Map *) &RedistMap),
00048     Transposer_(0),
00049     Replicate_(false),
00050     ConstructTranspose_(false),
00051     MakeDataContiguous_(false),
00052     MapGenerated_(false),
00053     RedistProblemCreated_(false),
00054     ptr_(0)
00055 {
00056 }
00057 //=============================================================================
00058 Epetra_LinearProblemRedistor::Epetra_LinearProblemRedistor(Epetra_LinearProblem * OrigProblem, 
00059                                                            int NumProc,
00060                                                            bool Replicate)
00061   : OrigProblem_(OrigProblem),
00062     NumProc_(NumProc),
00063     RedistProblem_(0),
00064     RedistMap_(0),
00065     Transposer_(0),
00066     Replicate_(Replicate),
00067     ConstructTranspose_(false),
00068     MakeDataContiguous_(false),
00069     MapGenerated_(false),
00070     RedistProblemCreated_(false),
00071     ptr_(0)
00072 {
00073 }
00074 //=============================================================================
00075 Epetra_LinearProblemRedistor::Epetra_LinearProblemRedistor(const Epetra_LinearProblemRedistor& Source)
00076   : OrigProblem_(Source.OrigProblem_),
00077     NumProc_(Source.NumProc_),
00078     RedistProblem_(Source.RedistProblem_),
00079     RedistMap_(Source.RedistMap_),
00080     Transposer_(Source.Transposer_),
00081     Replicate_(Source.Replicate_),
00082     ConstructTranspose_(Source.ConstructTranspose_),
00083     MakeDataContiguous_(Source.MakeDataContiguous_),
00084     RedistProblemCreated_(Source.RedistProblemCreated_),
00085     ptr_(0)
00086 {
00087 
00088   if (RedistProblem_!=0) RedistProblem_ = new Epetra_LinearProblem(*Source.RedistProblem_);
00089   if (Transposer_!=0) Transposer_ = new Epetra_RowMatrixTransposer(*Source.Transposer_);
00090 }
00091 //=========================================================================
00092 Epetra_LinearProblemRedistor::~Epetra_LinearProblemRedistor(){
00093 
00094   if (ptr_!=0) {delete [] ptr_; ptr_=0;}
00095   if (MapGenerated_ && RedistMap_!=0) {delete RedistMap_; RedistMap_=0;}
00096   if (RedistProblem_!=0) {
00097      // If no tranpose, then we must delete matrix (otherwise the transposer must).
00098     if (!ConstructTranspose_ && RedistProblem_->GetMatrix()!=0)
00099       delete RedistProblem_->GetMatrix();
00100     if (RedistProblem_->GetLHS()!=0) delete RedistProblem_->GetLHS();
00101     if (RedistProblem_->GetRHS()!=0) delete RedistProblem_->GetRHS();
00102     delete RedistProblem_; RedistProblem_=0;
00103   }
00104   if (RedistExporter_!=0) {delete RedistExporter_; RedistExporter_=0;}
00105   if (Transposer_!=0) {delete Transposer_; Transposer_ = 0;}
00106 
00107 }
00108 
00109 //=========================================================================
00110 int Epetra_LinearProblemRedistor::GenerateRedistMap() {
00111 
00112 
00113   if (MapGenerated_) return(0);
00114 
00115   const Epetra_Map & SourceMap = OrigProblem_->GetMatrix()->RowMatrixRowMap();
00116   const Epetra_Comm & Comm = SourceMap.Comm();
00117   int IndexBase = SourceMap.IndexBase();
00118 
00119   int NumProc = Comm.NumProc();
00120 
00121   if (NumProc_!=NumProc) return(-1); // Right now we are only supporting redistribution to all processors.
00122 
00123   // Build a list of contiguous GIDs that will be used in either case below.
00124   int NumMyRedistElements = 0;
00125   if ((Comm.MyPID()==0) || Replicate_) NumMyRedistElements = SourceMap.NumGlobalElements();
00126   
00127   // Now build the GID list to broadcast SourceMapGIDs to all processors that need it (or just to PE 0).
00128   int * ContigIDs = 0;
00129   if (NumMyRedistElements>0) ContigIDs = new int[NumMyRedistElements];
00130   for (int i=0; i<NumMyRedistElements; i++) ContigIDs[i] = IndexBase + i;
00131   
00132   // Case 1: If the map for the input matrix is not a linear contiguous map, then we have to collect
00133   // the map indices in order to construct a compatible map containing all GIDs.
00134   if (!SourceMap.LinearMap()) {
00135 
00136     // First generate a linear map of the same distribution as RowMatrixRowMap
00137     Epetra_Map SourceLinearMap(-1, SourceMap.NumMyElements(), IndexBase, Comm);
00138     // Generate Int vector containing GIDs of SourceMap
00139     Epetra_IntVector SourceMapGIDs(View, SourceLinearMap, SourceMap.MyGlobalElements());
00140     // Now Build target map for SourceMapGIDs and Importer to get IDs
00141     Epetra_Map GIDsTargetMap(-1, NumMyRedistElements, ContigIDs, IndexBase, Comm);
00142     if (NumMyRedistElements>0) delete [] ContigIDs;
00143 
00144     Epetra_Import GIDsImporter(GIDsTargetMap, SourceMap);
00145     Epetra_IntVector TargetMapGIDs(GIDsTargetMap);
00146 
00147     // Now Send SourceMap GIDs to PE 0, and all processors if Replicate is true
00148     EPETRA_CHK_ERR(TargetMapGIDs.Import(SourceMapGIDs, GIDsImporter, Insert));
00149 
00150     // Finally, create RedistMap containing all GIDs of SourceMap on PE 0, or all PEs of Replicate is true
00151     RedistMap_ = new Epetra_Map(-1, NumMyRedistElements, TargetMapGIDs.Values(), IndexBase, Comm);
00152 
00153   }
00154   // Case 2: If the map has contiguous IDs then we can simply build the map right away using the list
00155   // of contiguous GIDs directly
00156   else
00157     RedistMap_ = new Epetra_Map(-1, NumMyRedistElements, ContigIDs, IndexBase, Comm);
00158 
00159   MapGenerated_ = true;
00160 
00161   return(0);
00162 }
00163 
00164 //=========================================================================
00165 int Epetra_LinearProblemRedistor::CreateRedistProblem(const bool ConstructTranspose, 
00166                                                       const bool MakeDataContiguous, 
00167                                                       Epetra_LinearProblem *& RedistProblem) {
00168 
00169   if (RedistProblemCreated_) EPETRA_CHK_ERR(-1);  // This method can only be called once
00170 
00171   Epetra_RowMatrix * OrigMatrix = OrigProblem_->GetMatrix();
00172   Epetra_MultiVector * OrigLHS = OrigProblem_->GetLHS();
00173   Epetra_MultiVector * OrigRHS = OrigProblem_->GetRHS();
00174     
00175   if (OrigMatrix==0) EPETRA_CHK_ERR(-2); // There is no matrix associated with this Problem
00176 
00177 
00178   if (RedistMap_==0) {
00179     EPETRA_CHK_ERR(GenerateRedistMap());
00180   }
00181 
00182   RedistExporter_ = new Epetra_Export(OrigProblem_->GetMatrix()->RowMatrixRowMap(), *RedistMap_);
00183 
00184   RedistProblem_ = new Epetra_LinearProblem();
00185   Epetra_CrsMatrix * RedistMatrix;
00186 
00187   // Check if the tranpose should be create or not
00188   if (ConstructTranspose) {
00189     Transposer_ = new Epetra_RowMatrixTransposer(OrigMatrix);
00190     EPETRA_CHK_ERR(Transposer_->CreateTranspose(MakeDataContiguous, RedistMatrix, RedistMap_));
00191   }
00192   else {
00193     // If not, then just do the redistribution based on the the RedistMap
00194     RedistMatrix = new Epetra_CrsMatrix(Copy, *RedistMap_, 0);
00195     // need to do this next step until we generalize the Import/Export ops for CrsMatrix
00196     Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix); 
00197     EPETRA_CHK_ERR(RedistMatrix->Export(*OrigCrsMatrix, *RedistExporter_, Add));
00198     EPETRA_CHK_ERR(RedistMatrix->FillComplete());
00199   }
00200 
00201   RedistProblem_->SetOperator(RedistMatrix);
00202 
00203   // Now redistribute the RHS and LHS if non-zero
00204 
00205   Epetra_MultiVector * RedistLHS = 0;
00206   Epetra_MultiVector * RedistRHS = 0;
00207 
00208   int ierr = 0;
00209 
00210   if (OrigLHS!=0) {
00211     RedistLHS = new Epetra_MultiVector(*RedistMap_, OrigLHS->NumVectors());
00212     EPETRA_CHK_ERR(RedistLHS->Export(*OrigLHS, *RedistExporter_, Add));
00213   }
00214   else ierr = 1;
00215 
00216   if (OrigRHS!=0) {
00217     RedistRHS = new Epetra_MultiVector(*RedistMap_, OrigLHS->NumVectors());
00218     EPETRA_CHK_ERR(RedistRHS->Export(*OrigRHS, *RedistExporter_, Add));
00219   }
00220   else ierr ++;
00221 
00222   RedistProblem_->SetLHS(RedistLHS);
00223   RedistProblem_->SetRHS(RedistRHS);
00224 
00225   RedistProblemCreated_ = true;
00226 
00227   return(ierr);
00228 }
00229 
00230 //=========================================================================
00231 int Epetra_LinearProblemRedistor::UpdateRedistProblemValues(Epetra_LinearProblem * ProblemWithNewValues) {
00232   
00233   if (!RedistProblemCreated_) EPETRA_CHK_ERR(-1);  // This method can only be called after CreateRedistProblem()
00234 
00235   Epetra_RowMatrix * OrigMatrix = ProblemWithNewValues->GetMatrix();
00236   Epetra_MultiVector * OrigLHS = ProblemWithNewValues->GetLHS();
00237   Epetra_MultiVector * OrigRHS = ProblemWithNewValues->GetRHS();
00238     
00239   if (OrigMatrix==0) EPETRA_CHK_ERR(-2); // There is no matrix associated with this Problem
00240 
00241 
00242   Epetra_CrsMatrix * RedistMatrix = dynamic_cast<Epetra_CrsMatrix *>(RedistProblem_->GetMatrix());
00243 
00244   // Check if the tranpose should be create or not
00245   if (ConstructTranspose_) {
00246     EPETRA_CHK_ERR(Transposer_->UpdateTransposeValues(OrigMatrix));
00247   }
00248   else {
00249     // If not, then just do the redistribution based on the the RedistMap
00250     
00251     EPETRA_CHK_ERR(RedistMatrix->PutScalar(0.0));
00252     // need to do this next step until we generalize the Import/Export ops for CrsMatrix
00253     Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix); 
00254 
00255     if (OrigCrsMatrix==0) EPETRA_CHK_ERR(-3); // Broken for a RowMatrix at this point
00256     EPETRA_CHK_ERR(RedistMatrix->Export(*OrigCrsMatrix, *RedistExporter_, Add));
00257   }
00258 
00259   // Now redistribute the RHS and LHS if non-zero
00260 
00261 
00262   if (OrigLHS!=0) {
00263     EPETRA_CHK_ERR(RedistProblem_->GetLHS()->Export(*OrigLHS, *RedistExporter_, Add));
00264   }
00265 
00266   if (OrigRHS!=0) {
00267     EPETRA_CHK_ERR(RedistProblem_->GetRHS()->Export(*OrigRHS, *RedistExporter_, Add));
00268   }
00269 
00270   return(0);
00271 }
00272 
00273 //=========================================================================
00274 // NOTE: This method should be removed and replaced with calls to Epetra_Util_ExtractHbData()
00275 int Epetra_LinearProblemRedistor::ExtractHbData(int & M, int & N, int & nz, int * & ptr, 
00276                                                 int * & ind, double * & val, int & Nrhs, 
00277                                                 double * & rhs, int & ldrhs, 
00278                                                 double * & lhs, int & ldlhs) const {
00279 
00280   Epetra_CrsMatrix * RedistMatrix = dynamic_cast<Epetra_CrsMatrix *>(RedistProblem_->GetMatrix());
00281   
00282   if (RedistMatrix==0) EPETRA_CHK_ERR(-1); // This matrix is zero or not an Epetra_CrsMatrix
00283   if (!RedistMatrix->IndicesAreContiguous()) { // Data must be contiguous for this to work
00284     EPETRA_CHK_ERR(-2);
00285   }
00286 
00287   M = RedistMatrix->NumMyRows();
00288   N = RedistMatrix->NumMyCols();
00289   nz = RedistMatrix->NumMyNonzeros();
00290   val = (*RedistMatrix)[0];        // Dangerous, but cheap and effective way to access first element in 
00291 
00292   const Epetra_CrsGraph & RedistGraph = RedistMatrix->Graph();
00293   ind = RedistGraph[0];  // list of values and indices
00294 
00295   Epetra_MultiVector * LHS = RedistProblem_->GetLHS();
00296   Epetra_MultiVector * RHS = RedistProblem_->GetRHS();
00297   Nrhs = RHS->NumVectors();
00298   if (Nrhs>1) {
00299     if (!RHS->ConstantStride()) {EPETRA_CHK_ERR(-3)}; // Must have strided vectors
00300     if (!LHS->ConstantStride()) {EPETRA_CHK_ERR(-4)}; // Must have strided vectors
00301   }
00302   ldrhs = RHS->Stride();
00303   rhs = (*RHS)[0]; // Dangerous but effective (again)
00304   ldlhs = LHS->Stride();
00305   lhs = (*LHS)[0];
00306 
00307   // Finally build ptr vector
00308 
00309   if (ptr_==0) {
00310     ptr_ = new int[M+1];
00311     ptr_[0] = 0;
00312     for (int i=0; i<M; i++) ptr_[i+1] = ptr_[i] + RedistGraph.NumMyIndices(i);
00313   }
00314   ptr = ptr_;
00315   
00316   return(0);
00317 }
00318 
00319 //=========================================================================
00320 int Epetra_LinearProblemRedistor::UpdateOriginalLHS(Epetra_MultiVector * LHS) {
00321 
00322   if (LHS==0) {EPETRA_CHK_ERR(-1);}
00323   if (!RedistProblemCreated_) EPETRA_CHK_ERR(-2);
00324 
00325   EPETRA_CHK_ERR(LHS->Import(*RedistProblem_->GetLHS(), *RedistExporter_, Average));
00326 
00327   return(0);
00328 }
00329 
00330 //=========================================================================
00331 int Epetra_LinearProblemRedistor::UpdateRedistRHS(Epetra_MultiVector * RHS) {
00332 
00333   if (RHS==0) {EPETRA_CHK_ERR(-1);}
00334   if (!RedistProblemCreated_) EPETRA_CHK_ERR(-2);
00335 
00336   EPETRA_CHK_ERR(RedistProblem_->GetRHS()->Export(*RHS, *RedistExporter_, Insert));
00337 
00338   return(0);
00339 }
00340 
00341 
00342 

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7