EpetraExt Development
EpetraExt_Transpose_RowMatrix.cpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //     EpetraExt: Epetra Extended - Linear Algebra Services Package
00005 //                 Copyright (2001) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include <EpetraExt_Transpose_RowMatrix.h>
00030 
00031 #include <Epetra_Export.h>
00032 #include <Epetra_CrsGraph.h>
00033 #include <Epetra_CrsMatrix.h>
00034 #include <Epetra_Map.h>
00035 
00036 #include <vector>
00037 
00038 namespace EpetraExt {
00039 
00040 RowMatrix_Transpose::
00041 ~RowMatrix_Transpose()
00042 {
00043   if( TransposeMatrix_ ) delete TransposeMatrix_;
00044   if( TransposeExporter_ ) delete TransposeExporter_;
00045 
00046   if( !OrigMatrixIsCrsMatrix_ )
00047   {
00048     delete [] Indices_;
00049     delete [] Values_;
00050   }
00051 
00052   for( int i = 0; i < NumMyCols_; ++i )
00053     if( TransNumNz_[i] )
00054     {
00055       delete [] TransIndices_[i];
00056       delete [] TransValues_[i];
00057     }
00058 
00059   delete [] TransNumNz_;
00060   delete [] TransIndices_;
00061   delete [] TransValues_;
00062   delete [] TransMyGlobalEquations_;
00063 }
00064 
00065 //=========================================================================
00066 RowMatrix_Transpose::NewTypeRef
00067 RowMatrix_Transpose::
00068 operator()( OriginalTypeRef orig )
00069 {
00070   origObj_ = &orig;
00071 
00072   int i, j, err;
00073 
00074   if( !TransposeRowMap_ )
00075   {
00076     if( IgnoreNonLocalCols_ )
00077       TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorRangeMap()); // Should be replaced with refcount =
00078     else
00079       TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorDomainMap()); // Should be replaced with refcount =
00080   }
00081 
00082   // This routine will work for any RowMatrix object, but will attempt cast the matrix to a CrsMatrix if
00083   // possible (because we can then use a View of the matrix and graph, which is much cheaper).
00084 
00085   // First get the local indices to count how many nonzeros will be in the 
00086   // transpose graph on each processor
00087   Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix*>(&orig);
00088 
00089   OrigMatrixIsCrsMatrix_ = (OrigCrsMatrix!=0); // If this pointer is non-zero, the cast to CrsMatrix worked
00090 
00091   NumMyRows_ = orig.NumMyRows();
00092   NumMyCols_ = orig.NumMyCols();
00093   TransNumNz_ = new int[NumMyCols_];
00094   TransIndices_ = new int*[NumMyCols_];
00095   TransValues_ = new double*[NumMyCols_];
00096   TransMyGlobalEquations_ = new int[NumMyCols_];
00097 
00098   int NumIndices;
00099 
00100   if (OrigMatrixIsCrsMatrix_)
00101   {
00102     const Epetra_CrsGraph & OrigGraph = OrigCrsMatrix->Graph(); // Get matrix graph
00103 
00104     for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0;
00105     for (i=0; i<NumMyRows_; i++)
00106     {
00107       err = OrigGraph.ExtractMyRowView(i, NumIndices, Indices_); // Get view of ith row
00108       if (err != 0) throw OrigGraph.ReportError("ExtractMyRowView failed",err);
00109       for (j=0; j<NumIndices; j++) ++TransNumNz_[Indices_[j]];
00110     }
00111   }
00112   else // Original is not a CrsMatrix
00113   {
00114     MaxNumEntries_ = 0;
00115     int NumEntries;
00116     for (i=0; i<NumMyRows_; i++)
00117     {
00118       orig.NumMyRowEntries(i, NumEntries);
00119       MaxNumEntries_ = EPETRA_MAX(MaxNumEntries_, NumEntries);
00120     }
00121     Indices_ = new int[MaxNumEntries_];
00122     Values_ = new double[MaxNumEntries_];
00123 
00124     for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0;
00125     for (i=0; i<NumMyRows_; i++)
00126     {
00127       err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_); 
00128       if (err != 0) {
00129         std::cerr << "ExtractMyRowCopy failed."<<std::endl;
00130         throw err;
00131       }
00132       for (j=0; j<NumIndices; j++) ++TransNumNz_[Indices_[j]];
00133     }
00134   }
00135 
00136   // Most of remaining code is common to both cases
00137   for(i=0; i<NumMyCols_; i++)
00138   {
00139     NumIndices = TransNumNz_[i];
00140     if (NumIndices>0)
00141     {
00142       TransIndices_[i] = new int[NumIndices];
00143       TransValues_[i] = new double[NumIndices];
00144     }
00145   }
00146 
00147   // Now copy values and global indices into newly create transpose storage
00148 
00149   for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0; // Reset transpose NumNz counter
00150   for (i=0; i<NumMyRows_; i++)
00151   {
00152     if (OrigMatrixIsCrsMatrix_)
00153       err = OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_);
00154     else
00155       err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_);
00156     if (err != 0) {
00157       std::cerr << "ExtractMyRowCopy failed."<<std::endl;
00158       throw err;
00159     }
00160 
00161     int ii = orig.RowMatrixRowMap().GID(i);
00162     for (j=0; j<NumIndices; j++)
00163     {
00164       int TransRow = Indices_[j];
00165       int loc = TransNumNz_[TransRow];
00166       TransIndices_[TransRow][loc] = ii;
00167       TransValues_[TransRow][loc] = Values_[j];
00168       ++TransNumNz_[TransRow]; // increment counter into current transpose row
00169     }
00170   }
00171 
00172   //  Build Transpose matrix with some rows being shared across processors.
00173   //  We will use a view here since the matrix will not be used for anything else
00174 
00175   const Epetra_Map & TransMap = orig.RowMatrixColMap();
00176 
00177   Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz_);
00178   TransMap.MyGlobalElements(TransMyGlobalEquations_);
00179   
00180   for (i=0; i<NumMyCols_; i++) {
00181     err = TempTransA1.InsertGlobalValues(TransMyGlobalEquations_[i], 
00182                      TransNumNz_[i], TransValues_[i], TransIndices_[i]);
00183     if (err < 0) throw TempTransA1.ReportError("InsertGlobalValues failed.",err);
00184   }
00185  
00186   // Note: The following call to FillComplete is currently necessary because
00187   //      some global constants that are needed by the Export () are computed in this routine
00188   err = TempTransA1.FillComplete(orig.OperatorRangeMap(),*TransposeRowMap_, false);
00189   if (err != 0) {
00190     throw TempTransA1.ReportError("FillComplete failed.",err);
00191   }
00192 
00193   // Now that transpose matrix with shared rows is entered, create a new matrix that will
00194   // get the transpose with uniquely owned rows (using the same row distribution as A).
00195   if( IgnoreNonLocalCols_ )
00196     TransposeMatrix_ = new Epetra_CrsMatrix(Copy, *TransposeRowMap_, *TransposeRowMap_, 0);
00197   else
00198     TransposeMatrix_ = new Epetra_CrsMatrix(Copy, *TransposeRowMap_,0);
00199 
00200   // Create an Export object that will move TempTransA around
00201   TransposeExporter_ = new Epetra_Export(TransMap, *TransposeRowMap_);
00202 
00203   err = TransposeMatrix_->Export(TempTransA1, *TransposeExporter_, Add);
00204   if (err != 0) throw TransposeMatrix_->ReportError("Export failed.",err);
00205   
00206   err = TransposeMatrix_->FillComplete(orig.OperatorRangeMap(),*TransposeRowMap_);
00207   if (err != 0) throw TransposeMatrix_->ReportError("FillComplete failed.",err);
00208 
00209   if (MakeDataContiguous_) {
00210     err = TransposeMatrix_->MakeDataContiguous();
00211     if (err != 0) throw TransposeMatrix_->ReportError("MakeDataContiguous failed.",err);
00212   }
00213 
00214   newObj_ = TransposeMatrix_;
00215 
00216   return *newObj_;
00217 }
00218 
00219 //=========================================================================
00220 bool EpetraExt::RowMatrix_Transpose::fwd()
00221 {
00222   int i, j, NumIndices, err;
00223 
00224   Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix*>(origObj_);
00225 
00226   // Now copy values and global indices into newly create transpose storage
00227 
00228   for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0; // Reset transpose NumNz counter
00229   for (i=0; i<NumMyRows_; i++)
00230   {
00231     if(OrigMatrixIsCrsMatrix_)
00232       err = OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_);
00233     else
00234       err = origObj_->ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_);
00235     if (err != 0) {
00236       std::cerr << "ExtractMyRowCopy/View failed."<<std::endl;
00237       throw err;
00238     }
00239 
00240     int ii = origObj_->RowMatrixRowMap().GID(i);
00241     for (j=0; j<NumIndices; j++)
00242     {
00243       int TransRow = Indices_[j];
00244       int loc = TransNumNz_[TransRow];
00245       TransIndices_[TransRow][loc] = ii;
00246       TransValues_[TransRow][loc] = Values_[j];
00247       ++TransNumNz_[TransRow]; // increment counter into current transpose row
00248     }
00249   }
00250 
00251   //  Build Transpose matrix with some rows being shared across processors.
00252   //  We will use a view here since the matrix will not be used for anything else
00253   const Epetra_Map & TransMap = origObj_->RowMatrixColMap();
00254 
00255   Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz_);
00256   TransMap.MyGlobalElements(TransMyGlobalEquations_);
00257   
00258   for (i=0; i<NumMyCols_; i++)
00259       EPETRA_CHK_ERR(TempTransA1.InsertGlobalValues(TransMyGlobalEquations_[i], 
00260                      TransNumNz_[i], TransValues_[i], TransIndices_[i]));
00261  
00262   // Note: The following call to FillComplete is currently necessary because
00263   //     some global constants that are needed by the Export () are computed in this routine
00264   EPETRA_CHK_ERR(TempTransA1.FillComplete(false));
00265 
00266   // Now that transpose matrix with shared rows is entered, update values of target transpose matrix
00267   TransposeMatrix_->PutScalar(0.0);  // Zero out all values of the matrix
00268 
00269   EPETRA_CHK_ERR(TransposeMatrix_->Export(TempTransA1, *TransposeExporter_, Add));
00270 
00271   return(0);
00272 }
00273 
00274 //=========================================================================
00275 bool EpetraExt::RowMatrix_Transpose::rvs()
00276 {
00277   EPETRA_CHK_ERR(-1); //Not Implemented Yet
00278   return false;
00279 }
00280 
00281 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines