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