<commit_before>/** @file gsTrilinosSolvers.cpp

    @brief Wrappers for Trilinos solvers

    This file is part of the G+Smo library.

    This Source Code Form is subject to the terms of the Mozilla Public
    License, v. 2.0. If a copy of the MPL was not distributed with this
    file, You can obtain one at http://mozilla.org/MPL/2.0/.

    Author(s): <NAME>
*/

#include <gsTrilinos/gsTrilinosSolvers.h>
#include "gsTrilinosHeaders.h"

#include "Amesos.h"
#include "Amesos_BaseSolver.h"
#include "AztecOO.h"
#include "AztecOO_Version.h"

	 	 
//#include "Ifpack_ConfigDefs.h" 
//#include "Ifpack.h" 
//#include "Ifpack_AdditiveSchwarz.h" 

//#include "Amesos_Superlu.h"

#include "gsTrilinosHeaders.h"
#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosEpetraAdapter.hpp"
#include "BelosBlockCGSolMgr.hpp"

namespace gismo
{

namespace trilinos
{

namespace solver
{

struct AbstractSolverPrivate
{
    Epetra_LinearProblem Problem;
    Vector solution;
};

AbstractSolver::AbstractSolver() : my(NULL)
{

}


AbstractSolver::AbstractSolver(const SparseMatrix & A)
: my(new AbstractSolverPrivate)
{
    my->Problem.SetOperator(A.get());
    my->solution.setFrom(A); // i.e. A.get()->OperatorDomainMap()
    my->Problem.SetLHS(my->solution.get());
}

AbstractSolver::~AbstractSolver()
{
    delete my;
}

const Vector & AbstractSolver::solve(const Vector & b)
{
    my->Problem.SetRHS(b.get());
    solveProblem(); // virtual call
    return my->solution; 
}

void AbstractSolver::getSolution( gsVector<real_t> & sol, const int rank) const
{
    my->solution.copyTo(sol, rank);
}

void GMRES::solveProblem()
{
    /* 1. Preconditionner */
/*    
    // allocates an IFPACK factory. No data is associated 
    Ifpack Factory;

    // create the preconditioner. For valid PrecType values,
    // please check the documentation
    std::string PrecType = "ILU"; // incomplete LU
    int OverlapLevel = 1; // must be >= 0. ignored for Comm.NumProc() == 1
    Teuchos::RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*A, OverlapLevel) );
    assert(Prec != Teuchos::null);
    
    // specify parameters for ILU
    List.set("fact: drop tolerance", 1e-9);
    List.set("fact: level-of-fill", 1);
    // the combine mode is on the following:
    // "Add", "Zero", "Insert", "InsertAdd", "Average", "AbsMax"
    // Their meaning is as defined in file Epetra_CombineMode.h   
    List.set("schwarz: combine mode", "Add");
    // sets the parameters
    IFPACK_CHK_ERR(Prec->SetParameters(List));
    
    // initialize the preconditioner. At this point the matrix must
    // have been FillComplete()'d, but actual values are ignored.
    IFPACK_CHK_ERR(Prec->Initialize());
    
    // Builds the preconditioners, by looking for the values of 
    // the matrix.
    IFPACK_CHK_ERR(Prec->Compute());
    
*/
    /* 2. AztecOO solver / GMRES*/
    
    AztecOO Solver;
    Solver.SetProblem(my->Problem);
    Solver.SetAztecOption(AZ_solver, AZ_gmres);
    Solver.SetAztecOption(AZ_output,AZ_none);//32
    //Solver.SetPrecOperator(Prec);
    Solver.Iterate(m_maxIter, m_tolerance);
}

void KLU::solveProblem()
{
    static const char * SolverType = "Klu";
    static Amesos Factory;
    const bool IsAvailable = Factory.Query(SolverType);
    GISMO_ENSURE(IsAvailable, "Amesos KLU is not available.\n");
    
    Amesos_BaseSolver * Solver = Factory.Create(SolverType, my->Problem);
    Solver->SymbolicFactorization();
    Solver->NumericFactorization();
    Solver->Solve();
}

void SuperLU::solveProblem()
{
//    Amesos_Superlu Solver(my->Problem);
//    Solver.SymbolicFactorization();
//    Solver.NumericFactorization();
//    Solver.Solve();
}



/*   --- Belos--- */

struct BelosSolverPrivate
{
    Belos::LinearProblem<real_t,Epetra_MultiVector,Epetra_Operator> Problem;
};

BelosSolver::BelosSolver(const SparseMatrix & A
                         //, const std::string solver_teuchosUser
    )
: Base(A), myBelos(new BelosSolverPrivate), blocksize(1), maxiters(500)
{ 
    // Note: By default the string variable SolverTeuchosUser = "Belos". 
    // This can be adapted to get different values if other solvers with 
    // similar implementation structures are considered later on. 
    //SolverTeuchosUser = solver_teuchosUser; 

    myBelos->Problem.setOperator(A.getRCP());
    my->solution.setFrom(A); // i.e. A.get()->OperatorDomainMap()
    myBelos->Problem.setLHS(my->solution.getRCP());
    
	//!!!!!!// If the matrix is symmetric, specify this in the linear problem. 
	//!!!!!!my->Problem->setHermitian(); 
}

BelosSolver::~BelosSolver()
{
    delete myBelos;
}


void BelosSolver::solveProblem()
{
    // Grab right-hand side
    myBelos->Problem.setRHS( Teuchos::rcp(my->Problem.GetRHS(), false) );

    // Tell the program that setting of the linear problem is done. 
    // Throw an error if failed. 
    bool err_set = myBelos->Problem.setProblem(); 
    
    GISMO_ASSERT(true == err_set, "Error: Belos Problem couldn't be" 
                 " initialized."); 

    // Teuchos::ScalarTraits<double>::magnitudeType tol = 1.0e-5;

    double tol = 1.0e-5;

    Teuchos::ParameterList belosList;
    belosList.set( "Block Size", blocksize );         // Blocksize to be used by iterative solver
    belosList.set( "Maximum Iterations", maxiters );  // Maximum number of iterations allowed<commit_msg>minor<commit_after>/** @file gsTrilinosSolvers.cpp

    @brief Wrappers for Trilinos solvers

    This file is part of the G+Smo library.

    This Source Code Form is subject to the terms of the Mozilla Public
    License, v. 2.0. If a copy of the MPL was not distributed with this
    file, You can obtain one at http://mozilla.org/MPL/2.0/.

    Author(s): <NAME>
*/

#include <gsTrilinos/gsTrilinosSolvers.h>
#include "gsTrilinosHeaders.h"

#include "Amesos.h"
#include "Amesos_BaseSolver.h"
#include "AztecOO.h"
#include "AztecOO_Version.h"

	 	 
//#include "Ifpack_ConfigDefs.h" 
//#include "Ifpack.h" 
//#include "Ifpack_AdditiveSchwarz.h" 

//#include "Amesos_Superlu.h"

#include "gsTrilinosHeaders.h"
#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosEpetraAdapter.hpp"
#include "BelosBlockCGSolMgr.hpp"

namespace gismo
{

namespace trilinos
{

namespace solver
{

struct AbstractSolverPrivate
{
    Epetra_LinearProblem Problem;
    Vector solution;
};

AbstractSolver::AbstractSolver() : my(NULL)
{

}


AbstractSolver::AbstractSolver(const SparseMatrix & A)
: my(new AbstractSolverPrivate)
{
    my->Problem.SetOperator(A.get());
    my->solution.setFrom(A); // i.e. A.get()->OperatorDomainMap()
    my->Problem.SetLHS(my->solution.get());
}

AbstractSolver::~AbstractSolver()
{
    delete my;
}

const Vector & AbstractSolver::solve(const Vector & b)
{
    my->Problem.SetRHS(b.get());
    solveProblem(); // virtual call
    return my->solution; 
}

void AbstractSolver::getSolution( gsVector<real_t> & sol, const int rank) const
{
    my->solution.copyTo(sol, rank);
}

void GMRES::solveProblem()
{
    /* 1. Preconditionner */
/*    
    // allocates an IFPACK factory. No data is associated 
    Ifpack Factory;

    // create the preconditioner. For valid PrecType values,
    // please check the documentation
    std::string PrecType = "ILU"; // incomplete LU
    int OverlapLevel = 1; // must be >= 0. ignored for Comm.NumProc() == 1
    Teuchos::RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*A, OverlapLevel) );
    assert(Prec != Teuchos::null);
    
    // specify parameters for ILU
    List.set("fact: drop tolerance", 1e-9);
    List.set("fact: level-of-fill", 1);
    // the combine mode is on the following:
    // "Add", "Zero", "Insert", "InsertAdd", "Average", "AbsMax"
    // Their meaning is as defined in file Epetra_CombineMode.h   
    List.set("schwarz: combine mode", "Add");
    // sets the parameters
    IFPACK_CHK_ERR(Prec->SetParameters(List));
    
    // initialize the preconditioner. At this point the matrix must
    // have been FillComplete()'d, but actual values are ignored.
    IFPACK_CHK_ERR(Prec->Initialize());
    
    // Builds the preconditioners, by looking for the values of 
    // the matrix.
    IFPACK_CHK_ERR(Prec->Compute());
    
*/
    /* 2. AztecOO solver / GMRES*/
    
    AztecOO Solver;
    Solver.SetProblem(my->Problem);
    Solver.SetAztecOption(AZ_solver, AZ_gmres);
    Solver.SetAztecOption(AZ_output,AZ_none);//32
    //Solver.SetPrecOperator(Prec);
    Solver.Iterate(m_maxIter, m_tolerance);
}

void KLU::solveProblem()
{
    static const char * SolverType = "Klu";
    static Amesos Factory;
    const bool IsAvailable = Factory.Query(SolverType);
    GISMO_ENSURE(IsAvailable, "Amesos KLU is not available.\n");
    
    Amesos_BaseSolver * Solver = Factory.Create(SolverType, my->Problem);
    Solver->SymbolicFactorization();
    Solver->NumericFactorization();
    Solver->Solve();
}

void SuperLU::solveProblem()
{
//    Amesos_Superlu Solver(my->Problem);
//    Solver.SymbolicFactorization();
//    Solver.NumericFactorization();
//    Solver.Solve();
}



/*   --- Belos--- */

struct BelosSolverPrivate
{
    Belos::LinearProblem<real_t,Epetra_MultiVector,Epetra_Operator> Problem;
};

BelosSolver::BelosSolver(const SparseMatrix & A
                         //, const std::string solver_teuchosUser
    )
: Base(A), myBelos(new BelosSolverPrivate), blocksize(1), maxiters(500)
{ 
    // Note: By default the string variable SolverTeuchosUser = "Belos". 
    // This can be adapted to get different values if other solvers with 
    // similar implementation structures are considered later on. 
    //SolverTeuchosUser = solver_teuchosUser; 

    myBelos->Problem.setOperator(A.getRCP());
    myBelos->Problem.setLHS(my->solution.getRCP());
    
	//!!!!!!// If the matrix is symmetric, specify this in the linear problem. 
	//!!!!!!my->Problem->setHermitian(); 
}

BelosSolver::~BelosSolver()
{
    delete myBelos;
}


void BelosSolver::solveProblem()
{
    // Grab right-hand side
    myBelos->Problem.setRHS( Teuchos::rcp(my->Problem.GetRHS(), false) );

    // Tell the program that setting of the linear problem is done. 
    // Throw an error if failed. 
    bool err_set = myBelos->Problem.setProblem(); 
    
    GISMO_ASSERT(true == err_set, "Error: Belos Problem couldn't be" 
                 " initialized."); 

    // Teuchos::ScalarTraits<double>::magnitudeType tol = 1.0e-5;

    double tol = 1.0e-5;

    Teuchos::ParameterList belosList;
    belosList.set( "Block Size", blocksize );         // Blocksize to be used by iterative solver
    belosList.set( "Maximum Iterations", maxiters );  // Maximum number of iterations allowed
{'fraction_non_alphanumeric': 0.1321656050955414, 'fraction_numerical': 0.004496065942300487, 'mean_word_length': 3.4028865979381444, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 69, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5001288', 'n_tokens_mistral': 3783, 'n_tokens_neox': 3558, 'n_words': 842}