Epetra_VbrMatrix.h

Go to the documentation of this file.

/*
//@HEADER
// ************************************************************************
//
//               Epetra: Linear Algebra Services Package
//                 Copyright 2011 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER
*/

#ifndef EPETRA_VBRMATRIX_H
#define EPETRA_VBRMATRIX_H

// FIXME long long : whole file

#include <Epetra_ConfigDefs.h>
#include <Epetra_DistObject.h>
#include <Epetra_CompObject.h>
#include <Epetra_BLAS.h>
#include <Epetra_RowMatrix.h>
#include <Epetra_Operator.h>
#include <Epetra_CrsGraph.h>
#include <Epetra_SerialDenseMatrix.h>
class Epetra_BlockMap;
class Epetra_Map;
class Epetra_Import;
class Epetra_Export;
class Epetra_Vector;
class Epetra_MultiVector;


class EPETRA_LIB_DLL_EXPORT Epetra_VbrMatrix : public Epetra_DistObject,
          public Epetra_CompObject,
          public Epetra_BLAS,
          public virtual Epetra_RowMatrix {
 public:



  Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, int *NumBlockEntriesPerRow);


  Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, int NumBlockEntriesPerRow);


  Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, const Epetra_BlockMap& ColMap, int *NumBlockEntriesPerRow);


  Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_BlockMap& RowMap, const Epetra_BlockMap& ColMap, int NumBlockEntriesPerRow);


  Epetra_VbrMatrix(Epetra_DataAccess CV, const Epetra_CrsGraph & Graph);

  Epetra_VbrMatrix(const Epetra_VbrMatrix & Matrix);

  virtual ~Epetra_VbrMatrix();



  Epetra_VbrMatrix& operator=(const Epetra_VbrMatrix& src);


    int PutScalar(double ScalarConstant);


    int Scale(double ScalarConstant);


    int DirectSubmitBlockEntry(int GlobalBlockRow, int GlobalBlockCol,
                               const double *values, int LDA,
                               int NumRows, int NumCols, bool sum_into);


    int BeginInsertGlobalValues(int BlockRow,
        int NumBlockEntries,
        int * BlockIndices);


    int BeginInsertMyValues(int BlockRow, int NumBlockEntries, int * BlockIndices);


    int BeginReplaceGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices);


    int BeginReplaceMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices);


    int BeginSumIntoGlobalValues(int BlockRow, int NumBlockEntries, int *BlockIndices);


    int BeginSumIntoMyValues(int BlockRow, int NumBlockEntries, int *BlockIndices);

    /* Submit a block entry that will recorded in the block row that was initiated by one of the
       Begin routines listed above.  Once a one of the following routines: BeginInsertGlobalValues(),
       BeginInsertMyValues(), BeginReplaceGlobalValues(), BeginReplaceMyValues(), BeginSumIntoGlobalValues(),
       BeginSumIntoMyValues(), you \e must call SubmitBlockEntry() NumBlockEntries times to register the values
       corresponding to the block indices passed in to the Begin routine.  If the Epetra_VbrMatrix constuctor
       was called in Copy mode, the values will be copied.  However, no copying will be done until the EndSubmitEntries()
       function is call to complete submission of the current block row.  If the constructor was called in View mode, all
       block entries passed via SubmitBlockEntry() will not be copied, but a pointer will be set to point to the argument Values
       that was passed in by the user.

       For performance reasons, SubmitBlockEntry() does minimal processing of data.  Any processing that can be
       delayed is performed in EndSubmitEntries().

    \param In
           Values - The starting address of the values.
    \param In
           LDA - The stride between successive columns of Values.
    \param In
           NumRows - The number of rows passed in.
    \param In
           NumCols - The number of columns passed in.

    \return Integer error code, set to 0 if successful.
    */
    int SubmitBlockEntry(double *Values, int LDA, int NumRows, int NumCols);

    /* Submit a block entry that will recorded in the block row that was initiated by one of the
       Begin routines listed above.  Once a one of the following routines: BeginInsertGlobalValues(),
       BeginInsertMyValues(), BeginReplaceGlobalValues(), BeginReplaceMyValues(), BeginSumIntoGlobalValues(),
       BeginSumIntoMyValues(), you \e must call SubmitBlockEntry() NumBlockEntries times to register the values
       corresponding to the block indices passed in to the Begin routine.  If the Epetra_VbrMatrix constuctor
       was called in Copy mode, the values will be copied.  However, no copying will be done until the EndSubmitEntries()
       function is call to complete submission of the current block row.  If the constructor was called in View mode, all
       block entries passed via SubmitBlockEntry() will not be copied, but a pointer will be set to point to the argument Values
       that was passed in by the user.

       For performance reasons, SubmitBlockEntry() does minimal processing of data.  Any processing that can be
       delayed is performed in EndSubmitEntries().

    \param In
           Mat - Preformed dense matrix block.

    \return Integer error code, set to 0 if successful.
    */
    int SubmitBlockEntry( Epetra_SerialDenseMatrix &Mat );


    int EndSubmitEntries();


    int ReplaceDiagonalValues(const Epetra_Vector & Diagonal);

    /* This version of FillComplete assumes that the domain and range
       distributions are identical to the matrix row distributions.
       \return error code, 0 if successful. Returns a positive warning code of 3
        if the matrix is rectangular (meaning that the other overloading of
        FillComplete should have been called, with differen domain-map and
        range-map specified).
    */
    int FillComplete();

    /* This version of FillComplete requires the explicit specification of the domain
       and range distribution maps.  These maps are used for importing and exporting vector
       and multi-vector elements that are needed for distributed matrix computations.  For
       example, to compute y = Ax in parallel, we would specify the DomainMap as the distribution
       of the vector x and the RangeMap as the distribution of the vector y.
    \param In
           DomainMap - Map that describes the distribution of vector and multi-vectors in the
                 matrix domain.
    \param In
           RangeMap - Map that describes the distribution of vector and multi-vectors in the
                 matrix range.

    \return error code, 0 if successful. positive warning code of 2 if it is detected that the
    matrix-graph got out of sync since this matrix was constructed (for instance if
    graph.FillComplete() was called by another matrix that shares the graph)
    */
    int FillComplete(const Epetra_BlockMap& DomainMap, const Epetra_BlockMap& RangeMap);

    bool Filled() const {return(Graph_->Filled());};




    int ExtractGlobalBlockRowPointers(int BlockRow, int MaxNumBlockEntries,
              int & RowDim,  int & NumBlockEntries,
              int * BlockIndices,
              Epetra_SerialDenseMatrix ** & Values) const;


    int ExtractMyBlockRowPointers(int BlockRow, int MaxNumBlockEntries,
          int & RowDim, int & NumBlockEntries,
          int * BlockIndices,
          Epetra_SerialDenseMatrix** & Values) const;


    int BeginExtractGlobalBlockRowCopy(int BlockRow, int MaxNumBlockEntries,
               int & RowDim,  int & NumBlockEntries,
               int * BlockIndices, int * ColDims) const;


    int BeginExtractMyBlockRowCopy(int BlockRow, int MaxNumBlockEntries,
           int & RowDim, int & NumBlockEntries,
           int * BlockIndices, int * ColDims) const;


    int ExtractEntryCopy(int SizeOfValues, double * Values, int LDA, bool SumInto) const;


    int BeginExtractGlobalBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
               int * & BlockIndices) const;


    int BeginExtractMyBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
           int * & BlockIndices) const;



    int ExtractEntryView(Epetra_SerialDenseMatrix* & entry) const;


    int ExtractGlobalBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
          int * & BlockIndices,
          Epetra_SerialDenseMatrix** & Values) const;


    int ExtractMyBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
            int * & BlockIndices,
            Epetra_SerialDenseMatrix** & Values) const;



    int ExtractDiagonalCopy(Epetra_Vector & Diagonal) const;


    int BeginExtractBlockDiagonalCopy(int MaxNumBlockDiagonalEntries,
              int & NumBlockDiagonalEntries, int * RowColDims ) const;

    int ExtractBlockDiagonalEntryCopy(int SizeOfValues, double * Values, int LDA, bool SumInto) const;


    int BeginExtractBlockDiagonalView(int & NumBlockDiagonalEntries, int * & RowColDims ) const;


    int ExtractBlockDiagonalEntryView(double * & Values, int & LDA) const;





    int Multiply1(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) const;


    int Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;


    int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const;


    int Solve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;



    int InvRowSums(Epetra_Vector& x) const;


    int LeftScale(const Epetra_Vector& x);


    int InvColSums(Epetra_Vector& x) const ;


    int RightScale(const Epetra_Vector& x);



    int OptimizeStorage();

    bool StorageOptimized() const {return(StorageOptimized_);};

    bool IndicesAreGlobal() const {return(Graph_->IndicesAreGlobal());};

    bool IndicesAreLocal() const {return(Graph_->IndicesAreLocal());};

    bool IndicesAreContiguous() const {return(Graph_->IndicesAreContiguous());};

    bool LowerTriangular() const {return(Graph_->LowerTriangular());};

    bool UpperTriangular() const {return(Graph_->UpperTriangular());};

    bool NoDiagonal() const {return(Graph_->NoDiagonal());};




    /* Returns the quantity \f$ \| A \|_\infty\f$ such that
       \f[\| A \|_\infty = \max_{1\lei\lem} \sum_{j=1}^n |a_{ij}| \f].
     \warning The NormInf() method will not properly calculate the infinity norm for a matrix that has entries that are
     replicated on multiple processors.  */
    double NormInf() const;

    /* Returns the quantity \f$ \| A \|_1\f$ such that
       \f[\| A \|_1 = \max_{1\lej\len} \sum_{i=1}^m |a_{ij}| \f].
     \warning The NormOne() method will not properly calculate the one norm for a matrix that has entries that are
    */
    double NormOne() const;

    /* Returns the quantity \f[ \| A \|_{Frobenius} = \sqrt{\sum_{i=1}^m \sum_{j=1}^n\|a_{ij}\|^2}\f]
       \warning the NormFrobenius() method will not properly calculate the frobenius norm for a matrix that
       has entries which are replicated on multiple processors. In that case, the returned
       norm will be larger than the true norm.
    */
    double NormFrobenius() const;

    int MaxRowDim() const {return(Graph_->MaxRowDim());};

    int MaxColDim() const {return(Graph_->MaxColDim());};

    int GlobalMaxRowDim() const {return(Graph_->GlobalMaxRowDim());};

    int GlobalMaxColDim() const {return(Graph_->GlobalMaxColDim());};

    int NumMyRows() const {return(Graph_->NumMyRows());};
    int NumMyCols() const {return(Graph_->NumMyCols());};

    int NumMyNonzeros() const {return(Graph_->NumMyNonzeros());};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalRows() const {return(Graph_->NumGlobalRows());};
#endif
    long long NumGlobalRows64() const {return(Graph_->NumGlobalRows64());};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalCols() const {return(Graph_->NumGlobalCols());};
#endif
    long long NumGlobalCols64() const {return(Graph_->NumGlobalCols64());};

    /*
     Note that if maps are defined such that some nonzeros appear on
     multiple processors, then those nonzeros will be counted multiple
     times. If the user wishes to assemble a matrix from overlapping
     submatrices, they can use Epetra_FEVbrMatrix.
    */
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalNonzeros() const {return(Graph_->NumGlobalNonzeros());};
#endif
    long long NumGlobalNonzeros64() const {return(Graph_->NumGlobalNonzeros64());};

    int NumMyBlockRows() const {return(Graph_->NumMyBlockRows());};

    int NumMyBlockCols() const {return(Graph_->NumMyBlockCols());};

    int NumMyBlockEntries() const {return(Graph_->NumMyEntries());};

    int NumMyBlockDiagonals() const {return(Graph_->NumMyBlockDiagonals());};

    int NumMyDiagonals() const {return(Graph_->NumMyDiagonals());};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalBlockRows() const {return(Graph_->NumGlobalBlockRows());};
#endif
    long long NumGlobalBlockRows64() const {return(Graph_->NumGlobalBlockRows64());};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalBlockCols() const {return(Graph_->NumGlobalBlockCols());};
#endif
    long long NumGlobalBlockCols64() const {return(Graph_->NumGlobalBlockCols64());};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalBlockEntries() const {return(Graph_->NumGlobalEntries());};
#endif
    long long NumGlobalBlockEntries64() const {return(Graph_->NumGlobalEntries64());};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalBlockDiagonals() const {return(Graph_->NumGlobalBlockDiagonals());};
#endif
    long long NumGlobalBlockDiagonals64() const {return(Graph_->NumGlobalBlockDiagonals64());};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalDiagonals() const {return(Graph_->NumGlobalDiagonals());};
#endif
    long long NumGlobalDiagonals64() const {return(Graph_->NumGlobalDiagonals64());};

    int NumGlobalBlockEntries(int Row) const {return(Graph_->NumGlobalIndices(Row));};

    int NumAllocatedGlobalBlockEntries(int Row) const{return(Graph_->NumAllocatedGlobalIndices(Row));};

    int MaxNumBlockEntries() const {return(Graph_->MaxNumIndices());};

    int GlobalMaxNumBlockEntries() const {return(Graph_->GlobalMaxNumIndices());};

    int NumMyBlockEntries(int Row) const {return(Graph_->NumMyIndices(Row));};

    int NumAllocatedMyBlockEntries(int Row) const {return(Graph_->NumAllocatedMyIndices(Row));};


    int MaxNumNonzeros() const {return(Graph_->MaxNumNonzeros());};


    int GlobalMaxNumNonzeros() const {return(Graph_->GlobalMaxNumNonzeros());};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int  IndexBase() const {
    if(RowMap().GlobalIndicesInt())
      return (int) IndexBase64();
    throw "Epetra_VbrMatrix::IndexBase: GlobalIndices not int.";
  }
#endif
  long long  IndexBase64() const {return(Graph_->IndexBase64());};

    const Epetra_CrsGraph & Graph() const {return(*Graph_);};

    const Epetra_Import * Importer() const {return(Graph_->Importer());};

    const Epetra_Export * Exporter() const {return(Graph_->Exporter());};

    const Epetra_BlockMap & DomainMap() const {return(Graph_->DomainMap());};

    const Epetra_BlockMap & RangeMap() const  {return(Graph_->RangeMap());};

    const Epetra_BlockMap & RowMap() const {return(Graph_->RowMap());};

    const Epetra_BlockMap & ColMap() const {return(Graph_->ColMap());};


    const Epetra_Comm & Comm() const {return(Graph_->Comm());};



#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int LRID( int GRID_in) const {return(Graph_->LRID(GRID_in));};
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
    int LRID( long long GRID_in) const {return(Graph_->LRID(GRID_in));};
#endif

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int GRID(int LRID_in) const {return(Graph_->GRID(LRID_in));};
#endif
    long long GRID64( int LRID_in) const {return(Graph_->GRID64(LRID_in));};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int LCID( int GCID_in) const {return(Graph_->LCID(GCID_in));};
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
    int LCID( long long GCID_in) const {return(Graph_->LCID(GCID_in));};
#endif

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int GCID(int LCID_in) const {return(Graph_->GCID(LCID_in));};
#endif
    long long GCID64( int LCID_in) const {return(Graph_->GCID64(LCID_in));};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    bool  MyGRID(int GRID_in) const {return(Graph_->MyGRID(GRID_in));};
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
    bool  MyGRID(long long GRID_in) const {return(Graph_->MyGRID(GRID_in));};
#endif

    bool  MyLRID(int LRID_in) const {return(Graph_->MyLRID(LRID_in));};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    bool  MyGCID(int GCID_in) const {return(Graph_->MyGCID(GCID_in));};
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
    bool  MyGCID(long long GCID_in) const {return(Graph_->MyGCID(GCID_in));};
#endif

    bool  MyLCID(int LCID_in) const {return(Graph_->MyLCID(LCID_in));};

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    bool MyGlobalBlockRow(int GID) const {return(Graph_->MyGlobalRow(GID));};
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
    bool MyGlobalBlockRow(long long GID) const {return(Graph_->MyGlobalRow(GID));};
#endif




  virtual void Print(std::ostream & os) const;



    const char * Label() const {return(Epetra_Object::Label());};


  int SetUseTranspose(bool UseTranspose_in) {UseTranspose_ = UseTranspose_in; return(0);};


  int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;


  int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;

    bool HasNormInf() const {return(true);};

    bool UseTranspose() const {return(UseTranspose_);};

    const Epetra_Map & OperatorDomainMap() const
    {
      if (!HavePointObjects_) GeneratePointObjects();
      if (UseTranspose()) return(*OperatorRangeMap_);
      else return(*OperatorDomainMap_);
    }

    const Epetra_Map & OperatorRangeMap() const
    {
      if (!HavePointObjects_) GeneratePointObjects();
      if (UseTranspose()) return(*OperatorDomainMap_);
      else return(*OperatorRangeMap_);
    }




    int ExtractGlobalRowCopy(int GlobalRow, int Length, int & NumEntries, double *Values, int * Indices) const;


    int ExtractMyRowCopy(int MyRow, int Length, int & NumEntries, double *Values, int * Indices) const;


    int NumMyRowEntries(int MyRow, int & NumEntries) const;

    int MaxNumEntries() const;

    const Epetra_BlockMap& Map() const { return Epetra_DistObject::Map(); }

    const Epetra_Map & RowMatrixRowMap() const
      { if (!HavePointObjects_) GeneratePointObjects(); return(*RowMatrixRowMap_); };

    const Epetra_Map & RowMatrixColMap() const
      { if (!HavePointObjects_) GeneratePointObjects(); return(*RowMatrixColMap_); };

    const Epetra_Import * RowMatrixImporter() const
      { if (!HavePointObjects_) GeneratePointObjects(); return(RowMatrixImporter_); };




    const Epetra_BlockMap & BlockImportMap() const {return(Graph_->ImportMap());};

    int TransformToLocal();

    int TransformToLocal(const Epetra_BlockMap* DomainMap, const Epetra_BlockMap* RangeMap);

 protected:
    void DeleteMemory();
    bool Allocated() const {return(Allocated_);};
    int SetAllocated(bool Flag) {Allocated_ = Flag; return(0);};
    Epetra_SerialDenseMatrix *** Values() const {return(Entries_);};

  // Internal utilities

  int DoMultiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
  int DoSolve(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;
  void InitializeDefaults();
  int Allocate();
  int BeginInsertValues(int BlockRow, int NumBlockEntries,
      int * BlockIndices, bool IndicesAreLocal);
  int BeginReplaceValues(int BlockRow, int NumBlockEntries,
       int *BlockIndices, bool IndicesAreLocal);
  int BeginSumIntoValues(int BlockRow, int NumBlockEntries,
       int *BlockIndices, bool IndicesAreLocal);
  int SetupForSubmits(int BlockRow, int NumBlockEntries, int * BlockIndices,
          bool IndicesAreLocal, Epetra_CombineMode SubmitMode);
  int EndReplaceSumIntoValues();
  int EndInsertValues();

  int CopyMat(double * A, int LDA, int NumRows, int NumCols,
        double * B, int LDB, bool SumInto) const;
  int BeginExtractBlockRowCopy(int BlockRow, int MaxNumBlockEntries,
             int & RowDim, int & NumBlockEntries,
             int * BlockIndices, int * ColDims,
             bool IndicesAreLocal) const;
  int SetupForExtracts(int BlockRow, int & RowDim, int NumBlockEntries,
           bool ExtractView, bool IndicesAreLocal) const;
  int ExtractBlockDimsCopy(int NumBlockEntries, int * ColDims) const;
  int ExtractBlockRowPointers(int BlockRow, int MaxNumBlockEntries,
            int & RowDim, int & NumBlockEntries,
            int * BlockIndices,
            Epetra_SerialDenseMatrix ** & Values,
            bool IndicesAreLocal) const;
  int BeginExtractBlockRowView(int BlockRow, int & RowDim, int & NumBlockEntries,
             int * & BlockIndices,
             bool IndicesAreLocal) const;
  int CopyMatDiag(double * A, int LDA, int NumRows, int NumCols,
      double * Diagonal) const;
  int ReplaceMatDiag(double * A, int LDA, int NumRows, int NumCols,
         double * Diagonal);

  //This BlockRowMultiply accepts Alpha and Beta arguments. It is called
  //from within the 'solve' methods.
  void BlockRowMultiply(bool TransA, int RowDim, int NumEntries,
      int * BlockIndices, int RowOff,
      int * FirstPointInElementList, int * ElementSizeList,
      double Alpha, Epetra_SerialDenseMatrix** As,
      double ** X, double Beta, double ** Y, int NumVectors) const;

  //This BlockRowMultiply doesn't accept Alpha and Beta arguments, instead it
  //assumes that they are both 1.0. It is called from within the 'Multiply'
  //methods.
  void BlockRowMultiply(bool TransA, int RowDim, int NumEntries,
      int * BlockIndices, int RowOff,
      int * FirstPointInElementList,
      int * ElementSizeList,
      Epetra_SerialDenseMatrix** As,
      double ** X, double ** Y, int NumVectors) const;
  //
  // Assumes Alpha=Beta=1 and works only on storage optimized matrices
  //
  void FastBlockRowMultiply(bool TransA, int RowDim, int NumEntries,
          int * BlockIndices, int RowOff,
          int * FirstPointInElementList,
          int * ElementSizeList,
          Epetra_SerialDenseMatrix** As,
          double ** X, double ** Y, int NumVectors) const;

  int InverseSums(bool DoRows, Epetra_Vector& x) const;
  int Scale(bool DoRows, const Epetra_Vector& x);
  void BlockRowNormInf(int RowDim, int NumEntries,
           Epetra_SerialDenseMatrix** As,
           double * Y) const;
  void BlockRowNormOne(int RowDim, int NumEntries, int * BlockRowIndices,
           Epetra_SerialDenseMatrix** As,
           int * ColFirstPointInElementList, double * x) const;
  void SetStaticGraph(bool Flag) {StaticGraph_ = Flag;};

  int CheckSizes(const Epetra_SrcDistObject& A);

  int CopyAndPermute(const Epetra_SrcDistObject & Source,
                     int NumSameIDs,
                     int NumPermuteIDs,
                     int * PermuteToLIDs,
                     int *PermuteFromLIDs,
                     const Epetra_OffsetIndex * Indexor,
                     Epetra_CombineMode CombineMode = Zero);

  int PackAndPrepare(const Epetra_SrcDistObject & Source,
                     int NumExportIDs,
                     int * ExportLIDs,
                     int & LenExports,
                     char * & Exports,
                     int & SizeOfPacket,
                     int * Sizes,
                     bool & VarSizes,
                     Epetra_Distributor & Distor);

  int UnpackAndCombine(const Epetra_SrcDistObject & Source,
                       int NumImportIDs,
                       int * ImportLIDs,
                       int LenImports,
                       char * Imports,
                       int & SizeOfPacket,
                       Epetra_Distributor & Distor,
                       Epetra_CombineMode CombineMode,
                       const Epetra_OffsetIndex * Indexor);

  int SortEntries();

  bool Sorted() const {return(Graph_->Sorted());};

  int MergeRedundantEntries();

  bool NoRedundancies() const {return(Graph_->NoRedundancies());};

  bool StaticGraph() const {return(StaticGraph_);};

  int GeneratePointObjects() const;
  int BlockMap2PointMap(const Epetra_BlockMap & BlockMap, Epetra_Map * & PointMap) const;
  int UpdateOperatorXY(const Epetra_MultiVector& X, const Epetra_MultiVector& Y) const;

  Epetra_CrsGraph * Graph_;
  bool Allocated_;
  bool StaticGraph_;
  bool UseTranspose_;
  bool constructedWithFilledGraph_;
  bool matrixFillCompleteCalled_;
  bool StorageOptimized_;

  int NumMyBlockRows_;

  Epetra_DataAccess CV_;


  int * NumBlockEntriesPerRow_;
  int * NumAllocatedBlockEntriesPerRow_;
  int ** Indices_;
  int * ElementSizeList_;
  int * FirstPointInElementList_;

  Epetra_SerialDenseMatrix ***Entries_;

  double *All_Values_Orig_;
  double *All_Values_;

  mutable double NormInf_;
  mutable double NormOne_;
  mutable double NormFrob_;

  mutable Epetra_MultiVector * ImportVector_;
  mutable Epetra_MultiVector * ExportVector_;

  // State variables needed for constructing matrix entry-by-entry
  mutable int *TempRowDims_;
  mutable Epetra_SerialDenseMatrix **TempEntries_;
  mutable int LenTemps_;
  mutable int CurBlockRow_;
  mutable int CurNumBlockEntries_;
  mutable int * CurBlockIndices_;
  mutable int CurEntry_;
  mutable bool CurIndicesAreLocal_;
  mutable Epetra_CombineMode CurSubmitMode_;

  // State variables needed for extracting entries
  mutable int CurExtractBlockRow_;
  mutable int CurExtractEntry_;
  mutable int CurExtractNumBlockEntries_;
  mutable bool CurExtractIndicesAreLocal_;
  mutable bool CurExtractView_;
  mutable int CurRowDim_;

  // State variable for extracting block diagonal entries
  mutable int CurBlockDiag_;

  // Maps and importer that support the Epetra_RowMatrix interface
  mutable Epetra_Map * RowMatrixRowMap_;
  mutable Epetra_Map * RowMatrixColMap_;
  mutable Epetra_Import * RowMatrixImporter_;

  // Maps that support the Epetra_Operator interface
  mutable Epetra_Map * OperatorDomainMap_;
  mutable Epetra_Map * OperatorRangeMap_;
  mutable Epetra_MultiVector * OperatorX_;
  mutable Epetra_MultiVector * OperatorY_;

  // bool to indicate if above four point maps and importer have already been created
  mutable bool HavePointObjects_;

  bool squareFillCompleteCalled_;
};

#endif /* EPETRA_VBRMATRIX_H */