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