[Trilinos-Users] Complex support using Belos LSQRSolMgr

[Einar Otnes]

Dear experts,
I have been trying to use Belos LSQR solver for a complex valued system,
Ax=b, but I'm not able to compile it.

Using a scalar type that is real valued, typedef double ST, the example
attached compiles and runs OK. Setting the scalar type to be complex, i.e.
typedef std::complex<double> ST, does not compile. The example I've
provided is based on one of the examples presented by Heidi Thornquist in
TUG 2011.

So, in short, does Belos LSQR solver support complex types at this stage?

The version I'm running is Trilinos 10.12.2, compiled using gcc-4.6.1

Thanks for all your help.

Sincerely,

Einar Otnes
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20120903/637eea6f/attachment.html 
-------------- next part --------------

/***********************************************************************
 *
 * A program to test Belos with Thyra objects
 *
 ***********************************************************************/

#include "Teuchos_ParameterList.hpp"
#include "Thyra_TpetraThyraWrappers.hpp"
#include "Thyra_VectorBase.hpp"
#include "Thyra_LinearOpBase.hpp"
#include "Thyra_OperatorVectorClientSupport.hpp"
#include "BelosConfigDefs.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosThyraAdapter.hpp"
#include "BelosLSQRSolMgr.hpp"

//#define COMPLEX // Turn on to test with complex type


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

#ifdef COMPLEX
	// If I set ST as complex, example does not compile...
	typedef std::complex<double> ST;
#else
	// If I set ST as real, everything is fine...
	typedef double ST;
#endif






	typedef int OT;
	typedef Thyra::LinearOpBase<ST>           OP;
	typedef Thyra::VectorBase<ST>              V;
	typedef Thyra::MultiVectorBase<ST>        MV;
	typedef Teuchos::ScalarTraits<ST    >    SCT;
	typedef Belos::MultiVecTraits<ST,MV>     MVT;
	typedef Belos::OperatorTraits<ST,MV,OP>  OPT;
	typedef SCT::magnitudeType                MT;

	// Want to solve the system Ax=b

	// Set problem size
	OT m = 5;
	OT n = 4;

	// 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(domSp);
	Teuchos::RCP<V>  b = Thyra::createMember(rngSp);

	// Initialise input data

	Thyra::randomize(SCT::zero(),SCT::one(),A.ptr());

	// Set value for solution, x=1
	Thyra::put_scalar(SCT::one(),x.ptr());
	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...
	// We want to solve the problem b=Ax. What is x?

	bool verbose = true;
	bool debug = false;
	bool proc_verbose = false;
	int frequency = -1;

	int maxiters = 1000;         // maximum number of iterations allowed
	int maxsubspace = 50;       // maximum number of blocks
	int maxrestarts = 150;       // number of restarts allowed
	MT tol = 1.0e-10;            // relative residual tolerance
	MT lambda = 0.0;            // Preconditioner for LSQR


	Teuchos::ParameterList belosList;
	belosList.set( "Maximum Iterations", maxiters );       // Maximum number of iterations allowed
	belosList.set( "Lambda", lambda );                     // Relative convergence tolerance requested

	int verbosity = Belos::Errors + Belos::Warnings;
	if (verbose) {
		verbosity += Belos::TimingDetails + Belos::StatusTestDetails;
		if (frequency > 0)
			belosList.set( "Output Frequency", frequency );
	}

	if (debug) {
		verbosity += Belos::Debug;
	}

	belosList.set( "Verbosity", verbosity );


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

	Belos::LinearProblem<ST,MV,OP> problem( A, x, b );
	bool set = problem.setProblem();

	if (set == false) {
		std::cout << std::endl << "ERROR:  Belos::LinearProblem failed to set up correctly!" << std::endl;
		return EXIT_FAILURE;
	}



	// Create solver manager.
	Teuchos::RCP< Belos::SolverManager<ST,MV,OP> > newSolver =
			Teuchos::rcp(
					new Belos::LSQRSolMgr<ST,MV,OP>
	( Teuchos::rcp(&problem,false), Teuchos::rcp(&belosList,false) )
			);

	// Solve
	Belos::ReturnType ret = newSolver->solve();
	if (ret!=Belos::Converged) {
		std::cout << std::endl << "ERROR:  Belos did not converge!" << std::endl;
		return EXIT_FAILURE;
	}

	std::cout << std::endl << "SUCCESS:  Belos converged!" << std::endl;
	x->describe(*out,Teuchos::VERB_EXTREME);

	return EXIT_SUCCESS;
}