Epetra_CrsMatrix.h

Go to the documentation of this file.
00001 
00002 //@HEADER
00003 /*
00004 ************************************************************************
00005 
00006               Epetra: Linear Algebra Services Package 
00007                 Copyright (2001) Sandia Corporation
00008 
00009 Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00010 license for use of this work by or on behalf of the U.S. Government.
00011 
00012 This library is free software; you can redistribute it and/or modify
00013 it under the terms of the GNU Lesser General Public License as
00014 published by the Free Software Foundation; either version 2.1 of the
00015 License, or (at your option) any later version.
00016  
00017 This library is distributed in the hope that it will be useful, but
00018 WITHOUT ANY WARRANTY; without even the implied warranty of
00019 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 Lesser General Public License for more details.
00021  
00022 You should have received a copy of the GNU Lesser General Public
00023 License along with this library; if not, write to the Free Software
00024 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 USA
00026 Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00027 
00028 ************************************************************************
00029 */
00030 //@HEADER
00031 
00032 #ifndef EPETRA_CRSMATRIX_H
00033 #define EPETRA_CRSMATRIX_H
00034 
00035 #include "Epetra_DistObject.h" 
00036 #include "Epetra_CompObject.h" 
00037 #include "Epetra_BLAS.h"
00038 #include "Epetra_RowMatrix.h"
00039 #include "Epetra_Operator.h"
00040 #include "Epetra_CrsGraph.h"
00041 class Epetra_Map;
00042 class Epetra_Import;
00043 class Epetra_Export;
00044 class Epetra_Vector;
00045 class Epetra_MultiVector;
00046 
00047 // Define this to see a complete dump a an Epetra_CrsMatrix::Multiply(...) call
00048 //#define EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
00049 
00050 #ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
00051 extern bool Epetra_CrsMatrixTraceDumpMultiply;
00052 #endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
00053 
00055 
00152 class Epetra_CrsMatrix: public Epetra_DistObject, public Epetra_CompObject, public Epetra_BLAS, public virtual Epetra_RowMatrix {
00153  public:
00154 
00156 
00157 
00158 
00169   Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& RowMap, const int* NumEntriesPerRow, bool StaticProfile = false);
00170   
00172 
00184   Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& RowMap, int NumEntriesPerRow, bool StaticProfile = false);
00185 
00187 
00200   Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& RowMap, const Epetra_Map& ColMap, const int* NumEntriesPerRow, bool StaticProfile = false);
00201   
00203 
00217   Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& RowMap, const Epetra_Map& ColMap, int NumEntriesPerRow, bool StaticProfile = false);
00218 
00220 
00226   Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_CrsGraph& Graph);
00227   
00229   Epetra_CrsMatrix(const Epetra_CrsMatrix& Matrix);
00230   
00232   virtual ~Epetra_CrsMatrix();
00234   
00236 
00237 
00239   Epetra_CrsMatrix& operator=(const Epetra_CrsMatrix& src);
00240 
00242 
00249   int PutScalar(double ScalarConstant);
00250   
00252 
00259   int Scale(double ScalarConstant);
00260   
00262 
00288   virtual int InsertGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices);
00289 
00291 
00303   virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices);
00304   
00306 
00318   virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices);
00319 
00321 
00334   int InsertMyValues(int MyRow, int NumEntries, double* Values, int* Indices);
00335 
00337 
00349   int ReplaceMyValues(int MyRow, int NumEntries, double* Values, int* Indices);
00350 
00352 
00365   int SumIntoMyValues(int MyRow, int NumEntries, double* Values, int* Indices);
00366 
00368 
00381   int ReplaceDiagonalValues(const Epetra_Vector& Diagonal);
00382   
00384 
00386 
00387   
00388 
00390   /* This version of FillComplete assumes that the domain and range
00391      distributions are identical to the matrix row distributions.
00392     \param OptimizeDataStorage - (In) If true, storage will be packed for optimal performance.  Depending
00393            on how the matrix was constructed, optimizing the storage may have no impact on performance
00394      or one-time memory use, or may have a large impact.  If the user was careful in allocating memory
00395      for the matrix by setting StaticProfile to true in the matrix constructor, then no extra storage
00396      will be allocated in attempting to optimize storage.  If the user did not set StaticProfile to true,
00397      then optimizing the storage will temporarily use additional memory, will have a noticeable impact
00398      on performance and ultimately reduce the storage associated with the matrix.
00399 
00400      By default storage will be optimized.  If you cannot tolerate the increased temporary memory use,
00401      should set this value to false.
00402 
00403      \return error code, 0 if successful. Returns a positive warning code of 3
00404         if the matrix is rectangular (meaning that the other overloading of
00405         FillComplete should have been called, with differen domain-map and
00406         range-map specified).
00407   */
00408   int FillComplete(bool OptimizeDataStorage = true);
00409 
00411   /* This version of FillComplete requires the explicit specification of the domain
00412      and range distribution maps.  These maps are used for importing and exporting vector
00413      and multi-vector elements that are needed for distributed matrix computations.  For
00414      example, to compute y = Ax in parallel, we would specify the DomainMap as the distribution
00415      of the vector x and the RangeMap as the distribution of the vector y.
00416     \param DomainMap - (In) Map that describes the distribution of vector and multi-vectors in the
00417     matrix domain.
00418     \param RangeMap - (In) Map that describes the distribution of vector and multi-vectors in the
00419     matrix range.
00420 
00421     \param OptimizeDataStorage - (In) If true, storage will be packed for optimal performance.  Depending
00422            on how the matrix was constructed, optimizing the storage may have no impact on performance
00423      or one-time memory use, or may have a large impact.  If the user was careful in allocating memory
00424      for the matrix by setting StaticProfile to true in the matrix constructor, then no extra storage
00425      will be allocated in attempting to optimize storage.  If the user did not set StaticProfile to true,
00426      then optimizing the storage will temporarily use additional memory, will have a noticeable impact
00427      on performance and ultimately reduce the storage associated with the matrix.
00428 
00429      By default storage will be optimized.  If you cannot tolerate the increased temporary memory use,
00430      should set this value to false.
00431 
00432     \return error code, 0 if successful. positive warning code of 2 if it is detected that the
00433     matrix-graph got out of sync since this matrix was constructed (for instance if
00434     graph.FillComplete() was called by another matrix that shares the graph)
00435 
00436     \post IndicesAreLocal()==true
00437     */
00438   int FillComplete(const Epetra_Map& DomainMap, const Epetra_Map& RangeMap, bool OptimizeDataStorage = true);
00439     
00441 
00461   int OptimizeStorage();
00462 
00463   
00465   int MakeDataContiguous() {EPETRA_CHK_ERR(OptimizeStorage()); return(0);}
00467   
00469 
00470   
00472 
00482   int ExtractGlobalRowCopy(int GlobalRow, int Length, int& NumEntries, double* Values, int* Indices) const;
00483 
00485 
00496   int ExtractMyRowCopy(int MyRow, int Length, int& NumEntries, double* Values, int* Indices) const;
00497 
00499 
00507   int ExtractGlobalRowCopy(int GlobalRow, int Length, int& NumEntries, double* Values) const;
00508 
00510 
00518   int ExtractMyRowCopy(int MyRow, int Length, int& NumEntries, double* Values) const;
00519 
00521 
00528   int ExtractDiagonalCopy(Epetra_Vector& Diagonal) const;
00529   
00531 
00542   int ExtractGlobalRowView(int GlobalRow, int& NumEntries, double*& Values, int*& Indices) const;
00543   
00545 
00556   int ExtractMyRowView(int MyRow, int& NumEntries, double*& Values, int*& Indices) const;
00557   
00559 
00566   int ExtractGlobalRowView(int GlobalRow, int& NumEntries, double*& Values) const;
00567 
00569 
00576   int ExtractMyRowView(int MyRow, int& NumEntries, double*& Values) const;
00578   
00580 
00581   
00583 
00592   int Multiply(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const;
00593   int Multiply1(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const;
00594 
00595 
00597 
00606   int Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00607   int Multiply1(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00608 
00610 
00621   int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const;
00622   
00624 
00635   int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
00636   
00638 
00654   int InvRowSums(Epetra_Vector& x) const;
00655 
00657 
00669   int InvRowMaxs(Epetra_Vector& x) const;
00670   
00672 
00680   int LeftScale(const Epetra_Vector& x);
00681   
00683 
00700   int InvColSums(Epetra_Vector& x) const;
00701 
00703 
00715   int InvColMaxs(Epetra_Vector& x) const;
00716 
00718 
00726   int RightScale(const Epetra_Vector& x);
00728   
00730 
00731   
00732   
00734   bool Filled() const {return(Graph_.Filled());}
00735   
00737   bool StorageOptimized() const {return(StorageOptimized_);}
00738   
00740   bool IndicesAreGlobal() const {return(Graph_.IndicesAreGlobal());}
00741   
00743   bool IndicesAreLocal() const {return(Graph_.IndicesAreLocal());}
00744   
00746   bool IndicesAreContiguous() const {return(Graph_.IndicesAreContiguous());}
00747   
00749   bool LowerTriangular() const {return(Graph_.LowerTriangular());}
00750   
00752   bool UpperTriangular() const {return(Graph_.UpperTriangular());}
00753   
00755   bool NoDiagonal() const {return(Graph_.NoDiagonal());}
00756   
00758   
00760 
00761   
00763   /* Returns the quantity \f$ \| A \|_\infty\f$ such that
00764      \f[\| A \|_\infty = \max_{1\lei\lem} \sum_{j=1}^n |a_{ij}| \f]
00765      \warning The NormInf() method will not properly calculate the infinity norm for a matrix that has entries that are
00766      replicated on multiple processors. */ 
00767   double NormInf() const;
00768   
00770   /* Returns the quantity \f$ \| A \|_1\f$ such that
00771      \f[\| A \|_1= \max_{1\lej\len} \sum_{i=1}^m |a_{ij}| \f].
00772      \warning The NormOne() method will not properly calculate the one norm for a matrix that has entries that are
00773      replicated on multiple processors.
00774   */ 
00775   double NormOne() const;
00776 
00778   /* Returns the quantity \f[ \| A \|_{Frobenius} = \sqrt{\sum_{i=1}^m \sum_{j=1}^n\|a_{ij}\|^2}\f]
00779      \warning the NormFrobenius() method will not properly calculate the frobenius norm for a matrix that
00780      has entries which are replicated on multiple processors. In that case, the returned
00781      norm will be larger than the true norm.
00782    */
00783   double NormFrobenius() const;
00784 
00786   /*
00787     Note that if maps are defined such that some nonzeros appear on
00788     multiple processors, then those nonzeros will be counted multiple times.
00789     If the user wishes to assemble a matrix from overlapping submatrices,
00790     they can use Epetra_FECrsMatrix.
00791   */
00792   int NumGlobalNonzeros() const {return(Graph_.NumGlobalNonzeros());}
00793   
00795   int NumGlobalRows() const {return(Graph_.NumGlobalRows());}
00796   
00798   int NumGlobalCols() const {return(Graph_.NumGlobalCols());}
00799   
00801   int NumGlobalDiagonals() const {return(Graph_.NumGlobalDiagonals());}
00802   
00804   int NumMyNonzeros() const {return(Graph_.NumMyNonzeros());}
00805   
00807   int NumMyRows() const {return(Graph_.NumMyRows());}
00808   
00810 
00814   int NumMyCols() const {return(Graph_.NumMyCols());}
00815   
00817 
00820   int NumMyDiagonals() const {return(Graph_.NumMyDiagonals());}
00821   
00823   int NumGlobalEntries(int Row) const {return(Graph_.NumGlobalIndices(Row));}
00824   
00826   int NumAllocatedGlobalEntries(int Row) const{return(Graph_.NumAllocatedGlobalIndices(Row));}
00827   
00829 
00832   int MaxNumEntries() const {return(Graph_.MaxNumIndices());}
00833 
00835 
00838   int GlobalMaxNumEntries() const {return(Graph_.GlobalMaxNumIndices());}
00839   
00841   int NumMyEntries(int Row) const {return(Graph_.NumMyIndices(Row));}
00842   
00844   int NumAllocatedMyEntries(int Row) const {return(Graph_.NumAllocatedMyIndices(Row));}
00845   
00847   int IndexBase() const {return(Graph_.IndexBase());}
00848   
00849   
00851   bool StaticGraph() {return(StaticGraph_);}
00852 
00854   const Epetra_CrsGraph& Graph() const {return(Graph_);}
00855   
00857   const Epetra_Map& RowMap() const {return((Epetra_Map &)Graph_.RowMap());}
00858 
00860 
00866   int ReplaceRowMap(const Epetra_BlockMap& newmap);
00867 
00869 
00873   bool HaveColMap() const {return(Graph_.HaveColMap());}
00874 
00876 
00882   int ReplaceColMap(const Epetra_BlockMap& newmap);
00883 
00884 
00886 
00893   const Epetra_Map& ColMap() const {return((Epetra_Map &) Graph_.ColMap());}
00894   
00896 
00899   const Epetra_Map& DomainMap() const {return((Epetra_Map &)Graph_.DomainMap());}
00900   
00902 
00905   const Epetra_Map& RangeMap() const  {return((Epetra_Map &)Graph_.RangeMap());}
00906   
00908   const Epetra_Import* Importer() const {return(Graph_.Importer());}
00909   
00911   const Epetra_Export* Exporter() const {return(Graph_.Exporter());}
00912   
00914   const Epetra_Comm& Comm() const {return(Epetra_DistObject::Comm());}
00916   
00918 
00919 
00920   int LRID( int GRID_in) const {return(Graph_.LRID(GRID_in));}
00921   
00923   int GRID( int LRID_in) const {return(Graph_.GRID(LRID_in));}
00924   
00926 
00929   int LCID( int GCID_in) const {return(Graph_.LCID(GCID_in));}
00930   
00932 
00935   int GCID( int LCID_in) const {return(Graph_.GCID(LCID_in));}
00936   
00938   bool MyGRID(int GRID_in) const {return(Graph_.MyGRID(GRID_in));}
00939   
00941   bool MyLRID(int LRID_in) const {return(Graph_.MyLRID(LRID_in));}
00942   
00944 
00947   bool MyGCID(int GCID_in) const {return(Graph_.MyGCID(GCID_in));}
00948    
00950 
00953   bool MyLCID(int LCID_in) const {return(Graph_.MyLCID(LCID_in));}
00954 
00956   bool MyGlobalRow(int GID) const {return(Graph_.MyGlobalRow(GID));}
00958   
00959   
00961 
00962 
00964   virtual void Print(ostream& os) const;
00966 
00968 
00969 
00971   const char* Label() const {return(Epetra_Object::Label());}
00972   
00974 
00982   int SetUseTranspose(bool UseTranspose_in) {UseTranspose_ = UseTranspose_in; return(0);}
00983 
00985 
00993   int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
00994     return(Epetra_CrsMatrix::Multiply(Epetra_CrsMatrix::UseTranspose(), X, Y));}
00995   
00997 
01010   int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
01011     return(Solve(UpperTriangular(), Epetra_CrsMatrix::UseTranspose(), NoDiagonal(), X, Y));}
01012 
01014   bool HasNormInf() const {return(true);}
01015   
01017   bool UseTranspose() const {return(UseTranspose_);}
01018   
01020   const Epetra_Map& OperatorDomainMap() const
01021         {
01022           if (UseTranspose()) return(RangeMap());
01023           else return(DomainMap());
01024         }
01025   
01027   const Epetra_Map& OperatorRangeMap() const
01028         {
01029           if (UseTranspose()) return(DomainMap());
01030           else return(RangeMap());
01031         }
01032 
01034 
01035 
01036 
01038 
01047   int NumMyRowEntries(int MyRow, int& NumEntries) const;
01048 
01050   const Epetra_Map& RowMatrixRowMap() const {return(RowMap());}
01051   
01053   const Epetra_Map& RowMatrixColMap() const {return(ColMap());}
01054   
01056   const Epetra_Import* RowMatrixImporter() const {return(Importer());}
01057   
01059   
01061 
01062   
01064 
01069   inline double* operator[] (int Loc) { 
01070     if (StorageOptimized()){ int * ind = Graph().IndexOffset(); return(All_Values_+ind[Loc]);}
01071     else return Values_[Loc];}
01072   inline double* operator[] (int Loc) const {
01073     if (StorageOptimized()){ int * ind = Graph().IndexOffset(); return(All_Values_+ind[Loc]);}
01074     else return Values_[Loc];}
01076   
01078 
01079 
01080 
01093   int ExtractCrsDataPointers(int *& IndexOffset, int *& Indices, double *& Values_in) const {
01094     if (StorageOptimized()) { 
01095       IndexOffset = Graph().IndexOffset();
01096       Indices = Graph().All_Indices();
01097       Values_in  = All_Values(); 
01098       return (0);
01099     } 
01100     else { IndexOffset = 0; Indices = 0; Values_in  = 0; return (-1);} }
01101   
01103 
01113   int SortGhostsAssociatedWithEachProcessor(bool Flag)  {Graph_.SortGhostsAssociatedWithEachProcessor(Flag); return(0);}
01115   
01117 
01118   
01120   const Epetra_Map& ImportMap() const {return((Epetra_Map&) Graph_.ImportMap());}
01121 
01123   int TransformToLocal();
01124 
01126   int TransformToLocal(const Epetra_Map* DomainMap, const Epetra_Map* RangeMap);
01127 
01129 
01130 
01131  protected:
01132   bool Allocated() const {return(Allocated_);}
01133   int SetAllocated(bool Flag) {Allocated_ = Flag; return(0);}
01134   double** Values() const {
01135     if (StorageOptimized()) throw ReportError("This method: double** Values() cannot be called when StorageOptimized()==true", -1);
01136     else return(Values_);}
01137   double* All_Values() const {
01138     if (!StorageOptimized()) throw ReportError("This method: double* All_Values()cannot be called when StorageOptimized()==false", -1);
01139     else return(All_Values_);}
01140   double* Values(int LocalRow) const {
01141     if (StorageOptimized())
01142       if (Graph().StorageOptimized())
01143   return(All_Values_+Graph().IndexOffset()[LocalRow]);
01144       else throw ReportError("This method: double* Values()cannot be called when StorageOptimized()==true and Graph().StorageOptimized()==false", -1);
01145     else return(Values_[LocalRow]);}
01146   
01147   void InitializeDefaults();
01148   int Allocate();
01149 
01150   int InsertValues(int LocalRow, int NumEntries, double* Values, int* Indices);
01151 
01152   int InsertOffsetValues(int GlobalRow, int NumEntries, double *Values, int *Indices);
01153   int ReplaceOffsetValues(int GlobalRow, int NumEntries, double *Values, int *Indices);
01154   int SumIntoOffsetValues(int GlobalRow, int NumEntries, double *Values, int *Indices);
01155   void UpdateImportVector(int NumVectors) const;
01156   void UpdateExportVector(int NumVectors) const;
01157   void GeneralMV(double * x, double * y) const;
01158   void GeneralMTV(double * x, double * y) const;
01159   void GeneralMM(double ** X, int LDX, double ** Y, int LDY, int NumVectors) const;
01160   void GeneralMTM(double ** X, int LDX, double ** Y, int LDY, int NumVectors) const;
01161   void GeneralSV(bool Upper, bool Trans, bool UnitDiagonal, double * x, double * y) const;
01162   void GeneralSM(bool Upper, bool Trans, bool UnitDiagonal, double ** X, int LDX, double ** Y, int LDY, int NumVectors) const;
01163 
01164   void SetStaticGraph(bool Flag) {StaticGraph_ = Flag;}
01165 
01166   int CheckSizes(const Epetra_SrcDistObject& A);
01167 
01168   int CopyAndPermute(const Epetra_SrcDistObject& Source,
01169                      int NumSameIDs, 
01170                      int NumPermuteIDs,
01171                      int* PermuteToLIDs,
01172                      int* PermuteFromLIDs,
01173                      const Epetra_OffsetIndex * Indexor);
01174   int CopyAndPermuteCrsMatrix(const Epetra_CrsMatrix& A,
01175                               int NumSameIDs, 
01176                               int NumPermuteIDs,
01177                               int* PermuteToLIDs,
01178                               int* PermuteFromLIDs,
01179                               const Epetra_OffsetIndex * Indexor);
01180   int CopyAndPermuteRowMatrix(const Epetra_RowMatrix& A,
01181                               int NumSameIDs, 
01182                               int NumPermuteIDs,
01183                               int* PermuteToLIDs,
01184                               int* PermuteFromLIDs,
01185                               const Epetra_OffsetIndex * Indexor);
01186   
01187   int PackAndPrepare(const Epetra_SrcDistObject& Source,
01188                      int NumExportIDs,
01189                      int* ExportLIDs,
01190                      int& LenExports,
01191                      char*& Exports,
01192                      int& SizeOfPacket,
01193                      int* Sizes,
01194                      bool& VarSizes,
01195                      Epetra_Distributor& Distor);
01196 
01197   int UnpackAndCombine(const Epetra_SrcDistObject& Source, 
01198                        int NumImportIDs,
01199                        int* ImportLIDs, 
01200                        int LenImports,
01201                        char* Imports,
01202                        int& SizeOfPacket, 
01203                        Epetra_Distributor& Distor,
01204                        Epetra_CombineMode CombineMode,
01205                        const Epetra_OffsetIndex * Indexor);
01206 
01208   int SortEntries();
01209 
01211   bool Sorted() const {return(Graph_.Sorted());}
01212 
01214   int MergeRedundantEntries();
01215 
01217   bool NoRedundancies() const {return(Graph_.NoRedundancies());}
01218     
01219   void DeleteMemory();
01220 
01221   Epetra_CrsGraph Graph_;
01222   bool Allocated_;
01223   bool StaticGraph_;
01224   bool UseTranspose_;
01225   bool constructedWithFilledGraph_;
01226   bool matrixFillCompleteCalled_;
01227   bool StorageOptimized_;
01228 
01229   double** Values_;
01230   double* All_Values_;
01231   mutable double NormInf_;
01232   mutable double NormOne_;
01233   mutable double NormFrob_;
01234 
01235   int NumMyRows_;
01236   mutable Epetra_MultiVector* ImportVector_;
01237   mutable Epetra_MultiVector* ExportVector_;
01238 
01239   Epetra_DataAccess CV_;
01240 
01241   bool squareFillCompleteCalled_;
01242  private:
01243 
01244   // These are the pre-5.0 versions of solve.  They are still faster that generic 5.0 solves, so we keep them around
01245   int Solve1(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const;
01246   int Solve1(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
01247 
01248 };
01249 #endif /* EPETRA_CRSMATRIX_H */

Generated on Wed May 12 21:41:05 2010 for Epetra Package Browser (Single Doxygen Collection) by  doxygen 1.4.7