EpetraExt_PETScAIJMatrix.h

Go to the documentation of this file.

//@HEADER
// ***********************************************************************
//
//     EpetraExt: Epetra Extended - 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

/*#############################################################################
# CVS File Information
#    Current revision: $Revision$
#    Last modified:    $Date$
#    Modified by:      $Author$
#############################################################################*/

#ifndef _EPETRAEXT_PETSCAIJMATRIX_H_
#define _EPETRAEXT_PETSCAIJMATRIX_H_

#include "Epetra_Object.h"
#include "Epetra_CompObject.h"
#include "Epetra_RowMatrix.h"
#include "Epetra_Map.h"
#ifdef HAVE_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_SerialComm.h"
#endif
//Petsc headers.
//Note: Petsc internally hard-codes paths to headers, relative to the PETSC home
//      directory.  This means that --with-incdirs must contain the full path(s)
//      to the header below plus the PETSc home directory.
#include "src/mat/impls/aij/mpi/mpiaij.h"

class Epetra_Import;
class Epetra_Export;
class Epetra_Vector;
class Epetra_MultiVector;


class Epetra_PETScAIJMatrix: public Epetra_Object, public Epetra_CompObject, public virtual Epetra_RowMatrix  {
      
 public:



  Epetra_PETScAIJMatrix(Mat Amat);

  virtual ~Epetra_PETScAIJMatrix();
  



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


    int ExtractDiagonalCopy(Epetra_Vector & Diagonal) const;




    int Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& 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);




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

    bool LowerTriangular() const {return(false);};

    bool UpperTriangular() const {return(false);};

  


    Mat Amat() const {return(Amat_);};

    /* Returns the quantity \f$ \| A \|_\infty\f$ such that
       \f[\| A \|_\infty = \max_{1\lei\lem} \sum_{j=1}^n |a_{ij}| \f].
    */ 
    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].
    */ 
    double NormOne() const;

#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
    int NumGlobalNonzeros() const {return(NumGlobalNonzeros_);};

    int NumGlobalRows() const {return(OperatorRangeMap().NumGlobalPoints());};

    int NumGlobalCols() const {return(OperatorDomainMap().NumGlobalPoints());};

    int NumGlobalDiagonals() const{return(OperatorDomainMap().NumGlobalPoints());};
#endif
    long long NumGlobalNonzeros64() const {return(NumGlobalNonzeros_);};

    long long NumGlobalRows64() const {return(OperatorRangeMap().NumGlobalPoints64());};

    long long NumGlobalCols64() const {return(OperatorDomainMap().NumGlobalPoints64());};

    long long NumGlobalDiagonals64() const{return(OperatorDomainMap().NumGlobalPoints64());};
    
    int NumMyNonzeros() const {return(NumMyNonzeros_);};

    int NumMyRows() const {return(OperatorRangeMap().NumMyPoints());};

    int NumMyCols() const {return(RowMatrixColMap().NumMyPoints());};

    int NumMyDiagonals() const {return(OperatorRangeMap().NumMyPoints());};

    const Epetra_Map & OperatorDomainMap() const {return(*DomainMap_);};

    const Epetra_Map & OperatorRangeMap() const  {return(*DomainMap_);};

    const Epetra_BlockMap& Map() const {return(RowMatrixRowMap());}

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

    const Epetra_Map & RowMatrixColMap() const {return(*ColMap_);};

    virtual const Epetra_Import * RowMatrixImporter() const {return(Importer_);};

    const Epetra_Comm & Comm() const {return(*Comm_);};
  
  


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



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

  int SetUseTranspose(bool UseTranspose)
    {(void)UseTranspose; return(-1);}


  int Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    return(Epetra_PETScAIJMatrix::Multiply(Epetra_PETScAIJMatrix::UseTranspose(), X, Y));};


  int ApplyInverse(const Epetra_MultiVector& X,
                   Epetra_MultiVector& Y) const
    {(void)X; (void)Y; return(-1);}

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

    virtual bool UseTranspose() const {return(false);}




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

    int MaxNumEntries() const;

 private:

    int GetRow(int Row) const;
    Mat Amat_;                //general PETSc matrix type
    mutable double * Values_;
    mutable int * Indices_;
    mutable int MaxNumEntries_;
    
#ifdef HAVE_MPI
    Epetra_MpiComm * Comm_;
#else
    Epetra_SerialComm * Comm_;
#endif  
    Epetra_Map * DomainMap_;
    Epetra_Map * ColMap_;
    Epetra_Import * Importer_;
    mutable Epetra_MultiVector * ImportVector_;
 
    double NumGlobalNonzeros_;
    int NumMyNonzeros_;
    int NumMyRows_;
    int NumGlobalRows_;
    int NumMyCols_;
    int PetscRowStart_;
    int PetscRowEnd_;
    enum petscMatrixType {PETSC_SEQ_AIJ, PETSC_MPI_AIJ};
    MatType  MatType_;         //really const char*
    mutable double NormInf_;
    mutable double NormOne_;

    
  //FIXME we need a copy ctor
  //Epetra_PETScAIJMatrix(const Epetra_PETScAIJMatrix & Matrix) {(void)Matrix;}
};
#endif /* _EPETRAEXT_PETSCAIJMATRIX_H_ */