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::CreateTransposeLocal(OriginalTypeRef orig) 
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: CreateTransposeLocal 1");
00090   mtime->start();
00091 #endif
00092 
00093   int i,j,err;
00094   const Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<const Epetra_CrsMatrix*>(&orig);
00095   if(OrigCrsMatrix) OrigMatrixIsCrsMatrix_ = true;
00096   else OrigMatrixIsCrsMatrix_ = false;
00097 
00098   const Epetra_Map & TransMap      = orig.RowMatrixColMap();
00099   int TransNnz                     = orig.NumMyNonzeros();
00100   int NumIndices;
00101 
00102   Epetra_CrsMatrix *TempTransA1 = new Epetra_CrsMatrix(Copy, TransMap,orig.RowMatrixRowMap(),0);
00103   Epetra_IntSerialDenseVector & TransRowptr = TempTransA1->ExpertExtractIndexOffset();
00104   Epetra_IntSerialDenseVector & TransColind = TempTransA1->ExpertExtractIndices();
00105   double *&                     TransVals   = TempTransA1->ExpertExtractValues();  
00106   NumMyRows_ = orig.NumMyRows();
00107   NumMyCols_ = orig.NumMyCols();
00108 
00109   TransRowptr.Resize(NumMyCols_+1);
00110   TransColind.Resize(TransNnz);
00111   resize_doubles(0,TransNnz,TransVals);
00112   std::vector<int> CurrentStart(NumMyCols_,0);
00113 
00114   // Count up nnz using the Rowptr to count the number of non-nonzeros.
00115   if (OrigMatrixIsCrsMatrix_)
00116   {
00117     const Epetra_CrsGraph & OrigGraph = OrigCrsMatrix->Graph(); // Get matrix graph
00118 
00119     for (i=0; i<NumMyRows_; i++)
00120     {
00121       err = OrigGraph.ExtractMyRowView(i, NumIndices, Indices_); // Get view of ith row
00122       if (err != 0) throw OrigGraph.ReportError("ExtractMyRowView failed",err);
00123       for (j=0; j<NumIndices; j++) ++CurrentStart[Indices_[j]];
00124     }
00125   }
00126   else // Original is not a CrsMatrix
00127   {
00128     MaxNumEntries_ = 0;
00129     MaxNumEntries_ = orig.MaxNumEntries();
00130     delete [] Indices_; delete [] Values_;
00131     Indices_ = new int[MaxNumEntries_];
00132     Values_  = new double[MaxNumEntries_];
00133 
00134     for (i=0; i<NumMyRows_; i++)
00135     {
00136       err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_); 
00137       if (err != 0) {
00138         std::cerr << "ExtractMyRowCopy failed."<<std::endl;
00139         throw err;
00140       }
00141       for (j=0; j<NumIndices; j++) ++CurrentStart[Indices_[j]];
00142     }
00143   }
00144 
00145   // Scansum the rowptr; reset currentstart
00146   TransRowptr[0] = 0;
00147   for (i=1;i<NumMyCols_+1; i++)  TransRowptr[i]   = CurrentStart[i-1] + TransRowptr[i-1];
00148   for (i=0;i<NumMyCols_;   i++)  CurrentStart[i]  = TransRowptr[i];
00149 
00150   // Now copy values and global indices into newly create transpose storage
00151   for (i=0; i<NumMyRows_; i++)
00152   {
00153     if (OrigMatrixIsCrsMatrix_)
00154       err = OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_);
00155     else
00156       err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_);
00157     if (err != 0) {
00158       std::cerr << "ExtractMyRowCopy failed."<<std::endl;
00159       throw err;
00160     }
00161 
00162     for (j=0; j<NumIndices; j++)
00163     {
00164       int idx = CurrentStart[Indices_[j]];
00165       TransColind[idx] = i;
00166       TransVals[idx]    = Values_[j];
00167       ++CurrentStart[Indices_[j]];// increment the counter into the local row
00168     }
00169   }
00170 
00171 #ifdef ENABLE_TRANSPOSE_TIMINGS
00172   mtime->stop();
00173   mtime=MM.getNewTimer("Transpose: CreateTransposeLocal 2");
00174   mtime->start();
00175 #endif
00176 
00177   // Prebuild the importers and exporters the no-communication way, flipping the importers
00178   // and exporters around.
00179   Epetra_Import * myimport = 0;
00180   Epetra_Export * myexport = 0;
00181   if(OrigMatrixIsCrsMatrix_ && OrigCrsMatrix->Importer())
00182     myexport = new Epetra_Export(*OrigCrsMatrix->Importer());
00183   if(OrigMatrixIsCrsMatrix_ && OrigCrsMatrix->Exporter())
00184     myimport = new Epetra_Import(*OrigCrsMatrix->Exporter());
00185 
00186 #ifdef ENABLE_TRANSPOSE_TIMINGS
00187   mtime->stop();
00188   mtime=MM.getNewTimer("Transpose: CreateTransposeLocal 3");
00189   mtime->start();
00190 #endif
00191 
00192   // Call ExpertStaticFillComplete
00193   err = TempTransA1->ExpertStaticFillComplete(orig.OperatorRangeMap(),orig.OperatorDomainMap(),myimport,myexport);
00194   if (err != 0) {
00195     throw TempTransA1->ReportError("ExpertStaticFillComplete failed.",err);
00196   }
00197 
00198 #ifdef ENABLE_TRANSPOSE_TIMINGS
00199   mtime->stop();
00200 #endif
00201  
00202   return TempTransA1;
00203 }
00204 
00205 
00206 //=========================================================================
00207 RowMatrix_Transpose::NewTypeRef
00208 RowMatrix_Transpose::
00209 operator()( OriginalTypeRef orig )
00210 {
00211   origObj_ = &orig;
00212 
00213   if( !TransposeRowMap_ )
00214   {
00215     if( IgnoreNonLocalCols_ )
00216       TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorRangeMap()); // Should be replaced with refcount =
00217     else
00218       TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorDomainMap()); // Should be replaced with refcount =
00219   }
00220 
00221   NumMyRows_ = orig.NumMyRows();
00222   NumMyCols_ = orig.NumMyCols();
00223 
00224   // This routine will work for any RowMatrix object, but will attempt cast the matrix to a CrsMatrix if
00225   // possible (because we can then use a View of the matrix and graph, which is much cheaper).
00226 
00227   // First get the local indices to count how many nonzeros will be in the 
00228   // transpose graph on each processor
00229   Epetra_CrsMatrix * TempTransA1 = CreateTransposeLocal(orig);  
00230 
00231   if(!TempTransA1->Exporter()) {
00232     // Short Circuit: There is no need to make another matrix since TransposeRowMap_== TransMap
00233     newObj_ = TransposeMatrix_ = TempTransA1;
00234     return *newObj_;    
00235   }
00236 
00237 #ifdef ENABLE_TRANSPOSE_TIMINGS
00238   Teuchos::Time myTime("global");
00239   Teuchos::TimeMonitor MM(myTime);
00240   Teuchos::RCP<Teuchos::Time> mtime;
00241   mtime=MM.getNewTimer("Transpose: Final FusedExport");
00242   mtime->start();
00243 #endif
00244 
00245   // Now that transpose matrix with shared rows is entered, create a new matrix that will
00246   // get the transpose with uniquely owned rows (using the same row distribution as A).
00247   TransposeMatrix_ = new Epetra_CrsMatrix(*TempTransA1,*TempTransA1->Exporter(),0,TransposeRowMap_);
00248 
00249 #ifdef ENABLE_TRANSPOSE_TIMINGS
00250   mtime->stop();
00251 #endif
00252 
00253   newObj_ = TransposeMatrix_;
00254   delete TempTransA1;
00255   return *newObj_;
00256 }
00257 
00258 //=========================================================================
00259 bool EpetraExt::RowMatrix_Transpose::fwd()
00260 {
00261   Epetra_CrsMatrix * TempTransA1 = CreateTransposeLocal(*origObj_);  
00262   const Epetra_Export * TransposeExporter=0;
00263   bool DeleteExporter = false;
00264   
00265   if(TempTransA1->Exporter()) TransposeExporter = TempTransA1->Exporter();
00266   else {
00267     DeleteExporter=true;
00268     TransposeExporter = new Epetra_Export(TransposeMatrix_->DomainMap(),TransposeMatrix_->RowMap());
00269   }
00270 
00271   TransposeMatrix_->PutScalar(0.0);  // Zero out all values of the matrix
00272 
00273   EPETRA_CHK_ERR(TransposeMatrix_->Export(*TempTransA1, *TransposeExporter, Add));
00274 
00275   if(DeleteExporter) delete TransposeExporter;
00276   delete TempTransA1;
00277   return 0;
00278 }
00279 
00280 //=========================================================================
00281 bool EpetraExt::RowMatrix_Transpose::rvs()
00282 {
00283   EPETRA_CHK_ERR(-1); //Not Implemented Yet
00284   return false;
00285 }
00286 
00287 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines