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

Veltz Romain romain.veltz at inria.fr
Fri Dec 13 00:02:36 MST 2013 

I don't really get what you are saying.

Are you saying that shiftedlinsys should take another Epetra_Operator instead of MFJac?

If yes, I don't know how to make it the shifted operator (is it just the mass matrix, i.e. identity for me?)? 

How this new operator would be related to the interface function bool ProblemInterface::applyShiftedMatrix ? 

On Dec 13, 2013, at 12:00 AM, "Salinger, Andrew" <agsalin at sandia.gov> wrote:

> 
> It looks like your shifted linear system is with the Jacobian operator. I think
> it needs to be with the shifted operator. 
> 
> 
> From: Veltz Romain <romain.veltz at inria.fr>
> Date: Wednesday, December 11, 2013 2:38 PM
> To: Andy Salinger <agsalin 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
> 
> 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(newLOCA::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);
> returntrue;
> }
> 
> 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/20131213/3efda049/attachment.html

More information about the Trilinos-Users mailing list