src/eigsolv/HierarchicalBasisPrec.cpp

Go to the documentation of this file.

#include "rlog/rlog.h"

#include "Amesos_config.h"
#include "Amesos.h"
#include "Teuchos_ParameterList.hpp"

#include "Epetra_LinearProblem.h"
#include "Epetra_CrsMatrix.h"
#include "Epetra_MultiVector.h"
#include "Epetra_Map.h"
#include "Epetra_Comm.h"
#include "AztecOO.h"
#include "AztecOO_Operator.h"

#ifdef HAVE_AMESOS_SUPERLUDIST
//commented out by Tizo
//because of double definitions
//in Tuechos_LAPACK_wrappers.h
//#include <superlu_ddefs.h>
//#include <supermatrix.h>
#endif

#include "OutputCapturer.h"
#include "pbedefs.h"
#include "myrlog.h"

#include "BlockPCGSolver.h"
#include "identityoperator.h"
#include "diagonalpreconditioner.h"
#include "BlockDiagOperator.h"

#include "Ifpack.h"
#include "Ifpack_IC.h"

#include "HierarchicalBasisPrec.h"

HierarchicalBasisPrec::HierarchicalBasisPrec(Epetra_CrsMatrix* Block11,
                                             Epetra_CrsMatrix* TMatrix,
                                             Epetra_CrsMatrix* NodeMatrix,
                                             Epetra_CrsMatrix* Block12,
                                             Epetra_CrsMatrix* Block22,
                                             const Epetra_Map& DMap,
                                             const Epetra_Map& RMap,
                                             Teuchos::ParameterList& params)
    : K11(Block11),
      TMatrix11(TMatrix),
      NodeMatrix11(NodeMatrix),
      K12(Block12),
      K22(Block22),
      DomainMap(DMap),
      RangeMap(RMap),
      Abase(0),
      ProblemL(0),
      solver22_(0),
      operator22_(0),
      operator11_(0),
      _Comm(Block11->Comm()),
      _useTranspose(false),
      isFactorized(false),
      params_(params)
{
    pbe_start(110, "2-level preconditioner: construction");

    //**************************************************************************
    // Setup solver for (1,1)-block

    if (params_.get<string>("solver11") == "lu") {
        ProblemL = new Epetra_LinearProblem();
        ProblemL->SetOperator( K11 );
        Amesos Afactory;
        Abase = Afactory.Create( "Amesos_Superludist", *ProblemL );
        assert(Abase != 0);
        int ierr = Abase->SymbolicFactorization();
        assert(ierr == 0);
        ierr = Abase->NumericFactorization();
        assert(ierr == 0);
        isFactorized = true;

        /*
          ParamList = new Teuchos::ParameterList();
          ParamList->set( "Redistribute", false );
          ParamList->set( "AddZeroToDiag", false );
          Teuchos::ParameterList& SuperludistParams = ParamList->sublist("Superludist") ;
          SuperludistParams.set( "ReuseSymbolic", true );
          SuperludistParams.set( "MaxProcesses", 2 );
          Abase->SetParameters(*ParamList);
        */
    } else if (params_.get<string>("solver11") == "pcg_ml") {
        // pcg_ml is only suitable for H
        rAssert(TMatrix11 == 0 && NodeMatrix11 == 0);
        
        const int maxIterCG = 500;
        const int amg_H_smoother_degree = 2;
        operator11_ = new BlockPCGSolver(_Comm, K11, 1e-10, maxIterCG, 0);
        // Set an AMG preconditioner for PCG
        // param is the degree of polynomial for the smoother
        dynamic_cast<BlockPCGSolver*>(operator11_)->setAMGPreconditioner(1, amg_H_smoother_degree);
    } else if (params_.get<string>("solver11") == "ml") {
        OutputCapturer cap;
        cap.begin_capture();
        if ( TMatrix11 == 0 || NodeMatrix11 == 0 ) {
            operator11_ = new ML_Epetra::MultiLevelPreconditioner(*K11, params_.sublist("ml_params"), true);
        }
        else {
            operator11_ = new ML_Epetra::MultiLevelPreconditioner(*K11,
                                                                  *TMatrix11,
                                                                  *NodeMatrix11,
                                                                  params_.sublist("ml_params"),
                                                                  true);
        }
        cap.end_capture();
        cap.log("ML output:", RLOG_CHANNEL("debug"));

        ostringstream buf;
        buf << "Unused ML parameters:" << endl;
        dynamic_cast<ML_Epetra::MultiLevelPreconditioner*>(operator11_)->PrintUnused(buf);
        rDebug(buf.str().c_str());
    } else if (params_.get<string>("solver11") == "if") {
        // Approximately solve (1,1)-block by a singe application of a incomplete 
        // LDLT-factorisation.
    
        // allocates an IFPACK factory. No data is associated
        // to this object (only method Create()).
        Ifpack Factory;
    
        // Parameters for IFPACK
        Teuchos::ParameterList List;

        // Create the preconditioner. For valid PrecType values,
        // please check the documentation
#if 0
        // FIXME: ICT is not a LDL^T method and cannot be used for indefinite matrices
        string PrecType = "ICT";
        List.set("fact: ict level-of-fill", 2.0);
#endif        

#if 0
        // FIXME: ICT is not a LDL^T method and cannot be used for indefinite matrices
        string PrecType = "ICT stand-alone";
        List.set("fact: ict level-of-fill", 2.0);
#endif        

#if 1
        // LDL^T
        string PrecType = "IC stand-alone";
        List.set("fact: level-of-fill", 1);
#endif        
        
#if 0
        string PrecType = "ILU stand-alone";
        List.set("fact: level-of-fill", 2);
#endif        
        
        int OverlapLevel = 1; // must be >= 0. If Comm.NumProc() == 1,
                              // it is ignored.
        operator11_ = Factory.Create(PrecType, K11, OverlapLevel);
        assert(operator11_ != 0);
    
        dynamic_cast<Ifpack_Preconditioner*>(operator11_)->SetParameters(List);
    
        // Initialize the preconditioner. At this point the matrix must
        // have been FillComplete()'d, but actual values are ignored.
        dynamic_cast<Ifpack_Preconditioner*>(operator11_)->Initialize();
    
        // Builds the preconditioners, by looking for the values of
        // the matrix.
        dynamic_cast<Ifpack_Preconditioner*>(operator11_)->Compute();
        
        // Information about the preconditioner
        cout << *dynamic_cast<Ifpack_Preconditioner*>(operator11_);
        Ifpack_IC* ic_prec = dynamic_cast<Ifpack_IC*>(operator11_);
        if (ic_prec != 0)
            rInfo("#nonzeros of IF-preconditioner for K11: %d", ic_prec->NumGlobalNonzeros());
            
    } else {
        rExit("Unkown solver for (1,1)-block of hierarchical preconditioner: %s", 
              params_.get<string>("solver11").c_str());
    }

    //**************************************************************************
    // Setup solver for (2,2)-block

    if (params_.get<string>("solver22") == "gmres_incomplete_block_jacobi") {
        // CGS with block-Jacobi (with incomplete factorisation of the diagonal blocks)
        // preconditioner.
        solver22_ = new AztecOO();
        solver22_->SetUserMatrix(K22);
        solver22_->SetAztecOption(AZ_solver, AZ_cgs);
        solver22_->SetAztecOption(AZ_precond, AZ_dom_decomp);
        solver22_->SetAztecOption(AZ_subdomain_solve, AZ_ilu);
        solver22_->SetAztecOption(AZ_output, AZ_none);
        operator22_ = new AztecOO_Operator(solver22_, 2);
    } else if (params_.get<string>("solver22") == "pcg_jacobi") {
        // PCG with Point-Jacobi precondioner
        solver22_ = new AztecOO();
        solver22_->SetUserMatrix(K22);
        solver22_->SetAztecOption(AZ_solver, AZ_cg);
        solver22_->SetAztecOption(AZ_precond, AZ_Jacobi);
        solver22_->SetAztecOption(AZ_poly_ord, 1);
        solver22_->SetAztecOption(AZ_output, AZ_none);
        operator22_ = new AztecOO_Operator(solver22_, 2);
    } else if (params_.get<string>("solver22") == "diag") {
        // Diagonal preconditioner
        operator22_ = new DiagonalPreconditioner(*K22);
    } else if (params_.get<string>("solver22") == "blockdiag") {
        // Block Diagonal preconditioner
        operator22_ = new BlockDiagOperator(K22);
    } else if (params_.get<string>("solver22") == "none") {
        operator22_ = new IdentityOperator(K22->OperatorDomainMap(),
                                           K22->OperatorRangeMap(),
                                           K22->Comm());
    } else {
        rExit("Unkown solver for (2,2)-block of hierarchical preconditioner: %s",
              params_.get<string>("solver22").c_str());
    }
    pbe_stop(110);
}

HierarchicalBasisPrec::~HierarchicalBasisPrec() {
    if (Abase) 
        delete Abase;
    if (operator11_) 
        delete operator11_;
    if (solver22_)
        delete solver22_;
    if (operator22_)
        delete operator22_;
}

int HierarchicalBasisPrec::SetUseTranspose(bool UseTranspose) {
    int ierr = 0;
    _useTranspose = UseTranspose;
    return(ierr);
}

int HierarchicalBasisPrec::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    int ierr = 0;

    pbe_start(111, "2-level preconditioner: operator");

    if (!X.Map().SameAs(OperatorDomainMap()))
        EPETRA_CHK_ERR(-1);
    if (!Y.Map().SameAs(OperatorRangeMap()))
        EPETRA_CHK_ERR(-2);
    if (Y.NumVectors() != X.NumVectors())
        EPETRA_CHK_ERR(-3);

    Epetra_MultiVector X1(View, K11->DomainMap(), X.Values(), X.MyLength(), X.NumVectors());
    Epetra_MultiVector X2(View, K22->DomainMap(), X.Values() + X1.MyLength(), X.MyLength(), X.NumVectors());
    Epetra_MultiVector Y1(View, K11->RangeMap(), Y.Values(), Y.MyLength(), Y.NumVectors());
    Epetra_MultiVector Y2(View, K22->RangeMap(), Y.Values() + Y1.MyLength(), Y.MyLength(), Y.NumVectors());

    double *tt = new double[Y1.NumVectors()*Y1.MyLength()];
    Epetra_MultiVector Ytmp1(View, K11->RangeMap(), tt, Y1.MyLength(), Y1.NumVectors());
    K11->Multiply(false, X1, Y1);
    K12->Multiply(false, X2, Ytmp1);
    Y1.Update(1.0, Ytmp1, 1.0);
    delete[] tt;

    tt = new double[Y2.NumVectors()*Y2.MyLength()];
    Epetra_MultiVector Ytmp2(View, K22->RangeMap(), tt, Y2.MyLength(), Y2.NumVectors());
    K12->Multiply(true, X1, Ytmp2);
    K22->Multiply(false, X2, Y2);
    Y2.Update(1.0, Ytmp2, 1.0);
    delete[] tt;
    pbe_stop(111);

    return(ierr);
}

int HierarchicalBasisPrec::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {
    int ierr = 0;

    //We suppose :
    //K11->FillComplete(map1, map1);
    //K12->FillComplete(map2, map1);
    //K22->FillComplete(map2, map2);
    //map of X and Y is the ordered composition of (map1,map2)

    if (!X.Map().SameAs(OperatorDomainMap()))
        EPETRA_CHK_ERR(-1);
    if (!Y.Map().SameAs(OperatorRangeMap()))
        EPETRA_CHK_ERR(-2);
    if (Y.NumVectors() != X.NumVectors())
        EPETRA_CHK_ERR(-3);

    pbe_start(112, "2-level preconditioner: application");

    Epetra_MultiVector X1(View, K11->RangeMap(), X.Values(), X.MyLength(), X.NumVectors());
    Epetra_MultiVector X2(View, K22->RangeMap(), X.Values() + X1.MyLength(), X.MyLength(), X.NumVectors());
    Epetra_MultiVector Y1(View, K11->DomainMap(), Y.Values(), Y.MyLength(), Y.NumVectors());
    Epetra_MultiVector Y2(View, K22->DomainMap(), Y.Values() + Y1.MyLength(), Y.MyLength(), Y.NumVectors());
    Epetra_MultiVector Xtmp1(K11->RangeMap(), X.NumVectors(), false);
    Epetra_MultiVector Xtmp2(K22->RangeMap(), X.NumVectors(), false);

    //  This is necessary only in the case X and Y points to the same address in memory.
    Epetra_MultiVector* X11;
    if (X.Values() == Y.Values()) {
        X11 = new Epetra_MultiVector(X1);
    }
    else {
        X11 = &X1;
    }

    //
    //  Solve K11*Y1 = X1
    //
    pbe_start(113, "2-level preconditioner: (1,1)-solve");
    if (params_.get<string>("solver11") == "lu") {
        assert(Abase != 0);
        assert(isFactorized);
        ProblemL->SetLHS( &Y1 );
        ProblemL->SetRHS( &X1 );
        EPETRA_CHK_ERR( Abase->Solve(  ) );
    } else if (params_.get<string>("solver11") == "pcg_ml")  {
        operator11_->ApplyInverse(X1, Y1);
    } else if (params_.get<string>("solver11") == "ml")  {
        operator11_->ApplyInverse(X1, Y1);
    } else if (params_.get<string>("solver11") == "if")  {
        operator11_->ApplyInverse(X1, Y1);
    } else {
        assert(false);
    }
    pbe_stop(113);

    //
    // Solve K22*Y2 = X2 - K12^T * Y1
    //
    pbe_start(115, "2-level preconditioner: off-diagonal multiply");
    K12->Multiply(true, Y1, Xtmp2);
    pbe_stop(115);
    Xtmp2.Update(1.0, X2, -1.0);
    pbe_start(114, "2-level preconditioner: (2,2)-solve");
    operator22_->ApplyInverse(Xtmp2, Y2);
    pbe_stop(114);

    //
    // Solve K11*Y1 = X1 - K12*Y2
    //
    pbe_start(115, "2-level preconditioner: (1,2)/(2,1)-mult");
    K12->Multiply(false, Y2, Xtmp1);
    pbe_stop(115);
    Xtmp1.Update(1.0, X1, -1.0);
    pbe_start(113, "2-level preconditioner: (1,1)-solve");
    if (params_.get<string>("solver11") == "lu") {
        ProblemL->SetLHS( &Y1 );
        ProblemL->SetRHS( &Xtmp1 );
        EPETRA_CHK_ERR( Abase->Solve(  ) );
    } else if (params_.get<string>("solver11") == "pcg_ml")  {
        operator11_->ApplyInverse(Xtmp1, Y1);
    } else if (params_.get<string>("solver11") == "ml")  {
        operator11_->ApplyInverse(Xtmp1, Y1);
    } else if (params_.get<string>("solver11") == "if")  {
        operator11_->ApplyInverse(X1, Y1);
    } else {
        assert(false);
    }

    if (X11 != &X1)
        delete X11;

    pbe_stop(113);
    pbe_stop(112);

    return(ierr);
}

double HierarchicalBasisPrec::NormInf() const {
    cout << "Not implemented yet !!!!!!!!!!" << endl;
    return(1);
}

const char* HierarchicalBasisPrec::Label() const {
    return "Hierarchical preconditioner";
}

bool HierarchicalBasisPrec::UseTranspose() const {
    return(false);
}

bool HierarchicalBasisPrec::HasNormInf() const {
    return(false);
}

const Epetra_Comm& HierarchicalBasisPrec::Comm() const {
    return(_Comm);
}

const Epetra_Map& HierarchicalBasisPrec::OperatorDomainMap() const {
    return(DomainMap);
}

const Epetra_Map& HierarchicalBasisPrec::OperatorRangeMap() const {
    return(RangeMap);
}

---