[Trilinos-Users] LSQR Solver in BelosThyra

[Einar Otnes]

Dear experts,

Does the BelosThyra adapter support the recently implemented LSQR solver. I have been trying to make the LSQR Solver work on a simple 2x2 example, to no avail. The CG and GMRES solvers are working as expected, but not the LSQR. The example code I'm using is listed below. Any thoughts or comments ?

Thanks,

Einar Otnes
 




/***********************************************************************
 *
 * A program to test Belos with Thyra objects
 *
 ***********************************************************************/
#include <BelosTpetraAdapter.hpp>
#include <Teuchos_ParameterList.hpp>
#include <Teuchos_RCP.hpp>

#include <Thyra_VectorBase.hpp>
#include <Thyra_VectorSpaceBase.hpp>
#include <Thyra_MultiVectorBase.hpp>


#include "BelosThyraAdapter.hpp"
#include <BelosSolverFactory.hpp>
#include "BelosBlockCGSolMgr.hpp"
#include "BelosLSQRSolMgr.hpp"
#include "BelosBlockGmresSolMgr.hpp"

#include "Thyra_TpetraThyraWrappers.hpp"

int main(int argc, char* argv[]) {

	typedef int OT;
	typedef double ST;
	typedef Thyra::LinearOpBase<ST>           OP;
	typedef Thyra::VectorBase<ST>              V;
	typedef Thyra::MultiVectorBase<ST>        MV;
	typedef Teuchos::ScalarTraits<ST    >    SCT;
	typedef SCT::magnitudeType                MT;
	typedef Belos::LinearProblem<ST, MV, OP>  problem_type;


	// Want to solve the system Ax=b

	// Set problem size
	OT m = 2;
	OT n = 2;

	// Allocate trace arrays, input and output

	Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > rngSp = Thyra::defaultSpmdVectorSpace<ST>(m);
	Teuchos::RCP<const Thyra::VectorSpaceBase<ST> > domSp = Thyra::defaultSpmdVectorSpace<ST>(n);

	Teuchos::RCP<Teuchos::FancyOStream>
	out = Teuchos::VerboseObjectBase::getDefaultOStream();


	Teuchos::RCP<MV> A = Thyra::createMembers(rngSp,domSp);
	Teuchos::RCP<V>  x = Thyra::createMember(A->domain());
	Teuchos::RCP<V>  b = Thyra::createMember(A->range());

	// Initialize input data

	// Create Identity Matrix

	Thyra::assign(A.ptr(),SCT::zero());

	for(int i=0; i< std::min(m,n);++i){
		Thyra::set_ele(i,SCT::one(),A->col(i).ptr());
	}


	// Set value for solution, x=1
	Thyra::put_scalar(SCT::one(),x.ptr());

	std::cout << "Correct Solution" << std::endl;
	x->describe(*out,Teuchos::VERB_EXTREME);




	// Calculate b=Ax to obtain the proper input to the inversion.
	A->apply(Thyra::NOTRANS,*x,b.ptr(),SCT::one(),SCT::zero());

	// Insert Belos Solver here...


	MT tol = 1.0e-8;            // relative residual tolerance
	MT lambda = 0.000;            // Preconditioner for LSQR

	Teuchos::RCP<Teuchos::ParameterList>  belosList =  Teuchos::parameterList();

	belosList->set( "Lambda", lambda );                     // Relative convergence tolerance requested
	belosList->set( "Convergence Tolerance", tol );         // Relative convergence tolerance requested




	// Initial guess for solution, x=0. The solution should give x=1.
	Thyra::put_scalar(SCT::zero(),x.ptr());

	Teuchos::RCP<problem_type> problem = rcp (new problem_type(A,x,b));

	problem->setProblem();



	// Create solver manager using factory.
	Belos::SolverFactory<ST, MV, OP > factory;

	// Select a solver
//	Teuchos::RCP<Belos::SolverManager<ST, MV, OP> > solver = factory.create ("GMRES", belosList);
	Teuchos::RCP<Belos::SolverManager<ST, MV, OP> > solver = factory.create ("CG", belosList);
//	Teuchos::RCP<Belos::SolverManager<ST, MV, OP> > solver = factory.create ("LSQR", belosList);

	solver->setProblem(problem);

	// Solve
	solver->solve();

	std::cout << "Calculated Solution" << std::endl;
	x->describe(*out,Teuchos::VERB_EXTREME);

	return EXIT_SUCCESS;
}