[Trilinos-Users] [EXTERNAL] slow anasazi in NOX::Epetra

[Veltz Romain]

Andy,

My mass matrix is the identity matrix. I hope it is not wrong (see below).

If I use the MatrixFree Jacobian, I get the same results. They are not wrong per se as if the function transformEigenvalue were not working.

If I use the Shift-Invert settings, I have the following output which insure that the class is called

WARNING: LOCA::AnasaziOperator::ShiftInvert::apply() - 
Return Type = NotConverged

Hence, I have no clue what can be wrong. You suggested that my function ProblemInterface::applyShiftedMatrix is buggy and/or not called but I think it is called because of what is reported above.

Thank you again for your suggestions,

Romain

PS: To help, I copied the part of the code that might relevant. Here are the definitions for the interface/group

Teuchos::RCP<LOCA::Epetra::Interface::TimeDependentMatrixFree> iReq = interface;
Teuchos::RCP<FFTJac> MFJac = Teuchos::rcp(new FFTJac(nlPrintParams, interface, *noxSoln, Problem));
Teuchos::RCP<LOCA::Epetra::Interface::TimeDependentMatrixFree> iReq = interface;
Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = MFJac;
// Create the linear systems
Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linsys =
			Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(nlPrintParams, lsParams,
					iReq, iJac, MFJac, *soln));
Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> shiftedLinsys =
			Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(nlPrintParams, lsParams,
					iReq, iJac, MFJac, *soln));
// access Jacobian from interface
interface->setJacobianOperator(MFJac);

// Create the Group
	Teuchos::RCP<LOCA::MultiContinuation::AbstractGroup> grp =
			Teuchos::rcp(new LOCA::Epetra::Group(
					globalData,
					nlPrintParams,
					iReq, 			// interface
					locaSoln,		// initial guess
					linsys,			// linear system
					shiftedLinsys,	// shifted linear system
					pVector));		// parameters vector

And, in my interface, I have
bool setJacobianOperator(Teuchos::RCP<Epetra_Operator> Jac)
	{
		Jac_= Jac;
	}

bool ProblemInterface::applyShiftedMatrix(double alpha, double beta,
		const NOX::Epetra::Vector& input,
		NOX::Epetra::Vector& result) const
{
	// compute result = alpha * J * input + beta * input
	Epetra_Vector X(input.getEpetraVector());
	Epetra_Vector JX(result.getEpetraVector());
	if (alpha!=0)
		assert(Jac_->Apply(X,JX)==0);
	result.update(alpha,JX,beta,X,0.0);
	return true;
}

Romain

On Dec 11, 2013, at 5:40 PM, "Salinger, Andrew" <agsalin at sandia.gov> wrote:

> 
> Veltz,
> 
> I have two ideas:
> 
> The identity matrix, while seemingly simple, has all eigenvalues repeated. This can cause problems for the iterative eigensolver. Can you try again with diagonal entries 1, 1.1, 1.2, 1.3,..
> 
> There might be a problem with your shifted/mass matrix calculation, since that is not used in the Jacobian Inverse method, but is used in the ones with a wrong results. You can
> try some tests of applying it agains simple vectors to make sure it is correct.
> 
> Andy
> 
> From: Veltz Romain <romain.veltz at inria.fr>
> Date: Wednesday, December 11, 2013 5:59 AM
> To: Eric Phipps <etphipp at sandia.gov>
> Cc: "trilinos-users at software.sandia.gov" <trilinos-users at software.sandia.gov>
> Subject: Re: [Trilinos-Users] [EXTERNAL] slow anasazi in NOX::Epetra
> 
> I thought I would be done with issues but I found the following output from Anasazi rather puzzling.
> Note that, for some reasons, I am using trillions-10.8.4 (updating to the current version give me some trouble under macosx).
> 
> I set my continuation code in a configuration where the Jacobian is the operator: -Id
> 
> Using the LOCA with Epetra - Anasazi with:
> aList.set("Sorting Order", "LM");
> aList.set("Operator","Jacobian Inverse");
> 
> the following EV are returned
> 
> Untransformed eigenvalues (since the operator was Jacobian Inverse)
> Eigenvalue 0 : -1.000e+00  -0.000e+00 i    :  RQresid 3.331e-16  0.000e+00 i
> Eigenvalue 1 : -1.000e+00  -0.000e+00 i    :  RQresid 6.439e-15  0.000e+00 i
> Eigenvalue 2 : -1.000e+00  -0.000e+00 i    :  RQresid 7.327e-15  0.000e+00 i
> Eigenvalue 3 : -1.000e+00  -0.000e+00 i    :  RQresid 5.995e-15  0.000e+00 i
> 
> 
> Using the LOCA with Epetra - Anasazi with:
> aList.set("Sorting Order", "LM");
> aList.set("Operator","Shift-Invert");
> aList.set("Shift",0.001);
> 
> the following ev are returned (The RQresid is large too, equal to the shift…):
> 
> Untransformed eigenvalues (since the operator was Shift-Invert)
> Eigenvalue 0 : -9.990e-01  -0.000e+00 i    :  RQresid 1.000e-03  0.000e+00 i
> Eigenvalue 1 : -9.990e-01  -0.000e+00 i    :  RQresid 1.000e-03  0.000e+00 i
> 
> Note that I have checked that the solver effectively enter the following function
> ->LOCA::AnasaziOperator::ShiftInvert::apply()
> 
> I don't understand these results… I thought the entry  in
> anasaziOp->transformEigenvalue((*evals_r)[i], (*evals_i)[i]);
> from LOCA::Eigensolver::AnasaziStrategy::computeEigenvalues would do the job…
> 
> Please note that I have the same issues with the Cayley transform.
> 
> Thank you for your help and advices,
> 
> Romain
> 
> 
> On Dec 10, 2013, at 5:33 PM, "Phipps, Eric T" <etphipp at sandia.gov> wrote:
> 
>> I don't understand what you mean be "I don't have access to my Jacobian operator".  You need to be able to implement y = (a*J + b*M)*x for a give (multi-) vector x and scalars a and b, where J is the Jacobian operator and M the mass matrix.  If you have a matrix-free operator to implement J*x, you can use the same approach to implement this, e.g., via y = a*J*x + b*M*x.
>> 
>> -Eric
>> 
>> On Dec 10, 2013, at 5:03 AM, "Veltz Romain" <romain.veltz at inria.fr> wrote:
> 

-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20131211/059ec0f7/attachment.html