src/eigsolv/LinearEigsolvOperators.cpp

Go to the documentation of this file.

#include "myrlog.h"

#include "ml_include.h"
#include "EpetraExt_MatrixMatrix.h"
#include "matrix_matrix.h"

// for the ML Maxwell preconditioner...
#if !defined(HAVE_ML_EPETRA) || !defined(HAVE_ML_TEUCHOS)
#error Trilinos needs to be configured with --enable-epetra --enable-teuchos
#endif
#include "ml_epetra_preconditioner.h"
#include "Teuchos_ParameterList.hpp"
#include "Ifpack.h"
#include "Ifpack_IC.h"

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

#include "directsolveroperator.h"
#include "BlockPCGSolver.h"
#include "diagonalpreconditioner.h"
#include "LinearEigsolvOperators.h"

using namespace Teuchos;

LinearEigsolvOperators::LinearEigsolvOperators(FemaxMesh& mesh_mixed,
                                               Teuchos::ParameterList& params,
                                               const Epetra_Comm& MyComm) :
    mesh_mixed_(mesh_mixed),
    Comm(MyComm),
    A(0),
    M(0),
    K(0),
    H(0),
    Y(0),
    C(0),
    Prec(0),
    h_solver_op_(0),
    HPrec(0),
    _sigma_K(0.0),
    params_(params),
    Aprec(0),
    AprecSolver(0)
{
    matrixAssembly(mesh_mixed);
}

LinearEigsolvOperators::~LinearEigsolvOperators() {
    OutputCapturer capt;
    capt.begin_capture();
    
    delete Prec;
    delete h_solver_op_;
    delete HPrec;
    delete Aprec;
    delete AprecSolver;
    delete A;
    delete M;
    delete K;
    delete Y;
    delete H;
    delete C;
        
    capt.end_capture();
    // Logging empty strings causes a "[FMT ERROR]: failed to format msg" message in rlog
    if (capt.get_stdout().size() > 0)
        rDebug("%s", capt.get_stdout().c_str());
}

int LinearEigsolvOperators::buildPreconditioner(double sigma) {
    // construct K = A - sigma*M
    K = build_Asigma(*A, *M, sigma);

    if (params_.get<bool>("export_mtx")) {
        pbe_start(130, "MatrixMarket export");
        export_Epetra_CrsMatrix(*K, 15, "K.mtx");
        pbe_stop(130);
    }

    Teuchos::ParameterList& asigma_precon(params_.sublist("asigma_precon"));
    string asigma_precon_type = asigma_precon.get<string>("type");
    if (asigma_precon_type == "neumann") {

        // Setting the AZ_Neumann preconditioner
        Aprec = new Epetra_LinearProblem();
        Aprec->SetOperator(K);

        AprecSolver = new AztecOO(*Aprec);

        // Specify the preconditioner inside AztecOO
        AprecSolver->SetAztecOption(AZ_precond, AZ_Neumann);
        AprecSolver->SetAztecOption(AZ_poly_ord, 10);

        // Specify the output level
        AprecSolver->SetAztecOption(AZ_output, AZ_none);

        // Specify the iterative solver used
        AprecSolver->SetAztecOption(AZ_solver, AZ_gmres);

        // Define an operator for this preconditioner. Preconditioner is then called
        // using Prec->ApplyInverse(...).
        int iterAztec = 2;
        Prec = new AztecOO_Operator(AprecSolver, iterAztec);

    } else if (asigma_precon_type == "ml") {

        OutputCapturer capt;
        capt.begin_capture();
        Prec = new ML_Epetra::MultiLevelPreconditioner(*K, // A - sigma *M
                                                       *Y,
                                                       *H,
                                                       params_.sublist("asigma_precon").sublist("ml_params"),
                                                       false);
        dynamic_cast<ML_Epetra::MultiLevelPreconditioner*>(Prec)->ComputePreconditioner();
        capt.end_capture();
        rLog(RLOG_CHANNEL("debug/all"), "%s", capt.get_stdout().c_str());

        ostringstream buf;
        buf << "Unused ML parameters:" << endl;
        dynamic_cast<ML_Epetra::MultiLevelPreconditioner*>(Prec)->PrintUnused(buf);
        rDebug(buf.str().c_str());

        // The ML preconditioner is applied by Prec->ApplyInverse(...). The Apply() method is not
        // implemented.

    } else if (asigma_precon_type == "lu") {

        // Preconditioner is LU-factorisation of K = A -sigma*M
        Prec = new DirectSolverOperator(*K, params_.sublist("asigma_precon").sublist("amesos_params"));

    } else if (asigma_precon_type == "diag") {

        // Preconditioner is diagonal of K = A -sigma*M
        Prec = new DiagonalPreconditioner(*K);

    } else if (asigma_precon_type == "if") {

        // Preconditioner is incomplete LDLT-factorisation of K = A -sigma*M
    
        // allocates an IFPACK factory. No data is associated
        // to this object (only method Create()).
        Ifpack Factory;
    
        // Create the preconditioner. For valid PrecType values,
        // please check the documentation
        string PrecType = "IC stand-alone"; // incomplete LU
        // string PrecType = "IC"; // incomplete LU
        int OverlapLevel = 1; // must be >= 0. If Comm.NumProc() == 1,
                              // it is ignored.
        Prec = Factory.Create(PrecType, K, OverlapLevel);
        assert(Prec != 0);
    
        // Specify parameters for IC
        Teuchos::ParameterList List;
        List.set("fact: drop tolerance", 1e-2);
        List.set("fact: level-of-fill", 1);
        IFPACK_CHK_ERR(dynamic_cast<Ifpack_Preconditioner*>(Prec)->SetParameters(List));
    
        // Initialize the preconditioner. At this point the matrix must
        // have been FillComplete()'d, but actual values are ignored.
        IFPACK_CHK_ERR(dynamic_cast<Ifpack_Preconditioner*>(Prec)->Initialize());
    
        // Builds the preconditioners, by looking for the values of
        // the matrix.
        IFPACK_CHK_ERR(dynamic_cast<Ifpack_Preconditioner*>(Prec)->Compute());
        
        // Information about the preconditioner
        cout << *dynamic_cast<Ifpack_Preconditioner*>(Prec);
        Ifpack_IC* ic_prec = dynamic_cast<Ifpack_IC*>(Prec);
        if (ic_prec != 0)
            rInfo("#nonzeros of IF-preconditioner for K: %d", ic_prec->NumGlobalNonzeros());
    } else {
        rExit("Unknown preconditioner for A - sigma * M: %s", asigma_precon_type.c_str());
    }

    return(0);
}

void LinearEigsolvOperators::matrixAssembly(FemaxMesh& mesh_mixed) {
    
    NedelecMesh&  mesh_nedelec = mesh_mixed.get_nedelec_mesh();
    LagrangeMesh& mesh_lagrange = mesh_mixed.get_lagrange_mesh();

    // Create Epetra_Maps
#if 1
    std::vector<int> owned_dofs;    // Epetra_Map expects int array
    for (mesh::id_t ldof = 0; ldof < mesh_nedelec.get_num_my_dofs(0); ++ ldof) {
        mesh::id_t gdof = mesh_nedelec.gdof(ldof);
        if (mesh_nedelec.is_owned_gdof(gdof)) {
            mesh::id_t mgdof = mesh_nedelec.map_dof(ldof);
            if (mgdof != mesh::ID_NONE)
                owned_dofs.push_back(mgdof);
        }
    }
    Epetra_Map map_nedelec (-1, owned_dofs.size(), &owned_dofs[0], 0, Comm);
    
    owned_dofs.clear();
    for (mesh::id_t ldof = 0; ldof < mesh_lagrange.get_num_my_dofs(0); ++ ldof) {
        mesh::id_t gdof = mesh_lagrange.gdof(ldof);
        if (mesh_lagrange.is_owned_gdof(gdof)) {
            mesh::id_t mgdof = mesh_lagrange.map_dof(ldof);
            if (mgdof != mesh::ID_NONE)
                owned_dofs.push_back(mgdof);
        }
    }
    Epetra_Map map_lagrange (-1, owned_dofs.size(), &owned_dofs[0], 0, Comm);
#else
    //nof_mapped_elements
    vector<int> ned_dofs = mesh_nedelec.get_my_dofs();
    Epetra_Map map_nedelec (-1, ned_dofs.size(), &ned_dofs[0], 0, Comm);
//     cout << map_nedelec;
    
    //  Create lagrange map in same way as nedelec map.
    vector<int> lag_dofs = mesh_lagrange.get_my_dofs();
    Epetra_Map map_lagrange (-1, lag_dofs.size(), &lag_dofs[0], 0, Comm);
    // cout << map_lagrange;
#endif

    //
    // Matrix assembly
    //
    rInfo("Assemble finite element matrices...");

    // Allocate matrices A and M
    log_mem_footprint(Comm, "before matrix assembly");
    pbe_start(60, "assembly of A and M");
    A = new Epetra_FECrsMatrix(Copy, map_nedelec, false);
    M = new Epetra_FECrsMatrix(Copy, map_nedelec, false);
    mesh_nedelec.assembleAM(A, M);
    
    // Finish data entry
    // transforms to local index space, sorts column indices ans merges redundant entries
    A->FillComplete(map_nedelec, map_nedelec);
    M->FillComplete(map_nedelec, map_nedelec);
    pbe_stop(60);
    log_matrix_stats("Matrix A", A);
    log_matrix_stats("Matrix M", M);
    log_mem_footprint(Comm, "after assembly of A and M");

//     Epetra_FECrsMatrix* A_s = new Epetra_FECrsMatrix(Copy, map_nedelec, false);
//     Epetra_FECrsMatrix* M_s = new Epetra_FECrsMatrix(Copy, map_nedelec, false);
//     mesh_nedelec.generate_sorted_AM(A, M, A_s, M_s);
    
//     A_s->FillComplete(map_nedelec, map_nedelec);
//     M_s->FillComplete(map_nedelec, map_nedelec);    
//     log_matrix_stats("Matrix A_sorted", A_s);
//     log_matrix_stats("Matrix M_sorted", M_s);
    
    // Calculate Y and Y'
    pbe_start(61, "construction of Y");
    Y = new Epetra_FECrsMatrix(Copy, map_nedelec, false);
    mesh_mixed.constructY(*Y, map_lagrange, map_nedelec);

    Epetra_FECrsMatrix* Y_transp = new Epetra_FECrsMatrix (Copy, map_lagrange, false);
    mesh_mixed.constructY_transp(*Y_transp, map_nedelec, map_lagrange);

    pbe_stop(61);
    log_matrix_stats("Matrix Y", Y);
    log_matrix_stats("Matrix Y'", Y_transp);
    log_mem_footprint(Comm, "after constr of Y and Y'");
   
    // Passing true to arg2 of MatrixMatrix::Multiply should be avoided since it slows down the routine by a factor of ~10
    // This is the reason why Y_transp is constructed.

    // Calculate C = M * Y
    pbe_start(62, "construction of C and C^T");
    C = new Epetra_CrsMatrix (Copy, map_nedelec, 0);
    EpetraExt::MatrixMatrix::Multiply(*M, false, *Y, false, *C);
    pbe_stop(62);
    log_matrix_stats("Matrix C", C);
    log_mem_footprint(Comm, "after constr of C");

    // Calculate H = Y^T * C
    pbe_start(63, "construction of H");
    H = new Epetra_CrsMatrix (Copy, map_lagrange, 0);
#ifdef USE_PM_MATRIXMATRIX
    poor_man_matrix_matrix(*Y_transp, *C, *H);
#else
    EpetraExt::MatrixMatrix::Multiply(*Y_transp, false, *C, false, *H);
#endif
    pbe_stop(63);
    log_matrix_stats("Matrix H", H);
    log_mem_footprint(Comm, "after constr of H");

    if (params_.get<bool>("export_mtx")) {
        pbe_start(130, "MatrixMarket export");
        export_Epetra_CrsMatrix(*A, 15, "A.mtx");
        export_Epetra_CrsMatrix(*M, 15, "M.mtx");
        export_Epetra_CrsMatrix(*Y,  2, "Y.mtx");
        export_Epetra_CrsMatrix(*Y_transp,  2, "Y_transp.mtx");
        export_Epetra_CrsMatrix(*C, 15, "C.mtx");
        export_Epetra_CrsMatrix(*H, 15, "H.mtx");
        
        export_Epetra_Map(A->RowMap(), "A_rowmap.txt");
        export_Epetra_Map(H->RowMap(), "H_rowmap.txt");
        
//         export_Epetra_CrsMatrix(*A_s, 15, "A_s.mtx");
//         export_Epetra_CrsMatrix(*M_s, 15, "M_s.mtx");
        
        pbe_stop(130);
    }

    delete Y_transp;
    log_mem_footprint(Comm, "after deleting Y'");
}

const Epetra_Operator* LinearEigsolvOperators::getA() const {
    return(A);
}

const Epetra_Operator* LinearEigsolvOperators::getM() const {
    return(M);
}

const Epetra_Operator* LinearEigsolvOperators::getAsigmaPrec(double sigma) {
    if (Prec == 0) {
        rInfo("Calculate preconditioner for A-sigma*M...");
        pbe_start(70, "Preconditioner A-sigma*M construction");
        buildPreconditioner(sigma);
        pbe_stop(70);
    } else {
        assert(_sigma_K == sigma);
    }
    _sigma_K = sigma;
    return(Prec);
}

const Epetra_Operator* LinearEigsolvOperators::getHSolver() {
    if (h_solver_op_ == 0) {
        rInfo("Calculate solver for H...");
        pbe_start(72, "Solver H construction");
        Teuchos::ParameterList& h_solver(params_.sublist("h_solver"));
        string h_solver_type = h_solver.get<string>("type");
        if (h_solver_type == "pcg") {
            // block PCG solver
            // FIXME: These parameters should be included in the parameter list
            int max_iter_cg = 100;
            double h_solver_tol = 1e-10;
            BlockPCGSolver* pcg = new BlockPCGSolver(Comm, H, h_solver_tol, max_iter_cg, 0);
            pcg->setPreconditioner(const_cast<Epetra_Operator*>(getHPrec()));
            h_solver_op_ = pcg;
        } else if (h_solver_type == "lu") {
            // direct solver
            h_solver_op_ = new DirectSolverOperator(*H, params_.sublist("h_solver").sublist("amesos_params"));
        } else if (h_solver_type == "prec") {
            // just do one preconditioning step instead of a full solve
            h_solver_op_ = const_cast<Epetra_Operator*>(getHPrec());
        } else {
            rExit("Unknown solver for H: %s", h_solver_type.c_str());
        }
        pbe_stop(72);
    }
    return(h_solver_op_);
}

const Epetra_Operator* LinearEigsolvOperators::getHPrec() {
    if (HPrec == 0) {
        rInfo("Calculate preconditioner for H...");
        pbe_start(71, "Preconditioner H construction");
        Teuchos::ParameterList h_precon(params_.sublist("h_precon"));
        string h_precon_type = h_precon.get<string>("type");
        if (h_precon_type == "ml") {
            // Build ML preconditioner for positive-definite matrix H
            OutputCapturer cap;
            cap.begin_capture();
            HPrec = new ML_Epetra::MultiLevelPreconditioner(*H, params_.sublist("h_precon").sublist("ml_params"));
            cap.end_capture();
            cap.log("ML output:", RLOG_CHANNEL("debug"));
        } else if (h_precon_type == "diag") {
    
            // Preconditioner is diagonal of K = A -sigma*M
            HPrec = new DiagonalPreconditioner(*H);
    
        } else if (h_precon_type == "if") {
    
            // Preconditioner is incomplete LDLT-factorisation of K = A -sigma*M
        
            Teuchos::ParameterList List;
        
            // allocates an IFPACK factory. No data is associated
            // to this object (only method Create()).
            Ifpack Factory;
        
            // create the preconditioner. For valid PrecType values,
            // please check the documentation
            string PrecType = "IC stand-alone"; // incomplete LU
            int OverlapLevel = 1; // must be >= 0. If Comm.NumProc() == 1,
                                  // it is ignored.
        
            HPrec = Factory.Create(PrecType, H, OverlapLevel);
            assert(HPrec != 0);
        
            // specify parameters for IC
            List.set("fact: drop tolerance", 1e-9);
            List.set("fact: level-of-fill", 1);
        
            // sets the parameters
            int info;
            info = dynamic_cast<Ifpack_Preconditioner*>(HPrec)->SetParameters(List);
            rAssert(info >= 0);
        
            // initialize the preconditioner. At this point the matrix must
            // have been FillComplete()'d, but actual values are ignored.
            info = dynamic_cast<Ifpack_Preconditioner*>(HPrec)->Initialize();
            rAssert(info >= 0);
        
            // Builds the preconditioners, by looking for the values of
            // the matrix.
            info = dynamic_cast<Ifpack_Preconditioner*>(HPrec)->Compute();
            rAssert(info >= 0);
            
            // information about the preconditioner
            Ifpack_IC* ic_prec = dynamic_cast<Ifpack_IC*>(HPrec);
            if (ic_prec != 0)
                rInfo("#nonzeros of IF-preconditioner for H: %d", ic_prec->NumGlobalNonzeros());
            cout << *dynamic_cast<Ifpack_Preconditioner*>(HPrec);
        } else {
            rExit("Unknown preconditioner for H: %s", h_precon_type.c_str());
        }
        pbe_stop(71);
    }
    return(HPrec);
}

Epetra_Operator* LinearEigsolvOperators::getH() {
    return(H);
}

Epetra_Operator* LinearEigsolvOperators::getY() const {
    return(Y);
}

Epetra_Operator* LinearEigsolvOperators::getC() const {
    return(C);
}

---