Amesos_Klu.h

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                Amesos: Direct Sparse Solver Package
//                 Copyright (2004) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER

 
#ifndef AMESOS_KLU_H
#define AMESOS_KLU_H
 
#include "Amesos_ConfigDefs.h"
#include "Amesos_BaseSolver.h"
#include "Amesos_NoCopiable.h"
#include "Amesos_Utils.h"
#include "Amesos_Time.h"
#include "Amesos_Status.h"
#include "Amesos_Control.h"
#include "Epetra_LinearProblem.h"
#include "Epetra_Time.h"
#include "Epetra_Import.h"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_Comm.h"
#endif
#include "Epetra_CrsGraph.h"
#include "Epetra_CrsMatrix.h"
 
// class EpetraExt::MultiVector_Reindex ;
// class EpetraExt::CrsMatrix_Reindex ;

 
// Amesos_Klu_Pimpl contains a pointer to two structures defined in
// klu.h:  trilinos_klu_symbolic and trilinos_klu_numeric.  This prevents Amesos_Klu.h
// from having to include klu.h.
//
//  Doxygen does not handle forward class references well.
#ifndef DOXYGEN_SHOULD_SKIP_THIS
class Amesos_Klu_Pimpl ;
class Amesos_StandardIndex ;
#endif
 
class Amesos_Klu: public Amesos_BaseSolver,
                  private Amesos_Time,
                  private Amesos_NoCopiable,
                  private Amesos_Utils,
                  private Amesos_Control,
                  private Amesos_Status {
 
public:


  Amesos_Klu(const Epetra_LinearProblem& LinearProblem );

  ~Amesos_Klu(void);


 
  int SymbolicFactorization() ;
 
  int NumericFactorization() ;
 
  int Solve();



  const Epetra_LinearProblem *GetProblem() const { return(Problem_); };


  bool MatrixShapeOK() const ;


  int SetUseTranspose(bool UseTranspose_in) {UseTranspose_ = UseTranspose_in; return(0);};
 
  bool UseTranspose() const {return(UseTranspose_);};
 
  const Epetra_Comm & Comm() const {return(GetProblem()->GetOperator()->Comm());};
 
  int SetParameters( Teuchos::ParameterList &ParameterList );

  int NumSymbolicFact() const { return( Amesos_Status::NumSymbolicFact_ ); }

  int NumNumericFact() const { return( Amesos_Status::NumNumericFact_ ); }

  int NumSolve() const { return( Amesos_Status::NumSolve_ ); }

  void PrintTiming() const;

  void PrintStatus() const;

  void GetTiming( Teuchos::ParameterList &TimingParameterList ) const { Amesos_Time::GetTiming( TimingParameterList ); }
 
private:


 
  /*
  CreateLocalMatrixAndExporters - Prepare to convert matrix and vectors to serial
    Preconditions:
      Problem_ must be set
 
    Postconditions:
      UseDataInPlace_ is set to 1 if the input matrix can be used in place, i.e.
        1)  is entirely stored on process 0
        2)  range map and domain map are same as the row map
      The following are only set if (! UseDataInPlace_ )"
        SerialMap_
        ImportToSerial_
        SerialCrsMatrixA_
 
      SerialMatrix_
   */
  int CreateLocalMatrixAndExporters() ;
  /*
    ExportToSerial
    Preconditions:
       UseDataInPlace_ must be set
       ImportToSerial and SerialCrsMatrixA_ must be set if UseDataInPlace_ != 1
    Postconditions
       SerialMatrix_ points to a serial version of the matrix
   */
  int ExportToSerial() ;
  /*
    ConvertToKluCRS - Convert matrix to form expected by Klu: Ai, Ap, Aval
    Preconditions:
      numentries_, RowMatrixA_, ImportToSerial_, StdIndexMatrix_, Reindex_
    Postconditions:
      SerialCrsMatrixA_
  */
  int ConvertToKluCRS(bool firsttime);
 
  /*
    PerformSymbolicFactorization - Call Klu to perform symbolic factorization
    Preconditions:
      UseDataInPlace_ must be set to 1 if the input matrix is entirely stored on process 0
      Ap, Ai and Aval point to a compressed row storage version of the input matrix A.
    Postconditions:
      Symbolic points to an KLU internal opaque object containing the
        symbolic factorization and accompanying information.
      SymbolicFactorizationOK_ = true;
    Note:  All action is performed on process 0
 
    Returns non-zero if the symbolic factorization failed
  */
 
  int PerformSymbolicFactorization();
 
  /*
    PerformNumericFactorization - Call Klu to perform numeric factorization
    Preconditions:
      UseDataInPlace_ must be set
      Ap, Ai and Aval point to a compressed row storage version of the input matrix A.
      Symbolic must be set
    Postconditions:
      Numeric points to an KLU internal opaque object containing the
        numeric factorization and accompanying information.
      NumericFactorizationOK_ = true;
    Note:  All action is performed on process 0
  */
  int PerformNumericFactorization();
 
  // @}
 
  int SerialXlda_ ;
 
#ifdef Bug_8212
  int *lose_this_;
#endif
  //
  //  PrivateKluData_ contains pointers to data needed by klu whose
  //  data structures are defined by klu.h
  //
  Teuchos::RCP<Amesos_Klu_Pimpl> PrivateKluData_;
  Teuchos::RCP<Amesos_StandardIndex> StdIndex_;
  Teuchos::RCP<Amesos_StandardIndex> StdIndexRange_;
  Teuchos::RCP<Amesos_StandardIndex> StdIndexDomain_;

  std::vector <int> Ap;
  std::vector <int> VecAi;
  std::vector <double> VecAval;
  double* Aval;
  int *Ai;

  int UseDataInPlace_;
  long long numentries_;
  long long NumGlobalElements_;

  Epetra_RowMatrix* RowMatrixA_;
  Epetra_CrsMatrix* CrsMatrixA_;
#if 0
  Teuchos::RCP<EpetraExt::MultiVector_Reindex> VecTrans_;
  Teuchos::RCP<EpetraExt::CrsMatrix_Reindex> MatTrans_;
  Teuchos::RCP<Epetra_Map> ContiguousMap_;
#endif
  Teuchos::RCP<Epetra_Map> SerialMap_;
  Teuchos::RCP<Epetra_CrsMatrix> SerialCrsMatrixA_;
  Epetra_RowMatrix* StdIndexMatrix_ ;
  Epetra_RowMatrix* SerialMatrix_ ;

  //  serially, StorageOptimized, the LHS and RHS are assumed to be available
  //  when SymbolicFactorization is called and not to change (address or number
  //  of vectors) thereafter.
  bool TrustMe_;
  int NumVectors_;
  double *SerialXBvalues_ ;
  double *SerialBvalues_ ;
  Teuchos::RCP<Epetra_MultiVector> SerialB_ ;
  Teuchos::RCP<Epetra_MultiVector> SerialX_ ;
  Teuchos::RCP<Epetra_MultiVector> SerialXextract_;
  Teuchos::RCP<Epetra_MultiVector> SerialBextract_;

  bool UseTranspose_;
  const Epetra_LinearProblem * Problem_;

  std::vector<int> ColIndicesV_;
  std::vector<double> RowValuesV_;
  Teuchos::RCP<Epetra_Import> ImportToSerial_;
  Teuchos::RCP<Epetra_Import> ImportRangeToSerial_;
  Teuchos::RCP<Epetra_Import> ImportDomainToSerial_;
  int MtxRedistTime_, MtxConvTime_, VecRedistTime_;
  int SymFactTime_, NumFactTime_, SolveTime_, OverheadTime_;
 
};  // class Amesos_Klu
 
#endif /* AMESOS_KLU_H */