src/eigsolv/QuadraticEigsolvOperators.cpp

Go to the documentation of this file.

#define max(a,b) ((a) >= (b) ? (a) : (b))

#include "myrlog.h"

#include "ml_include.h"
#include "EpetraExt_MatrixMatrix.h"
#include "Teuchos_ParameterList.hpp"
#include "ml_epetra_preconditioner.h"

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

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

#include "QuadraticEigsolvOperators.h"

QuadraticEigsolvOperators::QuadraticEigsolvOperators(FemaxMesh& mesh_mixed, Teuchos::ParameterList& params, const Epetra_Comm& MyComm) :
    Comm(MyComm),
    A(0),
    M(0),
    K(0),
    K11(0),
    K12(0),
    K22(0),
    H(0),
    H11(0),
    H12(0),
    H22(0),
    Y(0), 
    Y11(0),
    C(0),
    Prec(0),
    h_solver_op_(0),
    HPrec(0),
    _sigma_K(0.0),
    params_(params),
    Aprec(0),
    AprecSolver(0),
    non(0),
    noe(0)
{
    non = mesh_mixed.get_lagrange_mesh().get_num_global_mdofs(1);
    noe = mesh_mixed.get_nedelec_mesh().get_num_global_mdofs(1);

    matrixAssembly(mesh_mixed);
}

QuadraticEigsolvOperators::~QuadraticEigsolvOperators() {
    OutputCapturer capt;
    capt.begin_capture();
    
    delete A;
    delete M;
    delete K;
    delete K11;
    delete K12;
    delete K22;
    delete Y;
    delete H;
    delete H11;
    delete H12;
    delete H22;
    delete C;
    delete h_solver_op_;
    delete HPrec;
    delete Prec;
    delete Aprec;
    delete AprecSolver;
        
    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 QuadraticEigsolvOperators::buildPreconditioner(double sigma) {
    // construct K = A - sigma*M
    assert(K == 0);
    K = build_Asigma(*A, *M, sigma);

    // 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);

    //  AprecSolver.SetAztecOption(AZ_precond, AZ_dom_decomp);
    //  AprecSolver.SetAztecOption(AZ_ilu, AZ_subdomain_solve);

    // 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
    int iterAztec = 2;
    Prec = new AztecOO_Operator(AprecSolver, iterAztec);

    return(0);
}

void QuadraticEigsolvOperators::matrixAssembly(FemaxMesh& mesh_mixed) {

    NedelecMesh& mesh_nedelec = mesh_mixed.get_nedelec_mesh();
    LagrangeMesh& mesh_lagrange = mesh_mixed.get_lagrange_mesh();

    // Create Epetra_Maps
    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);
    rInfoAll("map_nedelec: owned_dofs.size()=%u", static_cast<unsigned int>(owned_dofs.size()));
    
    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);

    //
    // 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");

    // Calculate 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");

        pbe_stop(130);
    }

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

int QuadraticEigsolvOperators::buildHierarchicalPrec(double sigma) {
    // We conservatively assume, that ML will be used to solve the (1,1)-block. In 
    // doing so we guarantee, that Y11 and H11 are constructed. This will make sure that 
    // Maxwell ML will be used in HierarchicalBasisPrec.
    bool is_ml = true;

    const Epetra_Map& map_nedelec = M->OperatorDomainMap();
    int sizeLocal = map_nedelec.NumMyElements();

    int nBlockRows1 = 0;
    int nBlockRows2 = 0;
    int i = 0;
    int gid = 0;

    for (i = 0; i < sizeLocal; ++i) {
        if (map_nedelec.GID(i) < noe)
            nBlockRows1++;
        else
            nBlockRows2++;
    }

    int *myBlockRow1 = new int[nBlockRows1];
    int *myBlockRow2 = new int[nBlockRows2];

    nBlockRows1 = 0;
    nBlockRows2 = 0;
    for (i = 0; i < sizeLocal; ++i) {
        gid = map_nedelec.GID(i);
        if (gid < noe)
            myBlockRow1[nBlockRows1++] = gid;
        else
            myBlockRow2[nBlockRows2++] = gid-noe;
    }

    Epetra_Map map1(noe, nBlockRows1, myBlockRow1, 0, Comm);
    Epetra_Map map2(map_nedelec.NumGlobalElements() - noe, 
                    nBlockRows2, myBlockRow2, 0, Comm);

    K11 = new Epetra_CrsMatrix(Copy, map1, 0);
    K12 = new Epetra_CrsMatrix(Copy, map1, 0);
    K22 = new Epetra_CrsMatrix(Copy, map2, 0);
    if (is_ml) 
        Y11 = new Epetra_CrsMatrix(Copy, map1, 0);

    int numEntries = 0;
    int *indices = 0;
    double *values = 0;
    int info = 0;
    int j = 0;
    int col = 0;
    double val = 0.0;

    int numEntriesY = 0;
    int *indicesY = 0;
    double *valuesY = 0;

    for (i = 0; i < sizeLocal; ++i) {
        gid = map_nedelec.GID(i);

        if (gid < noe) { 

            assert(A->ExtractMyRowView(i, numEntries, values, indices) >= 0 );
            for (j = 0; j < numEntries; ++j) {
                col = A->GCID(indices[j]);
                if (col < noe) {
                    assert(K11->InsertGlobalValues(gid, 1, values + j, &col) >=0);
                }
                else {
                    col -= noe;
                    assert(K12->InsertGlobalValues(gid, 1, values + j, &col) >=0);
                }
            }

            assert(M->ExtractMyRowView(i, numEntries, values, indices) >= 0);
            for (j = 0; j < numEntries; ++j) {
                col = M->GCID(indices[j]);
                val = - sigma*values[j];
                if (col < noe) {
                    assert((info = K11->SumIntoGlobalValues(gid, 1, &val, &col)) >= 0);
                    if (info > 0) 
                        assert(K11->InsertGlobalValues(gid, 1, &val, &col) >= 0);
                }
                else {
                    col -= noe;
                    assert((info = K12->SumIntoGlobalValues(gid, 1, &val, &col)) >= 0);
                    if (info > 0) 
                        assert(K12->InsertGlobalValues(gid, 1, &val, &col) >= 0);
                }
            }

            if (is_ml) {
                assert(Y->ExtractMyRowView(i, numEntriesY, valuesY, indicesY) >= 0 );
                for (j = 0; j < numEntriesY; ++j) {
                    col = Y->GCID(indicesY[j]);
                    if (col < noe) {
                        assert(Y11->InsertGlobalValues(gid, 1, valuesY + j, &col) >=0);
                    }
                }
            }

        }
        
        else {

            gid -= noe;

            assert(A->ExtractMyRowView(i, numEntries, values, indices) >=0);
            for (j = 0; j < numEntries; ++j) {
                col = A->GCID(indices[j]);
                //              if (col >= noe) {
                if (col >= noe && A->MyGRID(col)) {
                    col -= noe;
                    assert(K22->InsertGlobalValues(gid, 1, values + j, &col) >=0);
                }
            }

            assert(M->ExtractMyRowView(i, numEntries, values, indices) >= 0);
            for (j = 0; j < numEntries; ++j) {
                col = M->GCID(indices[j]);
                //              if (col >= noe) {
                // FIXME: Only local diagonal blocks stored in K22!
                if (col >= noe && A->MyGRID(col)) {
                    col -= noe;
                    val = - sigma*values[j];
                    assert((info = K22->SumIntoGlobalValues(gid, 1, &val, &col)) >= 0);
                    if (info > 0)
                        assert(K22->InsertGlobalValues(gid, 1, &val, &col) >= 0);
                }
            }
        }

    }

    K11->FillComplete(map1, map1);
    K12->FillComplete(map2, map1);
    K22->FillComplete(map2, map2);

    /*
      int numProc = Comm.NumProc();
      int myPid = Comm.MyPID();
      int* sizes = new int[numProc];
      sizes[myPid] = map2.NumMyElements();
      Comm.GatherAll(sizes+myPid, sizes, 1);
      int* delim = new int[numProc+1];
      delim[0] = 0;
      for (i=0; i<numProc; i++)
      delim[i+1] = delim[i]+sizes[i];

      if ( myPid==0 ) {
      cout << "Delimiters : ";
      for (i=0; i<=numProc; i++)
      cout << delim[i] << " ";
      cout << endl;
      }
      Comm.Barrier();

      int k = 0;
      int diag = 0;
      int numDiag = 0;
      int numOutDiag = 0;
      int in = 0;
      int out = 0;
      for (i=0; i<map2.NumMyElements(); i++) {
      gid = map2.GID(i);

      for (diag=1; gid >= delim[diag]; diag++);
      K22->ExtractMyRowView(i, numEntries, values, indices);
      for (j = 0; j < numEntries; ++j) {
      col = K22->GCID(indices[j]);
      for (k=1; col >= delim[k]; k++);
      if (diag == k) numDiag++; else numOutDiag++;
      }

      if ( gid >= delim[myPid] && gid < delim[myPid+1] )  
      in++;
      else
      out++;
      }

      cout << myPid << "  in = " << in << "   out = " << out << endl;

      cout << myPid << " numDiag = " << numDiag << "  numOutDiag = " << numOutDiag << endl;

      delete delim;
      delete sizes;
    */

    if (is_ml) {
        if (H11 == 0) 
            H11 = getH11();
        Y11->FillComplete(H11->RangeMap(), map1);
    }
    Prec = new HierarchicalBasisPrec(K11, Y11, H11, K12, K22, M->OperatorDomainMap(), M->OperatorRangeMap(), params_.sublist("asigma_precon"));


    /*
      K = new Epetra_CrsMatrix(Copy, map_nedelec, 0);
      sizeLocal = map_nedelec.NumMyElements();
      for (i = 0; i < sizeLocal; ++i) {
      int numEntries = 0;
      int *indices;
      double *values;
      int iMap = map_nedelec.GID(i);
      int info;
      info = A->ExtractMyRowView(i, numEntries, values, indices);
      assert(info >= 0);
      int jj;
      for (jj = 0; jj < numEntries; ++jj) {
      int col = A->GCID(indices[jj]);
      assert(K->InsertGlobalValues(iMap, 1, values + jj, &col) >=0);
      }
      info = M->ExtractMyRowView(i, numEntries, values, indices);
      assert(info >= 0);
      for (jj = 0; jj < numEntries; ++jj) {
      int col = M->GCID(indices[jj]);
      double val = - sigma*values[jj];
      info = K->SumIntoGlobalValues(iMap, 1, &val, &col);
      assert(info >= 0);
      if (info > 0) {
      assert(K->InsertGlobalValues(iMap, 1, &val, &col) >= 0);
      }
      }
      }
      K->FillComplete(map_nedelec, map_nedelec);
    */

    /*

    if (Comm.MyPID() == 0) {
    cout << "TEST K :......................" << endl;
    }
    Epetra_MultiVector X(map_nedelec, 1, false);
    Epetra_MultiVector Y1(map_nedelec, 1, false);
    Epetra_MultiVector Y2(map_nedelec, 1, false);
    X.Random();


    Epetra_LinearProblem Problem1;
    Problem1.SetOperator( K );
    Problem1.SetLHS( &Y1 );
    Problem1.SetRHS( &X );

    AztecOO Solver(Problem1);
    Solver.SetAztecOption( AZ_solver,  AZ_gmres );
    //    Solver.SetAztecOption( AZ_output, AZ_none );
    Solver.Iterate(100, 1e-10);

    if (Comm.MyPID() == 0) {
    cout << "TEST Prec :......................" << endl;
    }

    Prec->ApplyInverse(X, Y2);

    if (Comm.MyPID() == 0) {
    cout << "TEST Prec done :......................" << endl;
    }

    double nrm1 = 0.0;
    double nrm2 = 0.0;
    Y1.Norm2(&nrm1);
    Y2.Norm2(&nrm2);
    if (Comm.MyPID() == 0) {
    cout << "Norm(Y1) = " << nrm1 << endl;
    cout << "Norm(Y2) = " << nrm2 << endl;
    }    
    exit(0);
    */

    /*
      Epetra_Vector X(map_nedelec, false);
      Epetra_Vector Y1(map_nedelec, false);
      Epetra_Vector Y2(map_nedelec, false);
      X.Random();
      K->Multiply(false, X, Y1);
      Prec->Apply(X, Y2);
      double nrm1 = 0.0;
      double nrm2 = 0.0;
      Y1.Norm2(&nrm1);
      Y2.Norm2(&nrm2);
      if (Comm.MyPID() == 0) {
      cout << "Norm(Y1) = " << nrm1 << endl;
      cout << "Norm(Y2) = " << nrm2 << endl;
      }    
      exit(0);
    */    

    delete myBlockRow1;
    delete myBlockRow2;
    return(0);
}

int QuadraticEigsolvOperators::buildHierarchicalHPrec() {

    //
    // construct H = | H11 H12 |
    //               | H12 H22 |
    //
    // dim(H11) = non    dim(H22) = noe 

    const Epetra_Map& map_lagrange = H->OperatorDomainMap();
    int sizeLocal = map_lagrange.NumMyElements();

    int nBlockRows1 = 0;
    int nBlockRows2 = 0;
    int i = 0;
    int gid = 0;

    for (i = 0; i < sizeLocal; ++i) {
        if (map_lagrange.GID(i) < non)
            nBlockRows1++;
        else
            nBlockRows2++;
    }

    int *myBlockRow1 = new int[nBlockRows1];
    int *myBlockRow2 = new int[nBlockRows2];

    nBlockRows1 = 0;
    nBlockRows2 = 0;
    for (i = 0; i < sizeLocal; ++i) {
        gid = map_lagrange.GID(i);
        if (gid < non)
            myBlockRow1[nBlockRows1++] = gid;
        else
            myBlockRow2[nBlockRows2++] = gid-non;
    }

    Epetra_Map map1(non, nBlockRows1, myBlockRow1, 0, Comm);
    Epetra_Map map2(noe, nBlockRows2, myBlockRow2, 0, Comm);

    bool create_h11 = (H11 == 0) ? true : false;
    if (create_h11) 
        H11 = new Epetra_CrsMatrix(Copy, map1, 0);
    H12 = new Epetra_CrsMatrix(Copy, map1, 0);
    H22 = new Epetra_CrsMatrix(Copy, map2, 0);

    int numEntries = 0;
    int *indices = 0;
    double *values = 0;
    int j = 0;
    int col = 0;

    for (i = 0; i < sizeLocal; ++i) {
        gid = map_lagrange.GID(i);
        if (gid < non) {
            assert(H->ExtractMyRowView(i, numEntries, values, indices) >= 0 );
            for (j = 0; j < numEntries; ++j) {
                col = H->GCID(indices[j]);
                if (col < non) {
                    if (create_h11)
                        assert(H11->InsertGlobalValues(gid, 1, values + j, &col) >=0);
                }
                else {
                    col -= non;
                    assert(H12->InsertGlobalValues(gid, 1, values + j, &col) >=0);
                }
            }
        }
        else {
            gid -= non;
            assert(H->ExtractMyRowView(i, numEntries, values, indices) >=0);
            for (j = 0; j < numEntries; ++j) {
                col = H->GCID(indices[j]);
                //              if (col >= non) {
                // FIXME: Only local diagonal blocks stored in H22!
                if (col >= non && H->MyGRID(col)) {
                    col -= non;
                    assert(H22->InsertGlobalValues(gid, 1, values + j, &col) >=0);
                }
            }
        }
    }

    if (create_h11)
        H11->FillComplete(map1, map1);
    H12->FillComplete(map2, map1);
    H22->FillComplete(map2, map2);

    HPrec = new HierarchicalBasisPrec(H11, 0, 0, H12, H22, H->OperatorDomainMap(), H->OperatorRangeMap(), params_.sublist("h_precon"));
    
    delete myBlockRow1;
    delete myBlockRow2;
    return(0);
}


Epetra_CrsMatrix* QuadraticEigsolvOperators::getH11() {

    if (H11 == 0) {
        const Epetra_Map& map_lagrange = H->OperatorDomainMap();
        int sizeLocal = map_lagrange.NumMyElements();

        int nBlockRows1 = 0;
        int i = 0;
        int gid = 0;

        for (i = 0; i < sizeLocal; ++i) {
            if (map_lagrange.GID(i) < non)
                nBlockRows1++;
        }

        int *myBlockRow1 = new int[nBlockRows1];

        nBlockRows1 = 0;
        for (i = 0; i < sizeLocal; ++i) {
            gid = map_lagrange.GID(i);
            if (gid < non)
                myBlockRow1[nBlockRows1++] = gid;
        }

        Epetra_Map map1(non, nBlockRows1, myBlockRow1, 0, Comm);
        H11 = new Epetra_CrsMatrix(Copy, map1, 0);


        int numEntries = 0;
        int *indices = 0;
        double *values = 0;
        int j = 0;
        int col = 0;

        for (i = 0; i < sizeLocal; ++i) {
            gid = map_lagrange.GID(i);
            if (gid < non) {
                assert(H->ExtractMyRowView(i, numEntries, values, indices) >= 0 );
                for (j = 0; j < numEntries; ++j) {
                    col = H->GCID(indices[j]);
                    if (col < non) {
                        assert(H11->InsertGlobalValues(gid, 1, values + j, &col) >=0);
                    }
                }
            }
        }

        H11->FillComplete(map1, map1);
    }

    return(H11);
}

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

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

const Epetra_Operator* QuadraticEigsolvOperators::getAsigmaPrec(double sigma) {
    if (Prec == 0) {
        Teuchos::ParameterList& sublist(params_.sublist("asigma_precon"));
        string asigma_precon_type = sublist.get<string>("type");
        rInfo("Calculate preconditioner for A-sigma*M...");
        pbe_start(70, "Preconditioner A-sigma*M construction");
        if (asigma_precon_type == "neumann") {
            buildPreconditioner(sigma);
        } else if (asigma_precon_type == "2level") {
            buildHierarchicalPrec(sigma);
        } else if (asigma_precon_type == "lu") {
            // Preconditioner is LU-factorisation of K = A -sigma*M
            assert(K == 0);
            K = build_Asigma(*A, *M, sigma);
            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 {
            rExit("Unknown preconditioner for A - sigma*M: %s", asigma_precon_type.c_str());
        }
        pbe_stop(70);
    } else {
        assert(_sigma_K == sigma);
    }
    _sigma_K = sigma;
    return(Prec);
}

const Epetra_Operator* QuadraticEigsolvOperators::getHSolver() {
    if (h_solver_op_ == 0) {
        rInfo("Calculate solver for H...");
        pbe_start(72, "Solver H construction");
        Teuchos::ParameterList& sublist(params_.sublist("h_solver"));
        string h_solver_type = sublist.get<string>("type");
        if (h_solver_type == "pcg") {
            // block PCG solver
            // Fixme: parametrize this....
            int max_iter_cg = 100;
            double h_solver_tol = 1e-10;
            int verbosity = 0;
            BlockPCGSolver* pcg = new BlockPCGSolver(Comm, H, h_solver_tol, max_iter_cg, verbosity);
            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* QuadraticEigsolvOperators::getHPrec() {
    if (HPrec == 0) {
        rInfo("Calculate preconditioner for H...");
        pbe_start(71, "Preconditioner H construction");
        Teuchos::ParameterList& sublist(params_.sublist("h_precon"));
        string h_precon_type = sublist.get<string>("type");
        if (h_precon_type == "2level") {
            // 2-level hierarchical preconditioner
            buildHierarchicalHPrec();
        } else if (h_precon_type == "ml") {
            // 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"), true);
            cap.end_capture();
            cap.log("ML output:", RLOG_CHANNEL("debug"));
        } else {
            rExit("Unknown preconditioner for H: %s", h_precon_type.c_str());
        }
        pbe_stop(71);
    }
    return(HPrec);
}

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

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

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

---