Ifpack_IC.cpp

/*@HEADER
// ***********************************************************************
// 
//       Ifpack: Object-Oriented Algebraic Preconditioner Package
//                 Copyright (2002) 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
// 
// ***********************************************************************
//@HEADER
*/

#include "Ifpack_ConfigDefs.h"
#ifdef HAVE_IFPACK_TEUCHOS
#include "Ifpack_Preconditioner.h"
#include "Ifpack_IC.h"
#include "Ifpack_IC_Utils.h"
#include "Ifpack_Condest.h"
#include "Epetra_Comm.h"
#include "Epetra_Map.h"
#include "Epetra_RowMatrix.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_Vector.h"
#include "Epetra_MultiVector.h"
#include "Epetra_Util.h"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_RefCountPtr.hpp"
using namespace Teuchos;

//==============================================================================
Ifpack_IC::Ifpack_IC(Epetra_RowMatrix* A) :
  A_(*A),
  Comm_(A->Comm()),
  UseTranspose_(false),
  Condest_(-1.0),
  Athresh_(0.0),
  Rthresh_(1.0),
  Droptol_(0.0),
  Lfil_(0),
  Aict_(0),
  Lict_(0),
  Ldiag_(0),
  IsInitialized_(false),
  IsComputed_(false),
  NumInitialize_(0),
  NumCompute_(0),
  NumApplyInverse_(0),
  InitializeTime_(0.0),
  ComputeTime_(0),
  ApplyInverseTime_(0),
  ComputeFlops_(0.0),
  ApplyInverseFlops_(0.0)
{
#ifdef HAVE_IFPACK_TEUCHOS
  Teuchos::ParameterList List;
  SetParameters(List);
#endif

}
//==============================================================================
Ifpack_IC::~Ifpack_IC()
{
  if (Lict_ != 0) {
    Ifpack_AIJMatrix * Lict = (Ifpack_AIJMatrix *) Lict_;
    delete [] Lict->ptr;
    delete [] Lict->col;
    delete [] Lict->val;
    delete Lict;
  }
  if (Aict_ != 0) {
    Ifpack_AIJMatrix * Aict = (Ifpack_AIJMatrix *) Aict_;
    delete Aict;
  }
  if (Ldiag_ != 0) delete [] Ldiag_;

  IsInitialized_ = false;
  IsComputed_ = false;
}

//==========================================================================
int Ifpack_IC::SetParameters(Teuchos::ParameterList& List)
{

  Lfil_ = List.get("fact: level-of-fill",Lfil_);
  Athresh_ = List.get("fact: absolute threshold", Athresh_);
  Rthresh_ = List.get("fact: relative threshold", Rthresh_);
  Droptol_ = List.get("fact: drop tolerance", Droptol_);

  // set label
  sprintf(Label_, "IFPACK IC (fill=%d, drop=%f)",
      Lfil_, Droptol_);
  return(0);
}

//==========================================================================
int Ifpack_IC::Initialize()
{

  IsInitialized_ = false;
  // FIXME: construction of ILUK graph must be here
  
  IsInitialized_ = true;
  return(0);
}

//==========================================================================
int Ifpack_IC::ComputeSetup()
{
  // (re)allocate memory for ICT factors
  U_ = rcp(new Epetra_CrsMatrix(Copy, Matrix().RowMatrixRowMap(), 
                                Matrix().RowMatrixRowMap(), 0));
  D_ = rcp(new Epetra_Vector(Matrix().RowMatrixRowMap()));

  if (U_.get() == 0 || D_.get() == 0)
    IFPACK_CHK_ERR(-5); // memory allocation error

  int ierr = 0;
  int i, j;
  int* InI;
  int* UI;
  double* InV;
  double* UV;
  int NumIn, NumL, NumU;
  bool DiagFound;
  int NumNonzeroDiags = 0;

  // Get Maximun Row length
  int MaxNumEntries = Matrix().MaxNumEntries();

  InI = new int[MaxNumEntries]; // Allocate temp space
  UI  = new int[MaxNumEntries];
  InV = new double[MaxNumEntries];
  UV  = new double[MaxNumEntries];

  double *DV;
  ierr = D_->ExtractView(&DV); // Get view of diagonal

  // First we copy the user's matrix into diagonal vector and U, regardless of fill level

  int NumRows = Matrix().NumMyRows();

  for (i=0; i< NumRows; i++) {

    Matrix().ExtractMyRowCopy(i, MaxNumEntries, NumIn, InV, InI); // Get Values and Indices
    
    // Split into L and U (we don't assume that indices are ordered).
    NumL = 0; 
    NumU = 0; 
    DiagFound = false;
    
    for (j=0; j< NumIn; j++) {
      int k = InI[j];

      if (k==i) {
    DiagFound = true;
    // Store perturbed diagonal in Epetra_Vector D_
    DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; 
      }

      else if (k < 0) return(-1); // Out of range
      else if (i<k && k<NumRows) {
    UI[NumU] = k;
    UV[NumU] = InV[j];
    NumU++;
      }
    }
    
    // Check in things for this row of L and U

    if (DiagFound) NumNonzeroDiags++;
    if (NumU) U_->InsertMyValues(i, NumU, UV, UI);
    
  }

  delete [] UI;
  delete [] UV;
  delete [] InI;
  delete [] InV;

  U_->FillComplete(Matrix().OperatorDomainMap(), 
               Matrix().OperatorRangeMap());

  int ierr1 = 0;
  if (NumNonzeroDiags<U_->NumMyRows()) ierr1 = 1;
  Matrix().Comm().MaxAll(&ierr1, &ierr, 1);
  IFPACK_CHK_ERR(ierr);

  return(0);
}

//==========================================================================
int Ifpack_IC::Compute() {

  if (!IsInitialized()) 
    IFPACK_CHK_ERR(Initialize());

  IsComputed_ = false;

  // copy matrix into L and U.
  IFPACK_CHK_ERR(ComputeSetup());

  int i;

  int m, n, nz, Nrhs, ldrhs, ldlhs;
  int * ptr=0, * ind;
  double * val, * rhs, * lhs;

  int ierr = Epetra_Util_ExtractHbData(U_.get(), 0, 0, m, n, nz, ptr, ind,
                       val, Nrhs, rhs, ldrhs, lhs, ldlhs);
  if (ierr < 0) 
    IFPACK_CHK_ERR(ierr);

  Ifpack_AIJMatrix * Aict;
  if (Aict_==0) {
    Aict = new Ifpack_AIJMatrix;
    Aict_ = (void *) Aict;
  }
  else Aict = (Ifpack_AIJMatrix *) Aict_;
  Ifpack_AIJMatrix * Lict;
  if (Lict_==0) {
    Lict = new Ifpack_AIJMatrix;
    Lict_ = (void *) Lict;
  }
  else Lict = (Ifpack_AIJMatrix *) Lict_;
  Aict->val = val;
  Aict->col = ind;
  Aict->ptr = ptr;
  double *DV;
  EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal
    
  crout_ict(m, Aict, DV, Droptol_, Lfil_, Lict, &Ldiag_);

  // Get rid of unnecessary data
  delete [] ptr;

  // Create Epetra View of L from crout_ict
  U_ = rcp(new Epetra_CrsMatrix(View, A_.RowMatrixRowMap(), A_.RowMatrixRowMap(),0));
  D_ = rcp(new Epetra_Vector(View, A_.RowMatrixRowMap(), Ldiag_));

  ptr = Lict->ptr;
  ind = Lict->col;
  val = Lict->val;
    
  for (i=0; i< m; i++) {
    int NumEntries = ptr[i+1]-ptr[i];
    int * Indices = ind+ptr[i];
    double * Values = val+ptr[i];
    U_->InsertMyValues(i, NumEntries, Values, Indices);
  }

  U_->FillComplete(A_.OperatorDomainMap(), A_.OperatorRangeMap());
  D_->Reciprocal(*D_); // Put reciprocal of diagonal in this vector
 
  double current_flops = 2 * nz; // Just an estimate
  double total_flops = 0;
    
  A_.Comm().SumAll(&current_flops, &total_flops, 1); // Get total madds across all PEs

  ComputeFlops_ += total_flops; 
  // Now count the rest. NOTE: those flops are *global*
  ComputeFlops_ += (double) U_->NumGlobalNonzeros(); // Accounts for multiplier above
  ComputeFlops_ += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal

  IsComputed_ = true;

  return(0);

}

//=============================================================================
// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
int Ifpack_IC::ApplyInverse(const Epetra_MultiVector& X, 
                 Epetra_MultiVector& Y) const
{

  if (!IsComputed())
    IFPACK_CHK_ERR(-3); // compute preconditioner first

  if (X.NumVectors() != Y.NumVectors()) 
    IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size

  bool Upper = true;
  bool UnitDiagonal = true;

  // AztecOO gives X and Y pointing to the same memory location,
  // need to create an auxiliary vector, Xcopy
  const Epetra_MultiVector* Xcopy;
  if (X.Pointers()[0] == Y.Pointers()[0])
    Xcopy = new Epetra_MultiVector(X);
  else
    Xcopy = &X;

  U_->Solve(Upper, true, UnitDiagonal, *Xcopy, Y);
  Y.Multiply(1.0, *D_, Y, 0.0); // y = D*y (D_ has inverse of diagonal)
  U_->Solve(Upper, false, UnitDiagonal, Y, Y); // Solve Uy = y

  if (Xcopy != &X)
    delete Xcopy;

  ++NumApplyInverse_;
  ApplyInverseFlops_ += 4.0 * U_->NumGlobalNonzeros();
  ApplyInverseFlops_ += D_->GlobalLength();
  return(0);

}

//=============================================================================
// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
int Ifpack_IC::Apply(const Epetra_MultiVector& X, 
              Epetra_MultiVector& Y) const 
{

  if (X.NumVectors() != Y.NumVectors()) 
    IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size

  Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;
  Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;

  U_->Multiply(false, *X1, *Y1);
  Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)
  Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)
  Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1
  U_->Multiply(true, Y1temp, *Y1);
  Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)
  return(0);
}

//=============================================================================
double Ifpack_IC::Condest(const Ifpack_CondestType CT, 
                            const int MaxIters, const double Tol,
                            Epetra_RowMatrix* Matrix)
{
  if (!IsComputed()) // cannot compute right now
    return(-1.0);

  if (Condest_ == -1.0)
    Condest_ = Ifpack_Condest(*this, CT, MaxIters, Tol, Matrix);

  return(Condest_);
}

//=============================================================================
std::ostream&
Ifpack_IC::Print(std::ostream& os) const
{
  if (!Comm().MyPID()) {
    os << endl;
    os << "================================================================================" << endl;
    os << "Ifpack_IC: " << Label() << endl << endl;
    os << "Level-of-fill      = " << LevelOfFill() << endl;
    os << "Absolute threshold = " << AbsoluteThreshold() << endl;
    os << "Relative threshold = " << RelativeThreshold() << endl;
    os << "Drop tolerance     = " << DropTolerance() << endl;
    os << "Condition number estimate = " << Condest() << endl;
    os << "Global number of rows            = " << A_.NumGlobalRows() << endl;
    if (IsComputed_) {
      os << "Number of nonzeros of H         = " << U_->NumGlobalNonzeros() << endl;
      os << "nonzeros / rows                 = " 
         << 1.0 * U_->NumGlobalNonzeros() / U_->NumGlobalRows() << endl;
    }
    os << endl;
    os << "Phase           # calls   Total Time (s)       Total MFlops     MFlops/s" << endl;
    os << "-----           -------   --------------       ------------     --------" << endl;
    os << "Initialize()    "   << std::setw(5) << NumInitialize() 
       << "  " << std::setw(15) << InitializeTime() 
       << "               0.0            0.0" << endl;
    os << "Compute()       "   << std::setw(5) << NumCompute() 
       << "  " << std::setw(15) << ComputeTime()
       << "  " << std::setw(15) << 1.0e-6 * ComputeFlops();
    if (ComputeTime() != 0.0)
      os << "  " << std::setw(15) << 1.0e-6 * ComputeFlops() / ComputeTime() << endl;
    else
      os << "  " << std::setw(15) << 0.0 << endl;
    os << "ApplyInverse()  "   << std::setw(5) << NumApplyInverse() 
       << "  " << std::setw(15) << ApplyInverseTime()
       << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops();
    if (ApplyInverseTime() != 0.0) 
      os << "  " << std::setw(15) << 1.0e-6 * ApplyInverseFlops() / ApplyInverseTime() << endl;
    else
      os << "  " << std::setw(15) << 0.0 << endl;
    os << "================================================================================" << endl;
    os << endl;
  }

  
  return(os);
} 
#endif // HAVE_IFPACK_TEUCHOS

---