Epetra Package Browser (Single Doxygen Collection) Development
Epetra_RowMatrixTransposer.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_RowMatrixTransposer.h"
00044 #include "Epetra_RowMatrix.h"
00045 #include "Epetra_CrsMatrix.h"
00046 #include "Epetra_CrsGraph.h"
00047 #include "Epetra_Map.h"
00048 #include "Epetra_Export.h"
00049 //=============================================================================
00050 Epetra_RowMatrixTransposer::Epetra_RowMatrixTransposer(Epetra_RowMatrix * OrigMatrix)
00051   : OrigMatrix_(OrigMatrix),
00052     TransposeMatrix_(0),
00053     TransposeExporter_(0),
00054     TransposeRowMap_(0),
00055     TransposeCreated_(false),
00056     MakeDataContiguous_(false),
00057     NumMyRows_(0),
00058     NumMyCols_(0),
00059     MaxNumEntries_(0),
00060     Indices_(NULL),
00061     Values_(NULL),
00062     TransNumNz_(NULL),
00063     TransIndices_(NULL),
00064     TransValues_(NULL),
00065     TransMyGlobalEquations_(NULL),
00066     OrigMatrixIsCrsMatrix_(false)
00067 {
00068 }
00069 //=============================================================================
00070 Epetra_RowMatrixTransposer::Epetra_RowMatrixTransposer(const Epetra_RowMatrixTransposer& Source)
00071   :OrigMatrix_(Source.OrigMatrix_),
00072    TransposeMatrix_(0),
00073    TransposeExporter_(0),
00074    TransposeRowMap_(0),
00075    TransposeCreated_(Source.TransposeCreated_),
00076    MakeDataContiguous_(Source.MakeDataContiguous_),
00077    NumMyRows_(0),
00078    NumMyCols_(0),
00079    MaxNumEntries_(0),
00080    Indices_(NULL),
00081    Values_(NULL),
00082    TransNumNz_(NULL),
00083    TransIndices_(NULL),
00084    TransValues_(NULL),
00085    TransMyGlobalEquations_(NULL),
00086    OrigMatrixIsCrsMatrix_(false)
00087 {
00088   TransposeMatrix_ = new Epetra_CrsMatrix(*Source.TransposeMatrix_);
00089   if (MakeDataContiguous_) TransposeMatrix_->MakeDataContiguous();
00090   TransposeExporter_ = new Epetra_Export(*Source.TransposeExporter_);
00091 }
00092 //=========================================================================
00093 Epetra_RowMatrixTransposer::~Epetra_RowMatrixTransposer(){
00094 
00095   DeleteData();
00096 
00097 }
00098 
00099 //=========================================================================
00100 void Epetra_RowMatrixTransposer::DeleteData (){
00101 
00102   int i;
00103 
00104   if (TransposeExporter_!=0) {delete TransposeExporter_; TransposeExporter_=0;}
00105 
00106   // Delete any intermediate storage
00107 
00108   if (!OrigMatrixIsCrsMatrix_) {
00109     delete [] Indices_;
00110     delete [] Values_;
00111   }
00112   
00113   
00114   for(i=0; i<NumMyCols_; i++) {
00115     int NumIndices = TransNumNz_[i];
00116     if (NumIndices>0) {
00117       delete [] TransIndices_[i];
00118       delete [] TransValues_[i];
00119     }
00120   }
00121   delete [] TransNumNz_;
00122   delete [] TransIndices_;
00123   delete [] TransValues_;
00124   delete [] TransMyGlobalEquations_;
00125 }
00126 
00127 //=========================================================================
00128 int Epetra_RowMatrixTransposer::CreateTranspose (const bool MakeDataContiguous, 
00129              Epetra_CrsMatrix *& TransposeMatrix, 
00130              Epetra_Map * TransposeRowMap_in) {
00131 
00132   int i, j;
00133 
00134   if (TransposeCreated_) DeleteData(); // Get rid of existing data first
00135 
00136   if (TransposeRowMap_in==0)
00137     TransposeRowMap_ = (Epetra_Map *) &(OrigMatrix_->OperatorDomainMap()); // Should be replaced with refcount =
00138   else
00139     TransposeRowMap_ = TransposeRowMap_in; 
00140 
00141   // This routine will work for any RowMatrix object, but will attempt cast the matrix to a CrsMatrix if
00142   // possible (because we can then use a View of the matrix and graph, which is much cheaper).
00143 
00144   // First get the local indices to count how many nonzeros will be in the 
00145   // transpose graph on each processor
00146 
00147 
00148   Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix_);
00149 
00150   OrigMatrixIsCrsMatrix_ = (OrigCrsMatrix!=0); // If this pointer is non-zero, the cast to CrsMatrix worked
00151 
00152   NumMyRows_ = OrigMatrix_->NumMyRows();
00153   NumMyCols_ = OrigMatrix_->NumMyCols();
00154   NumMyRows_ = OrigMatrix_->NumMyRows();
00155   TransNumNz_ = new int[NumMyCols_];
00156   TransIndices_ = new int*[NumMyCols_];
00157   TransValues_ = new double*[NumMyCols_];
00158 
00159 
00160   int NumIndices;
00161 
00162   if (OrigMatrixIsCrsMatrix_) {
00163 
00164 
00165     const Epetra_CrsGraph & OrigGraph = OrigCrsMatrix->Graph(); // Get matrix graph
00166 
00167     for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0;
00168     for (i=0; i<NumMyRows_; i++) {
00169       EPETRA_CHK_ERR(OrigGraph.ExtractMyRowView(i, NumIndices, Indices_)); // Get view of ith row
00170       for (j=0; j<NumIndices; j++) ++TransNumNz_[Indices_[j]];
00171     }
00172   }
00173   else { // OrigMatrix is not a CrsMatrix
00174 
00175     MaxNumEntries_ = 0;
00176     int NumEntries;
00177     for (i=0; i<NumMyRows_; i++) {
00178       OrigMatrix_->NumMyRowEntries(i, NumEntries);
00179       MaxNumEntries_ = EPETRA_MAX(MaxNumEntries_, NumEntries);
00180     }
00181     Indices_ = new int[MaxNumEntries_];
00182     Values_ = new double[MaxNumEntries_];
00183 
00184     for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0;
00185     for (i=0; i<NumMyRows_; i++) {
00186       // Get ith row
00187       EPETRA_CHK_ERR(OrigMatrix_->ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_)); 
00188       for (j=0; j<NumIndices; j++) ++TransNumNz_[Indices_[j]];
00189     }
00190   }
00191 
00192 
00193   // Most of remaining code is common to both cases
00194 
00195   for(i=0; i<NumMyCols_; i++) {
00196     NumIndices = TransNumNz_[i];
00197     if (NumIndices>0) {
00198       TransIndices_[i] = new int[NumIndices];
00199       TransValues_[i] = new double[NumIndices];
00200     }
00201   }
00202 
00203   // Now copy values and global indices into newly created transpose storage
00204 
00205   for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0; // Reset transpose NumNz counter
00206   for (i=0; i<NumMyRows_; i++) {
00207     if (OrigMatrixIsCrsMatrix_) {
00208       EPETRA_CHK_ERR(OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_));
00209     }
00210     else {
00211       EPETRA_CHK_ERR(OrigMatrix_->ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_));
00212     }
00213 
00214     int ii = OrigMatrix_->RowMatrixRowMap().GID(i);
00215     for (j=0; j<NumIndices; j++) {
00216       int TransRow = Indices_[j];
00217       int loc = TransNumNz_[TransRow];
00218       TransIndices_[TransRow][loc] = ii;
00219       TransValues_[TransRow][loc] = Values_[j];
00220       ++TransNumNz_[TransRow]; // increment counter into current transpose row
00221     }
00222   }
00223 
00224   //  Build Transpose matrix with some rows being shared across processors.
00225   //  We will use a view here since the matrix will not be used for anything else
00226 
00227   const Epetra_Map & TransMap = OrigMatrix_->RowMatrixColMap();
00228 
00229   Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz_);
00230   TransMyGlobalEquations_ = new int[NumMyCols_];
00231   TransMap.MyGlobalElements(TransMyGlobalEquations_);
00232   
00233   /* Add  rows one-at-a-time */
00234 
00235   for (i=0; i<NumMyCols_; i++)
00236     {
00237       EPETRA_CHK_ERR(TempTransA1.InsertGlobalValues(TransMyGlobalEquations_[i], 
00238                 TransNumNz_[i], TransValues_[i], TransIndices_[i]));
00239     }
00240  
00241   // Note: The following call to FillComplete is currently necessary because
00242   //       some global constants that are needed by the Export () are computed in this routine
00243 
00244   const Epetra_Map& domain_map = OrigMatrix_->OperatorDomainMap();
00245   const Epetra_Map& range_map = OrigMatrix_->OperatorRangeMap();
00246 
00247   EPETRA_CHK_ERR(TempTransA1.FillComplete(range_map, domain_map, false));
00248 
00249   // Now that transpose matrix with shared rows is entered, create a new matrix that will
00250   // get the transpose with uniquely owned rows (using the same row distribution as A).
00251 
00252   TransposeMatrix_ = new Epetra_CrsMatrix(Copy, *TransposeRowMap_,0);
00253 
00254   // Create an Export object that will move TempTransA around
00255 
00256   TransposeExporter_ = new Epetra_Export(TransMap, *TransposeRowMap_);
00257 
00258   EPETRA_CHK_ERR(TransposeMatrix_->Export(TempTransA1, *TransposeExporter_, Add));
00259   
00260   EPETRA_CHK_ERR(TransposeMatrix_->FillComplete(range_map, domain_map));
00261 
00262   if (MakeDataContiguous) {
00263     EPETRA_CHK_ERR(TransposeMatrix_->MakeDataContiguous());
00264   }
00265 
00266   TransposeMatrix = TransposeMatrix_;
00267   TransposeCreated_ = true;
00268 
00269   return(0);
00270 }
00271 //=========================================================================
00272 int Epetra_RowMatrixTransposer::UpdateTransposeValues(Epetra_RowMatrix * MatrixWithNewValues){
00273 
00274   int i, j, NumIndices;
00275 
00276   if (!TransposeCreated_) EPETRA_CHK_ERR(-1); // Transpose must be already created
00277 
00278   // Sanity check of incoming matrix.  Perform some tests to see if it is compatible with original input matrix
00279   if (OrigMatrix_!=MatrixWithNewValues) { // Check if pointer of new matrix is same as previous input matrix
00280     OrigMatrix_ = MatrixWithNewValues; // Reset this pointer if not, then check for other attributes
00281     if (NumMyRows_ != OrigMatrix_->NumMyRows() ||
00282   NumMyCols_ != OrigMatrix_->NumMyCols() ||
00283   NumMyRows_ != OrigMatrix_->NumMyRows()) {
00284       EPETRA_CHK_ERR(-2); // New matrix not compatible with previous
00285     }
00286   }
00287 
00288   Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(MatrixWithNewValues);
00289 
00290   
00291   OrigMatrixIsCrsMatrix_ = (OrigCrsMatrix!=0); // If this pointer is non-zero, the cast to CrsMatrix worked
00292 
00293 
00294   // Now copy values and global indices into newly create transpose storage
00295 
00296   for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0; // Reset transpose NumNz counter
00297   for (i=0; i<NumMyRows_; i++) {
00298     if (OrigMatrixIsCrsMatrix_) {
00299       EPETRA_CHK_ERR(OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_));
00300     }
00301     else {
00302       EPETRA_CHK_ERR(OrigMatrix_->ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_));
00303     }
00304 
00305     int ii = OrigMatrix_->RowMatrixRowMap().GID(i);
00306     for (j=0; j<NumIndices; j++) {
00307       int TransRow = Indices_[j];
00308       int loc = TransNumNz_[TransRow];
00309       TransIndices_[TransRow][loc] = ii;
00310       TransValues_[TransRow][loc] = Values_[j];
00311       ++TransNumNz_[TransRow]; // increment counter into current transpose row
00312     }
00313   }
00314 
00315   //  Build Transpose matrix with some rows being shared across processors.
00316   //  We will use a view here since the matrix will not be used for anything else
00317 
00318   const Epetra_Map & TransMap = OrigMatrix_->RowMatrixColMap();
00319 
00320   Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz_);
00321   TransMap.MyGlobalElements(TransMyGlobalEquations_);
00322   
00323   /* Add  rows one-at-a-time */
00324 
00325   for (i=0; i<NumMyCols_; i++)
00326     {
00327       EPETRA_CHK_ERR(TempTransA1.InsertGlobalValues(TransMyGlobalEquations_[i], 
00328                 TransNumNz_[i], TransValues_[i], TransIndices_[i]));
00329     }
00330  
00331   // Note: The following call to FillComplete is currently necessary because
00332   //       some global constants that are needed by the Export () are computed in this routine
00333   const Epetra_Map& domain_map = OrigMatrix_->OperatorDomainMap();
00334   const Epetra_Map& range_map = OrigMatrix_->OperatorRangeMap();
00335 
00336   EPETRA_CHK_ERR(TempTransA1.FillComplete(range_map, domain_map, false));
00337 
00338   // Now that transpose matrix with shared rows is entered, update values of target transpose matrix
00339 
00340   TransposeMatrix_->PutScalar(0.0);  // Zero out all values of the matrix
00341 
00342   EPETRA_CHK_ERR(TransposeMatrix_->Export(TempTransA1, *TransposeExporter_, Add));
00343 
00344   return(0);
00345 }
00346 
00347 Epetra_RowMatrixTransposer&
00348 Epetra_RowMatrixTransposer::operator=(const Epetra_RowMatrixTransposer& src)
00349 {
00350   (void)src;//prevents unused variable compiler warning
00351 
00352   //not currently supported
00353   bool throw_error = true;
00354   if (throw_error) {
00355     std::cerr << std::endl
00356         << "Epetra_RowMatrixTransposer::operator= not supported."
00357         <<std::endl;
00358     throw -1;
00359   }
00360 
00361   return(*this);
00362 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines