Epetra_CrsMatrix.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_CRSMATRIX_H
#define EPETRA_CRSMATRIX_H

#if defined(Epetra_SHOW_DEPRECATED_WARNINGS)
#ifdef __GNUC__
#warning "The Epetra package is deprecated"
#endif
#endif



#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_Map.h"

#ifdef Epetra_ENABLE_CASK
#include "cask.h"
#endif

class Epetra_Map;
class Epetra_Import;
class Epetra_Export;
class Epetra_Vector;
class Epetra_MultiVector;
class Epetra_IntSerialDenseVector;


// Define this to see a complete dump a an Epetra_CrsMatrix::Multiply(...) call
//#define EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY

#ifdef EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY
extern bool Epetra_CrsMatrixTraceDumpMultiply;
#endif // EPETRA_CRS_MATRIX_TRACE_DUMP_MULTIPLY


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



  Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& RowMap, const int* NumEntriesPerRow, bool StaticProfile = false);


  Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& RowMap, int NumEntriesPerRow, bool StaticProfile = false);


  Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& RowMap, const Epetra_Map& ColMap, const int* NumEntriesPerRow, bool StaticProfile = false);


  Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_Map& RowMap, const Epetra_Map& ColMap, int NumEntriesPerRow, bool StaticProfile = false);


  Epetra_CrsMatrix(Epetra_DataAccess CV, const Epetra_CrsGraph& Graph);


  Epetra_CrsMatrix(const Epetra_CrsMatrix & SourceMatrix, const Epetra_Import & RowImporter, const Epetra_Map * DomainMap=0, const Epetra_Map * RangeMap=0, bool RestrictCommunicator = false);


  Epetra_CrsMatrix(const Epetra_CrsMatrix & SourceMatrix, const Epetra_Import & RowImporter, const Epetra_Import * DomainImporter, const Epetra_Map * DomainMap, const Epetra_Map * RangeMap, bool RestrictCommunicator);


  Epetra_CrsMatrix(const Epetra_CrsMatrix & SourceMatrix, const Epetra_Export & RowExporter, const Epetra_Map * DomainMap=0, const Epetra_Map * RangeMap=0, bool RestrictCommunicator = false);


  Epetra_CrsMatrix(const Epetra_CrsMatrix & SourceMatrix, const Epetra_Export & RowExporter, const Epetra_Export * DomainExporter, const Epetra_Map * DomainMap, const Epetra_Map * RangeMap, bool RestrictCommunicator);


  Epetra_CrsMatrix(const Epetra_CrsMatrix& Matrix);

  virtual ~Epetra_CrsMatrix();



  Epetra_CrsMatrix& operator=(const Epetra_CrsMatrix& src);


  int PutScalar(double ScalarConstant);


  int Scale(double ScalarConstant);


#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  virtual int InsertGlobalValues(int GlobalRow, int NumEntries, const double* Values, const int* Indices);
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  virtual int InsertGlobalValues(long long GlobalRow, int NumEntries, const double* Values, const long long* Indices);
#endif


#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  virtual int InsertGlobalValues(int GlobalRow, int NumEntries, double* Values, int* Indices);
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  virtual int InsertGlobalValues(long long GlobalRow, int NumEntries, double* Values, long long* Indices);
#endif


#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  virtual int ReplaceGlobalValues(int GlobalRow, int NumEntries, const double* Values, const int* Indices);
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  virtual int ReplaceGlobalValues(long long GlobalRow, int NumEntries, const double* Values, const long long* Indices);
#endif


#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  virtual int SumIntoGlobalValues(int GlobalRow, int NumEntries, const double* Values, const int* Indices);
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  virtual int SumIntoGlobalValues(long long GlobalRow, int NumEntries, const double* Values, const long long* Indices);
#endif


  int InsertMyValues(int MyRow, int NumEntries, const double* Values, const int* Indices);


  int InsertMyValues(int MyRow, int NumEntries, double* Values, int* Indices);


  int ReplaceMyValues(int MyRow, int NumEntries, const double* Values, const int* Indices);


  int SumIntoMyValues(int MyRow, int NumEntries, const double* Values, const int* Indices);


  int ReplaceDiagonalValues(const Epetra_Vector& Diagonal);





  /* This version of FillComplete assumes that the domain and range
     distributions are identical to the matrix row distributions.
    \param OptimizeDataStorage - (In) If true, storage will be packed for optimal performance.  Depending
           on how the matrix was constructed, optimizing the storage may have no impact on performance
     or one-time memory use, or may have a large impact.  If the user was careful in allocating memory
     for the matrix by setting StaticProfile to true in the matrix constructor, then no extra storage
     will be allocated in attempting to optimize storage.  If the user did not set StaticProfile to true,
     then optimizing the storage will temporarily use additional memory, will have a noticeable impact
     on performance and ultimately reduce the storage associated with the matrix.

     By default storage will be optimized.  If you cannot tolerate the increased temporary memory use,
     should set this value to false.

     \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(bool OptimizeDataStorage = true);

  /* 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 DomainMap - (In) Map that describes the distribution of vector and multi-vectors in the
    matrix domain.
    \param RangeMap - (In) Map that describes the distribution of vector and multi-vectors in the
    matrix range.

    \param OptimizeDataStorage - (In) If true, storage will be packed for optimal performance.  Depending
           on how the matrix was constructed, optimizing the storage may have no impact on performance
     or one-time memory use, or may have a large impact.  If the user was careful in allocating memory
     for the matrix by setting StaticProfile to true in the matrix constructor, then no extra storage
     will be allocated in attempting to optimize storage.  If the user did not set StaticProfile to true,
     then optimizing the storage will temporarily use additional memory, will have a noticeable impact
     on performance and ultimately reduce the storage associated with the matrix.

     By default storage will be optimized.  If you cannot tolerate the increased temporary memory use,
     should set this value to false.

    \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)

    \post IndicesAreLocal()==true
    */
  int FillComplete(const Epetra_Map& DomainMap, const Epetra_Map& RangeMap, bool OptimizeDataStorage = true);


  int OptimizeStorage();


  int MakeDataContiguous() {EPETRA_CHK_ERR(OptimizeStorage()); return(0);}




#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int ExtractGlobalRowCopy(int GlobalRow, int Length, int& NumEntries, double* Values, int* Indices) const;
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  int ExtractGlobalRowCopy(long long GlobalRow, int Length, int& NumEntries, double* Values, long long* Indices) const;
#endif


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


#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int ExtractGlobalRowCopy(int GlobalRow, int Length, int& NumEntries, double* Values) const;
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  int ExtractGlobalRowCopy(long long GlobalRow, int Length, int& NumEntries, double* Values) const;
#endif


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


  int ExtractDiagonalCopy(Epetra_Vector& Diagonal) const;


#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int ExtractGlobalRowView(int GlobalRow, int& NumEntries, double*& Values, int*& Indices) const;
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  int ExtractGlobalRowView(long long GlobalRow, int& NumEntries, double*& Values, long long*& Indices) const;
#endif


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


#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int ExtractGlobalRowView(int GlobalRow, int& NumEntries, double*& Values) const;
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  int ExtractGlobalRowView(long long GlobalRow, int& NumEntries, double*& Values) const;
#endif


  int ExtractMyRowView(int MyRow, int& NumEntries, double*& Values) const;




  int Multiply(bool TransA, const Epetra_Vector& x, Epetra_Vector& y) 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 Multiply1(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 InvRowMaxs(Epetra_Vector& x) const;


  int LeftScale(const Epetra_Vector& x);


  int InvColSums(Epetra_Vector& x) const;


  int InvColMaxs(Epetra_Vector& x) const;


  int RightScale(const Epetra_Vector& x);




  bool Filled() const {return(Graph_.Filled());}

  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
     replicated on multiple processors.
  */
  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;

  /*
    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_FECrsMatrix.
  */
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int  NumGlobalNonzeros() const {
      if(RowMap().GlobalIndicesInt())
        return (int) NumGlobalNonzeros64();
      throw "Epetra_CrsMatrix::NumGlobalNonzeros: GlobalIndices not int.";
    }
#endif
  long long NumGlobalNonzeros64() const {return(Graph_.NumGlobalNonzeros64());}

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int  NumGlobalRows() const {
      if(RowMap().GlobalIndicesInt())
        return (int) NumGlobalRows64();
      throw "Epetra_CrsMatrix::NumGlobalRows: GlobalIndices not int.";
    }
#endif
  long long NumGlobalRows64() const {return(Graph_.NumGlobalRows64());}

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int  NumGlobalCols() const {
      if(RowMap().GlobalIndicesInt())
        return (int) NumGlobalCols64();
      throw "Epetra_CrsMatrix::NumGlobalCols: GlobalIndices not int.";
    }
#endif
  long long NumGlobalCols64() const {return(Graph_.NumGlobalCols64());}

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int  NumGlobalDiagonals() const {
      if(RowMap().GlobalIndicesInt())
        return (int) NumGlobalDiagonals64();
      throw "Epetra_CrsMatrix::NumGlobalDiagonals: GlobalIndices not int.";
    }
#endif
  long long NumGlobalDiagonals64() const {return(Graph_.NumGlobalDiagonals64());}

  int NumMyNonzeros() const {return(Graph_.NumMyNonzeros());}

  int NumMyRows() const {return(Graph_.NumMyRows());}


  int NumMyCols() const {return(Graph_.NumMyCols());}


  int NumMyDiagonals() const {return(Graph_.NumMyDiagonals());}

  int NumGlobalEntries(long long Row) const {return(Graph_.NumGlobalIndices(Row));}

  int NumAllocatedGlobalEntries(int Row) const{return(Graph_.NumAllocatedGlobalIndices(Row));}


  int MaxNumEntries() const {return(Graph_.MaxNumIndices());}


  int GlobalMaxNumEntries() const {return(Graph_.GlobalMaxNumIndices());}

  int NumMyEntries(int Row) const {return(Graph_.NumMyIndices(Row));}

  int NumAllocatedMyEntries(int Row) const {return(Graph_.NumAllocatedMyIndices(Row));}

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


  bool StaticGraph() {return(StaticGraph_);}

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

  const Epetra_Map& RowMap() const {return((Epetra_Map &)Graph_.RowMap());}


  int ReplaceRowMap(const Epetra_BlockMap& newmap);


  bool HaveColMap() const {return(Graph_.HaveColMap());}


  int ReplaceColMap(const Epetra_BlockMap& newmap);


  int ReplaceDomainMapAndImporter(const Epetra_Map & NewDomainMap, const Epetra_Import * NewImporter);


  int RemoveEmptyProcessesInPlace(const Epetra_BlockMap * NewMap);


  const Epetra_Map& ColMap() const {return((Epetra_Map &) Graph_.ColMap());}


  const Epetra_Map& DomainMap() const {return((Epetra_Map &)Graph_.DomainMap());}


  const Epetra_Map& RangeMap() const  {return((Epetra_Map &)Graph_.RangeMap());}

  const Epetra_Import* Importer() const {return(Graph_.Importer());}

  const Epetra_Export* Exporter() const {return(Graph_.Exporter());}

  const Epetra_Comm& Comm() const {return(Epetra_DistObject::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

#if defined(EPETRA_NO_32BIT_GLOBAL_INDICES) && defined(EPETRA_NO_64BIT_GLOBAL_INDICES)
  // default implementation so that no compiler/linker error in case neither 32 nor 64
  // bit indices present.
  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 {
      if(RowMap().GlobalIndicesInt())
        return (int) GRID64(LRID_in);
      throw "Epetra_CrsMatrix::GRID: GlobalIndices not int.";
    }
#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 {
      if(RowMap().GlobalIndicesInt())
        return (int) GCID64(LCID_in);
      throw "Epetra_CrsMatrix::GCID: GlobalIndices not int.";
    }
#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 MyGlobalRow(int GID) const {return(Graph_.MyGlobalRow(GID));}
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  bool MyGlobalRow(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 {
    return(Epetra_CrsMatrix::Multiply(Epetra_CrsMatrix::UseTranspose(), X, Y));}


  int ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    return(Solve(UpperTriangular(), Epetra_CrsMatrix::UseTranspose(), NoDiagonal(), X, Y));}

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

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

  const Epetra_Map& OperatorDomainMap() const
        {
          if (UseTranspose()) return(RangeMap());
          else return(DomainMap());
        }

  const Epetra_Map& OperatorRangeMap() const
        {
          if (UseTranspose()) return(DomainMap());
          else return(RangeMap());
        }




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

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

  const Epetra_Map& RowMatrixRowMap() const {return(RowMap());}

  const Epetra_Map& RowMatrixColMap() const {return(ColMap());}

  const Epetra_Import* RowMatrixImporter() const {return(Importer());}





  inline double* operator[] (int Loc) {
    if (StorageOptimized()){ int * ind = Graph().IndexOffset(); return(All_Values_+ind[Loc]);}
    else return Values_[Loc];}
  inline double* operator[] (int Loc) const {
    if (StorageOptimized()){ int * ind = Graph().IndexOffset(); return(All_Values_+ind[Loc]);}
    else return Values_[Loc];}



  int ExtractCrsDataPointers(int *& IndexOffset, int *& Indices, double *& Values_in) const {
    if (StorageOptimized()) {
      IndexOffset = Graph().IndexOffset();
      Indices = Graph().All_Indices();
      Values_in  = All_Values();
      return (0);
    }
    else { IndexOffset = 0; Indices = 0; Values_in  = 0; return (-1);} }



    Epetra_IntSerialDenseVector& ExpertExtractIndexOffset();


    Epetra_IntSerialDenseVector& ExpertExtractIndices();


    double *& ExpertExtractValues() {return All_Values_;}



    int ExpertStaticFillComplete(const Epetra_Map & DomainMap,const Epetra_Map & RangeMap, const Epetra_Import * Importer=0, const Epetra_Export * Exporter=0, int NumMyDiagonals=-1);



    int ExpertMakeUniqueCrsGraphData();


  int SortGhostsAssociatedWithEachProcessor(bool Flag)  {Graph_.SortGhostsAssociatedWithEachProcessor(Flag); return(0);}



  const Epetra_Map& ImportMap() const {return((Epetra_Map&) Graph_.ImportMap());}

  int TransformToLocal();

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



 protected:
  bool Allocated() const {return(Allocated_);}
  int SetAllocated(bool Flag) {Allocated_ = Flag; return(0);}
  double** Values() const {
    if (StorageOptimized()) throw ReportError("This method: double** Values() cannot be called when StorageOptimized()==true", -1);
    else return(Values_);}
  double* All_Values() const {
    if (!StorageOptimized()) throw ReportError("This method: double* All_Values()cannot be called when StorageOptimized()==false", -1);
    else return(All_Values_);}
  double* Values(int LocalRow) const {
    if (StorageOptimized())
      if (Graph().StorageOptimized())
  return(All_Values_+Graph().IndexOffset()[LocalRow]);
      else throw ReportError("This method: double* Values()cannot be called when StorageOptimized()==true and Graph().StorageOptimized()==false", -1);
    else return(Values_[LocalRow]);}

  void InitializeDefaults();
  int Allocate();

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  int InsertValues(int LocalRow, int NumEntries, double* Values, int* Indices);
  int InsertValues(int LocalRow, int NumEntries, const double* Values, const int* Indices);
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  int InsertValues(int LocalRow, int NumEntries, double* Values, long long* Indices);
  int InsertValues(int LocalRow, int NumEntries, const double* Values, const long long* Indices);
#endif

  int InsertOffsetValues(long long GlobalRow, int NumEntries, double *Values, int *Indices);
  int InsertOffsetValues(long long GlobalRow, int NumEntries, const double *Values, const int *Indices);
  int ReplaceOffsetValues(long long GlobalRow, int NumEntries, const double *Values, const int *Indices);
  int SumIntoOffsetValues(long long GlobalRow, int NumEntries, const double *Values, const int *Indices);
  void UpdateImportVector(int NumVectors) const;
  void UpdateExportVector(int NumVectors) const;
  void GeneralMV(double * x, double * y) const;
  void GeneralMTV(double * x, double * y) const;
  void GeneralMM(double ** X, int LDX, double ** Y, int LDY, int NumVectors) const;
  void GeneralMTM(double ** X, int LDX, double ** Y, int LDY, int NumVectors) const;
  void GeneralSV(bool Upper, bool Trans, bool UnitDiagonal, double * x, double * y) const;
  void GeneralSM(bool Upper, bool Trans, bool UnitDiagonal, double ** X, int LDX, double ** Y, int LDY, int NumVectors) 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 CopyAndPermuteCrsMatrix(const Epetra_CrsMatrix& A,
                              int NumSameIDs,
                              int NumPermuteIDs,
                              int* PermuteToLIDs,
                              int* PermuteFromLIDs,
                              const Epetra_OffsetIndex * Indexor,
                              Epetra_CombineMode CombineMode);
  int CopyAndPermuteRowMatrix(const Epetra_RowMatrix& A,
                              int NumSameIDs,
                              int NumPermuteIDs,
                              int* PermuteToLIDs,
                              int* PermuteFromLIDs,
                              const Epetra_OffsetIndex * Indexor,
                              Epetra_CombineMode CombineMode);

  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());}

  void DeleteMemory();

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

  double** Values_;
  int* Values_alloc_lengths_;
  double* All_Values_;
  mutable double NormInf_;
  mutable double NormOne_;
  mutable double NormFrob_;

  int NumMyRows_;
  mutable Epetra_MultiVector* ImportVector_;
  mutable Epetra_MultiVector* ExportVector_;

  Epetra_DataAccess CV_;

  bool squareFillCompleteCalled_;
#ifdef Epetra_ENABLE_CASK
  caskHandle_t cask;
#endif
 private:

  // These are the pre-5.0 versions of solve.  They are still faster that generic 5.0 solves, so we keep them around
  int Solve1(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_Vector& x, Epetra_Vector& y) const;
  int Solve1(bool Upper, bool Trans, bool UnitDiagonal, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const;

private:

  template<typename int_type>
  int TInsertGlobalValues(int_type Row, int NumEntries, const double* values, const int_type* Indices);

  template<typename int_type>
  int TInsertGlobalValues(int_type Row, int NumEntries, double* values, int_type* Indices);

  template<typename int_type>
  int InsertValues(int Row, int NumEntries, const double* values, const int_type* Indices);

  template<typename int_type>
  int InsertValues(int Row, int NumEntries, double* values, int_type* Indices);

  template<typename int_type>
  int TReplaceGlobalValues(int_type Row, int NumEntries, const double * srcValues, const int_type *Indices);

  template<typename int_type>
  int TSumIntoGlobalValues(int_type Row, int NumEntries, const double * srcValues, const int_type *Indices);

  template<typename int_type>
  int ExtractGlobalRowCopy(int_type Row, int Length, int & NumEntries, double * values, int_type * Indices) const;

  template<typename int_type>
  int ExtractGlobalRowCopy(int_type Row, int Length, int & NumEntries, double * values) const;

  template<typename int_type>
  int ExtractGlobalRowView(int_type Row, int & NumEntries, double *& values, int_type *& Indices) const;

  template<typename int_type>
  int ExtractGlobalRowView(int_type Row, int & NumEntries, double *& values) const;

  template<typename int_type>
    int TCopyAndPermuteCrsMatrix(const Epetra_CrsMatrix& A,
                              int NumSameIDs,
                              int NumPermuteIDs,
                              int* PermuteToLIDs,
                              int* PermuteFromLIDs,
                              const Epetra_OffsetIndex * Indexor,
                              Epetra_CombineMode CombineMode);

  template<typename int_type>
    int TCopyAndPermuteRowMatrix(const Epetra_RowMatrix& A,
                              int NumSameIDs,
                              int NumPermuteIDs,
                              int* PermuteToLIDs,
                              int* PermuteFromLIDs,
                              const Epetra_OffsetIndex * Indexor,
                              Epetra_CombineMode CombineMode);

  template<typename int_type>
    int TUnpackAndCombine(const Epetra_SrcDistObject& Source,
                       int NumImportIDs,
                       int* ImportLIDs,
                       int LenImports,
                       char* Imports,
                       int& SizeOfPacket,
                       Epetra_Distributor& Distor,
                       Epetra_CombineMode CombineMode,
                       const Epetra_OffsetIndex * Indexor);

  // Used for fused[import|export] constructors
  template<class TransferType>
  void FusedTransfer(const Epetra_CrsMatrix & SourceMatrix,
          const TransferType & RowTransfer,
          const TransferType* DomainTransfer,
          const Epetra_Map * DomainMap,
          const Epetra_Map * RangeMap,
          bool RestrictCommunicator);

 public:

  void FusedImport(const Epetra_CrsMatrix & SourceMatrix,
        const Epetra_Import & RowImporter,
        const Epetra_Map * DomainMap,
        const Epetra_Map * RangeMap,
        bool RestrictCommunicator);

  void FusedExport(const Epetra_CrsMatrix & SourceMatrix,
        const Epetra_Export & RowExporter,
        const Epetra_Map * DomainMap,
        const Epetra_Map * RangeMap,
        bool RestrictCommunicator);

  void FusedImport(const Epetra_CrsMatrix & SourceMatrix,
        const Epetra_Import & RowImporter,
        const Epetra_Import * DomainImporter,
        const Epetra_Map * DomainMap,
        const Epetra_Map * RangeMap,
        bool RestrictCommunicator);

  void FusedExport(const Epetra_CrsMatrix & SourceMatrix,
        const Epetra_Export & RowExporter,
        const Epetra_Export * DomainExporter,
        const Epetra_Map * DomainMap,
        const Epetra_Map * RangeMap,
        bool RestrictCommunicator);


};
#endif /* EPETRA_CRSMATRIX_H */