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

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