Amesos_Scalapack.h

Go to the documentation of this file.

/*
  TODO:
  1)  Eliminate mb_ and grid_mb_  DONE 

  1)  Allow user level control over whether to use a 1D or 2D data 
      distribution.
  2)  Figure out how to redistribute the vectors.  (Allow no more than 
      nb vectors at a time for now).  Done - I think 
  3)  Move the code from software:MyExamples/TwodMap/TwoDMap.cpp in to
      this routine (for the 2D case) Done - I think
  4)  Create the ScaLAPACK 2D dense matrix.  Done - I think 
 */

// @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_SCALAPACK_H_
#define _AMESOS_SCALAPACK_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"
#ifdef EPETRA_MPI
#include "Epetra_MpiComm.h"
#else
#include "Epetra_Comm.h"
#endif
#include "Epetra_CrsGraph.h"


class Amesos_Scalapack: public Amesos_BaseSolver,
  private Amesos_Time,
  private Amesos_NoCopiable,
  private Amesos_Utils,
  private Amesos_Control, 
  private Amesos_Status 
 { 

public: 


  Amesos_Scalapack( const Epetra_LinearProblem& LinearProblem );


  ~Amesos_Scalapack(void);



    int SymbolicFactorization() ;


    int NumericFactorization() ;


    int Solve();

  

#if 0
  char * Label() const {return(Epetra_Object::Label());};
#endif
    
  const Epetra_LinearProblem *GetProblem() const { return(Problem_); };


  bool MatrixShapeOK() const ;


  int SetUseTranspose(bool UseTranspose) {UseTranspose_ = UseTranspose; 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:  
  /*
  RedistributeA - Convert matrix to a dense ScaLAPACK matrix
    Preconditions:
      Problem_ must be set 
      SetParameters() 
      ScaLAPACK1DMap and ScaLAPACK1DMatrix must either be 0 or be pointers to 
        appropriatly allocate objects.  If they are non-zero, those objects
  will be deleted and recreated.  
  
    Postconditions:
      nprow_, npcol_, DescA_ 
      RowMatrixA_ 
      ScaLAPACK1DMap_ 
      ScaLAPACK1DMatrix_ 
      ImportToScaLAPACK1D_
      ImportBackToOriginal_

   */
  int RedistributeA();

  /*
    ConvertToScalapack - Convert matirx to form expected by Scalapack: Ai, Ap, Aval
    Preconditions:
      Problem_ 
    Postconditions:
      nprow_, npcol_, 
  */
  int ConvertToScalapack();     

  /*
    PerformNumericFactorization - Call Scalapack to perform numeric factorization
    Preconditions:

    Postconditions:
      DenseA_, DescA_
  */
  int PerformNumericFactorization(); 

 protected:

  int iam_;                              //  Process number (i.e. Comm().MyPID() 
  
  int NumGlobalElements_;                //  Number of rows and columns in the Problem_->GetOperator()
  int NumGlobalNonzeros_; 

  //
  //  The following variables are required for the ScaLAPACK interface:
  //
  int nprow_ ;                           //  number of process rows: 1 for now
  int npcol_ ;                           //  number of process columns
  int ictxt_ ;                           //  BLACS context
  int m_per_p_;                          //  Number of columns per process
  int DescA_[10];                        //  ScaLAPACK array descriptor 

  Epetra_Map *ScaLAPACK1DMap_ ;          //  Points to a 1D Map which matches a ScaLAPACK 1D
                                         //  blocked (not block cyclic) distribution
  Epetra_CrsMatrix *ScaLAPACK1DMatrix_ ; //  Points to a  ScaLAPACK 1D
                                         //  blocked (not block cyclic) distribution
  Epetra_Map *VectorMap_ ;               //  Points to a Map for vectors X and B
  std::vector<double> DenseA_;                //  The data in a ScaLAPACK 1D blocked format
  std::vector<int> Ipiv_ ;                    //  ScaLAPACK pivot information
  int NumOurRows_ ;
  int NumOurColumns_ ;
 

  bool UseTranspose_;     
  const Epetra_LinearProblem * Problem_;
  
  
  // some timing internal to ScaLAPACK
  double ConTime_;                        // time to convert to ScaLAPACKformat
  double SymTime_;                        // time for symbolic factorization
  double NumTime_;                        // time for numeric factorization
  double SolTime_;                        // time for solution
  double VecTime_;                        // time to redistribute vectors
  double MatTime_;                        // time to redistribute matrix
  
  //
  //  Control of the data distribution
  //
  bool TwoD_distribution_;  // True if 2D data distribution is used
  int grid_nb_;             // Row and Column blocking factor (only used in 2D distribution)  
  int mypcol_;              // Process column in the ScaLAPACK2D grid
  int myprow_;              // Process row in the ScaLAPACK2D grid
  Epetra_CrsMatrix* FatOut_;//

  //
  //  Blocking factors (For both 1D and 2D data distributions)
  //
  int nb_;
  int lda_;

  Epetra_Time * Time_;
  
};  // End of  class Amesos_Scalapack  
#endif /* _AMESOS_SCALAPACK_H_ */