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 (2011) Sandia Corporation
00006 //
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ***********************************************************************
00040 //@HEADER
00041 
00042 #include <EpetraExt_Transpose_RowMatrix.h>
00043 
00044 #include <Epetra_Export.h>
00045 #include <Epetra_CrsGraph.h>
00046 #include <Epetra_CrsMatrix.h>
00047 #include <Epetra_Map.h>
00048 #include <Epetra_Import.h>
00049 #include <Epetra_Export.h>
00050 
00051 #include <Teuchos_TimeMonitor.hpp>
00052 #include <vector>
00053 
00054 //#define ENABLE_TRANSPOSE_TIMINGS
00055 
00056 namespace EpetraExt {
00057 
00058 // Provide a "resize" operation for double*'s. 
00059 inline void resize_doubles(int nold,int nnew,double*& d){
00060   if(nnew > nold){
00061     double *tmp = new double[nnew];
00062     for(int i=0; i<nold; i++)
00063       tmp[i]=d[i];
00064     delete [] d;
00065     d=tmp;
00066   }
00067 }
00068 
00069 
00070 RowMatrix_Transpose::
00071 ~RowMatrix_Transpose()
00072 {
00073   if( TransposeMatrix_ ) delete TransposeMatrix_;
00074 
00075   if( !OrigMatrixIsCrsMatrix_ )
00076   {
00077     delete [] Indices_;
00078     delete [] Values_;
00079   }
00080 }
00081 
00082 //=========================================================================
00083   Epetra_CrsMatrix* EpetraExt::RowMatrix_Transpose::BuildTempTrans() 
00084 {
00085 #ifdef ENABLE_TRANSPOSE_TIMINGS
00086   Teuchos::Time myTime("global");
00087   Teuchos::TimeMonitor MM(myTime);
00088   Teuchos::RCP<Teuchos::Time> mtime;
00089   mtime=MM.getNewTimer("Transpose: BuildTempTrans 1");
00090   mtime->start();
00091 #endif
00092 
00093   int i,j,err;
00094   const Epetra_RowMatrix & orig    = *origObj_;
00095   const Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<const Epetra_CrsMatrix*>(&orig);
00096   const Epetra_Map & TransMap      = orig.RowMatrixColMap();
00097   int TransNnz                     = orig.NumMyNonzeros();
00098   int NumIndices;
00099 
00100   Epetra_CrsMatrix *TempTransA1 = new Epetra_CrsMatrix(Copy, TransMap,orig.RowMatrixRowMap(),0);
00101   Epetra_IntSerialDenseVector & TransRowptr = TempTransA1->ExpertExtractIndexOffset();
00102   Epetra_IntSerialDenseVector & TransColind = TempTransA1->ExpertExtractIndices();
00103   double *&                     TransVals   = TempTransA1->ExpertExtractValues();  
00104 
00105   TransRowptr.Resize(NumMyCols_+1);
00106   TransColind.Resize(TransNnz);
00107   resize_doubles(0,TransNnz,TransVals);
00108   std::vector<int> CurrentStart(NumMyCols_,0);
00109 
00110   // Count up nnz using the Rowptr to count the number of non-nonzeros.
00111   if (OrigMatrixIsCrsMatrix_)
00112   {
00113     const Epetra_CrsGraph & OrigGraph = OrigCrsMatrix->Graph(); // Get matrix graph
00114 
00115     for (i=0; i<NumMyRows_; i++)
00116     {
00117       err = OrigGraph.ExtractMyRowView(i, NumIndices, Indices_); // Get view of ith row
00118       if (err != 0) throw OrigGraph.ReportError("ExtractMyRowView failed",err);
00119       for (j=0; j<NumIndices; j++) ++CurrentStart[Indices_[j]];
00120     }
00121   }
00122   else // Original is not a CrsMatrix
00123   {
00124     MaxNumEntries_ = 0;
00125     MaxNumEntries_ = orig.MaxNumEntries();
00126     delete [] Indices_; delete [] Values_;
00127     Indices_ = new int[MaxNumEntries_];
00128     Values_  = new double[MaxNumEntries_];
00129 
00130     for (i=0; i<NumMyRows_; i++)
00131     {
00132       err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_); 
00133       if (err != 0) {
00134         std::cerr << "ExtractMyRowCopy failed."<<std::endl;
00135         throw err;
00136       }
00137       for (j=0; j<NumIndices; j++) ++CurrentStart[Indices_[j]];
00138     }
00139   }
00140 
00141   // Scansum the rowptr; reset currentstart
00142   TransRowptr[0] = 0;
00143   for (i=1;i<NumMyCols_+1; i++)  TransRowptr[i]   = CurrentStart[i-1] + TransRowptr[i-1];
00144   for (i=0;i<NumMyCols_;   i++)  CurrentStart[i]  = TransRowptr[i];
00145 
00146   // Now copy values and global indices into newly create transpose storage
00147   for (i=0; i<NumMyRows_; i++)
00148   {
00149     if (OrigMatrixIsCrsMatrix_)
00150       err = OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_);
00151     else
00152       err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_);
00153     if (err != 0) {
00154       std::cerr << "ExtractMyRowCopy failed."<<std::endl;
00155       throw err;
00156     }
00157 
00158     for (j=0; j<NumIndices; j++)
00159     {
00160       int idx = CurrentStart[Indices_[j]];
00161       TransColind[idx] = i;
00162       TransVals[idx]    = Values_[j];
00163       ++CurrentStart[Indices_[j]];// increment the counter into the local row
00164     }
00165   }
00166 
00167 #ifdef ENABLE_TRANSPOSE_TIMINGS
00168   mtime->stop();
00169   mtime=MM.getNewTimer("Transpose: BuildTempTrans 2");
00170   mtime->start();
00171 #endif
00172 
00173   // Prebuild the importers and exporters the no-communication way, flipping the importers
00174   // and exporters around.
00175   Epetra_Import * myimport = 0;
00176   Epetra_Export * myexport = 0;
00177   if(OrigMatrixIsCrsMatrix_ && OrigCrsMatrix->Importer())
00178     myexport = new Epetra_Export(*OrigCrsMatrix->Importer());
00179   if(OrigMatrixIsCrsMatrix_ && OrigCrsMatrix->Exporter())
00180     myimport = new Epetra_Import(*OrigCrsMatrix->Exporter());
00181 
00182 #ifdef ENABLE_TRANSPOSE_TIMINGS
00183   mtime->stop();
00184   mtime=MM.getNewTimer("Transpose: BuildTempTrans 3");
00185   mtime->start();
00186 #endif
00187 
00188   // Call ExpertStaticFillComplete
00189   err = TempTransA1->ExpertStaticFillComplete(orig.OperatorRangeMap(),*TransposeRowMap_,myimport,myexport);
00190   if (err != 0) {
00191     throw TempTransA1->ReportError("ExpertStaticFillComplete failed.",err);
00192   }
00193 
00194 #ifdef ENABLE_TRANSPOSE_TIMINGS
00195   mtime->stop();
00196 #endif
00197  
00198   return TempTransA1;
00199 }
00200 
00201 
00202 //=========================================================================
00203 RowMatrix_Transpose::NewTypeRef
00204 RowMatrix_Transpose::
00205 operator()( OriginalTypeRef orig )
00206 {
00207   origObj_ = &orig;
00208 
00209   if( !TransposeRowMap_ )
00210   {
00211     if( IgnoreNonLocalCols_ )
00212       TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorRangeMap()); // Should be replaced with refcount =
00213     else
00214       TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorDomainMap()); // Should be replaced with refcount =
00215   }
00216 
00217   NumMyRows_ = orig.NumMyRows();
00218   NumMyCols_ = orig.NumMyCols();
00219 
00220   // This routine will work for any RowMatrix object, but will attempt cast the matrix to a CrsMatrix if
00221   // possible (because we can then use a View of the matrix and graph, which is much cheaper).
00222 
00223   // First get the local indices to count how many nonzeros will be in the 
00224   // transpose graph on each processor
00225   Epetra_CrsMatrix * TempTransA1 = BuildTempTrans();  
00226 
00227   if(!TempTransA1->Exporter()) {
00228     // Short Circuit: There is no need to make another matrix since TransposeRowMap_== TransMap
00229     newObj_ = TransposeMatrix_ = TempTransA1;
00230     return *newObj_;    
00231   }
00232 
00233 #ifdef ENABLE_TRANSPOSE_TIMINGS
00234   Teuchos::Time myTime("global");
00235   Teuchos::TimeMonitor MM(myTime);
00236   Teuchos::RCP<Teuchos::Time> mtime;
00237   mtime=MM.getNewTimer("Transpose: Final FusedExport");
00238   mtime->start();
00239 #endif
00240 
00241   // Now that transpose matrix with shared rows is entered, create a new matrix that will
00242   // get the transpose with uniquely owned rows (using the same row distribution as A).
00243   TransposeMatrix_ = new Epetra_CrsMatrix(*TempTransA1,*TempTransA1->Exporter(),0,TransposeRowMap_);
00244 
00245 #ifdef ENABLE_TRANSPOSE_TIMINGS
00246   mtime->stop();
00247 #endif
00248 
00249   newObj_ = TransposeMatrix_;
00250   delete TempTransA1;
00251   return *newObj_;
00252 }
00253 
00254 //=========================================================================
00255 bool EpetraExt::RowMatrix_Transpose::fwd()
00256 {
00257   Epetra_CrsMatrix * TempTransA1 = BuildTempTrans();  
00258   const Epetra_Export * TransposeExporter=0;
00259   bool DeleteExporter = false;
00260   
00261   if(TempTransA1->Exporter()) TransposeExporter = TempTransA1->Exporter();
00262   else {
00263     DeleteExporter=true;
00264     TransposeExporter = new Epetra_Export(TransposeMatrix_->DomainMap(),TransposeMatrix_->RowMap());
00265   }
00266 
00267   TransposeMatrix_->PutScalar(0.0);  // Zero out all values of the matrix
00268 
00269   EPETRA_CHK_ERR(TransposeMatrix_->Export(*TempTransA1, *TransposeExporter, Add));
00270 
00271   if(DeleteExporter) delete TransposeExporter;
00272   delete TempTransA1;
00273   return 0;
00274 }
00275 
00276 //=========================================================================
00277 bool EpetraExt::RowMatrix_Transpose::rvs()
00278 {
00279   EPETRA_CHK_ERR(-1); //Not Implemented Yet
00280   return false;
00281 }
00282 
00283 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines