EpetraExt Package Browser (Single Doxygen Collection) 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 
00049 #include <vector>
00050 
00051 namespace EpetraExt {
00052 
00053 RowMatrix_Transpose::
00054 ~RowMatrix_Transpose()
00055 {
00056   if( TransposeMatrix_ ) delete TransposeMatrix_;
00057   if( TransposeExporter_ ) delete TransposeExporter_;
00058 
00059   if( !OrigMatrixIsCrsMatrix_ )
00060   {
00061     delete [] Indices_;
00062     delete [] Values_;
00063   }
00064 
00065   for( int i = 0; i < NumMyCols_; ++i )
00066     if( TransNumNz_[i] )
00067     {
00068       delete [] TransIndices_[i];
00069       delete [] TransValues_[i];
00070     }
00071 
00072   delete [] TransNumNz_;
00073   delete [] TransIndices_;
00074   delete [] TransValues_;
00075   delete [] TransMyGlobalEquations_;
00076 }
00077 
00078 //=========================================================================
00079 RowMatrix_Transpose::NewTypeRef
00080 RowMatrix_Transpose::
00081 operator()( OriginalTypeRef orig )
00082 {
00083   origObj_ = &orig;
00084 
00085   int i, j, err;
00086 
00087   if( !TransposeRowMap_ )
00088   {
00089     if( IgnoreNonLocalCols_ )
00090       TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorRangeMap()); // Should be replaced with refcount =
00091     else
00092       TransposeRowMap_ = (Epetra_Map *) &(orig.OperatorDomainMap()); // Should be replaced with refcount =
00093   }
00094 
00095   // This routine will work for any RowMatrix object, but will attempt cast the matrix to a CrsMatrix if
00096   // possible (because we can then use a View of the matrix and graph, which is much cheaper).
00097 
00098   // First get the local indices to count how many nonzeros will be in the 
00099   // transpose graph on each processor
00100   Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix*>(&orig);
00101 
00102   OrigMatrixIsCrsMatrix_ = (OrigCrsMatrix!=0); // If this pointer is non-zero, the cast to CrsMatrix worked
00103 
00104   NumMyRows_ = orig.NumMyRows();
00105   NumMyCols_ = orig.NumMyCols();
00106   TransNumNz_ = new int[NumMyCols_];
00107   TransIndices_ = new int*[NumMyCols_];
00108   TransValues_ = new double*[NumMyCols_];
00109   TransMyGlobalEquations_ = new int[NumMyCols_];
00110 
00111   int NumIndices;
00112 
00113   if (OrigMatrixIsCrsMatrix_)
00114   {
00115     const Epetra_CrsGraph & OrigGraph = OrigCrsMatrix->Graph(); // Get matrix graph
00116 
00117     for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0;
00118     for (i=0; i<NumMyRows_; i++)
00119     {
00120       err = OrigGraph.ExtractMyRowView(i, NumIndices, Indices_); // Get view of ith row
00121       if (err != 0) throw OrigGraph.ReportError("ExtractMyRowView failed",err);
00122       for (j=0; j<NumIndices; j++) ++TransNumNz_[Indices_[j]];
00123     }
00124   }
00125   else // Original is not a CrsMatrix
00126   {
00127     MaxNumEntries_ = 0;
00128     int NumEntries;
00129     for (i=0; i<NumMyRows_; i++)
00130     {
00131       orig.NumMyRowEntries(i, NumEntries);
00132       MaxNumEntries_ = EPETRA_MAX(MaxNumEntries_, NumEntries);
00133     }
00134     Indices_ = new int[MaxNumEntries_];
00135     Values_ = new double[MaxNumEntries_];
00136 
00137     for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0;
00138     for (i=0; i<NumMyRows_; i++)
00139     {
00140       err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_); 
00141       if (err != 0) {
00142         std::cerr << "ExtractMyRowCopy failed."<<std::endl;
00143         throw err;
00144       }
00145       for (j=0; j<NumIndices; j++) ++TransNumNz_[Indices_[j]];
00146     }
00147   }
00148 
00149   // Most of remaining code is common to both cases
00150   for(i=0; i<NumMyCols_; i++)
00151   {
00152     NumIndices = TransNumNz_[i];
00153     if (NumIndices>0)
00154     {
00155       TransIndices_[i] = new int[NumIndices];
00156       TransValues_[i] = new double[NumIndices];
00157     }
00158   }
00159 
00160   // Now copy values and global indices into newly create transpose storage
00161 
00162   for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0; // Reset transpose NumNz counter
00163   for (i=0; i<NumMyRows_; i++)
00164   {
00165     if (OrigMatrixIsCrsMatrix_)
00166       err = OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_);
00167     else
00168       err = orig.ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_);
00169     if (err != 0) {
00170       std::cerr << "ExtractMyRowCopy failed."<<std::endl;
00171       throw err;
00172     }
00173 
00174     int ii = orig.RowMatrixRowMap().GID(i);
00175     for (j=0; j<NumIndices; j++)
00176     {
00177       int TransRow = Indices_[j];
00178       int loc = TransNumNz_[TransRow];
00179       TransIndices_[TransRow][loc] = ii;
00180       TransValues_[TransRow][loc] = Values_[j];
00181       ++TransNumNz_[TransRow]; // increment counter into current transpose row
00182     }
00183   }
00184 
00185   //  Build Transpose matrix with some rows being shared across processors.
00186   //  We will use a view here since the matrix will not be used for anything else
00187 
00188   const Epetra_Map & TransMap = orig.RowMatrixColMap();
00189 
00190   Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz_);
00191   TransMap.MyGlobalElements(TransMyGlobalEquations_);
00192   
00193   for (i=0; i<NumMyCols_; i++) {
00194     err = TempTransA1.InsertGlobalValues(TransMyGlobalEquations_[i], 
00195                      TransNumNz_[i], TransValues_[i], TransIndices_[i]);
00196     if (err < 0) throw TempTransA1.ReportError("InsertGlobalValues failed.",err);
00197   }
00198  
00199   // Note: The following call to FillComplete is currently necessary because
00200   //      some global constants that are needed by the Export () are computed in this routine
00201   err = TempTransA1.FillComplete(orig.OperatorRangeMap(),*TransposeRowMap_, false);
00202   if (err != 0) {
00203     throw TempTransA1.ReportError("FillComplete failed.",err);
00204   }
00205 
00206   // Now that transpose matrix with shared rows is entered, create a new matrix that will
00207   // get the transpose with uniquely owned rows (using the same row distribution as A).
00208   if( IgnoreNonLocalCols_ )
00209     TransposeMatrix_ = new Epetra_CrsMatrix(Copy, *TransposeRowMap_, *TransposeRowMap_, 0);
00210   else
00211     TransposeMatrix_ = new Epetra_CrsMatrix(Copy, *TransposeRowMap_,0);
00212 
00213   // Create an Export object that will move TempTransA around
00214   TransposeExporter_ = new Epetra_Export(TransMap, *TransposeRowMap_);
00215 
00216   err = TransposeMatrix_->Export(TempTransA1, *TransposeExporter_, Add);
00217   if (err != 0) throw TransposeMatrix_->ReportError("Export failed.",err);
00218   
00219   err = TransposeMatrix_->FillComplete(orig.OperatorRangeMap(),*TransposeRowMap_);
00220   if (err != 0) throw TransposeMatrix_->ReportError("FillComplete failed.",err);
00221 
00222   if (MakeDataContiguous_) {
00223     err = TransposeMatrix_->MakeDataContiguous();
00224     if (err != 0) throw TransposeMatrix_->ReportError("MakeDataContiguous failed.",err);
00225   }
00226 
00227   newObj_ = TransposeMatrix_;
00228 
00229   return *newObj_;
00230 }
00231 
00232 //=========================================================================
00233 bool EpetraExt::RowMatrix_Transpose::fwd()
00234 {
00235   int i, j, NumIndices, err;
00236 
00237   Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix*>(origObj_);
00238 
00239   // Now copy values and global indices into newly create transpose storage
00240 
00241   for (i=0;i<NumMyCols_; i++) TransNumNz_[i] = 0; // Reset transpose NumNz counter
00242   for (i=0; i<NumMyRows_; i++)
00243   {
00244     if(OrigMatrixIsCrsMatrix_)
00245       err = OrigCrsMatrix->ExtractMyRowView(i, NumIndices, Values_, Indices_);
00246     else
00247       err = origObj_->ExtractMyRowCopy(i, MaxNumEntries_, NumIndices, Values_, Indices_);
00248     if (err != 0) {
00249       std::cerr << "ExtractMyRowCopy/View failed."<<std::endl;
00250       throw err;
00251     }
00252 
00253     int ii = origObj_->RowMatrixRowMap().GID(i);
00254     for (j=0; j<NumIndices; j++)
00255     {
00256       int TransRow = Indices_[j];
00257       int loc = TransNumNz_[TransRow];
00258       TransIndices_[TransRow][loc] = ii;
00259       TransValues_[TransRow][loc] = Values_[j];
00260       ++TransNumNz_[TransRow]; // increment counter into current transpose row
00261     }
00262   }
00263 
00264   //  Build Transpose matrix with some rows being shared across processors.
00265   //  We will use a view here since the matrix will not be used for anything else
00266   const Epetra_Map & TransMap = origObj_->RowMatrixColMap();
00267 
00268   Epetra_CrsMatrix TempTransA1(View, TransMap, TransNumNz_);
00269   TransMap.MyGlobalElements(TransMyGlobalEquations_);
00270   
00271   for (i=0; i<NumMyCols_; i++)
00272       EPETRA_CHK_ERR(TempTransA1.InsertGlobalValues(TransMyGlobalEquations_[i], 
00273                      TransNumNz_[i], TransValues_[i], TransIndices_[i]));
00274  
00275   // Note: The following call to FillComplete is currently necessary because
00276   //     some global constants that are needed by the Export () are computed in this routine
00277   EPETRA_CHK_ERR(TempTransA1.FillComplete(false));
00278 
00279   // Now that transpose matrix with shared rows is entered, update values of target transpose matrix
00280   TransposeMatrix_->PutScalar(0.0);  // Zero out all values of the matrix
00281 
00282   EPETRA_CHK_ERR(TransposeMatrix_->Export(TempTransA1, *TransposeExporter_, Add));
00283 
00284   return(0);
00285 }
00286 
00287 //=========================================================================
00288 bool EpetraExt::RowMatrix_Transpose::rvs()
00289 {
00290   EPETRA_CHK_ERR(-1); //Not Implemented Yet
00291   return false;
00292 }
00293 
00294 } // namespace EpetraExt
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines