[Trilinos-Users] [EXTERNAL] slow anasazi in NOX::Epetra
Veltz Romain
romain.veltz at inria.fr
Wed Jan 22 09:22:26 MST 2014
I still don't understand my mistake.
Anyone has an idea please?
On Dec 22, 2013, at 4:30 PM, Veltz Romain <romain.veltz at inria.fr> wrote:
> Andrew,
>
> Are you saying that shiftedlinsys should take another Epetra_Operator instead of MFJac?
>
> 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:
>>>
>>
>
> _______________________________________________
> Trilinos-Users mailing list
> Trilinos-Users at software.sandia.gov
> http://software.sandia.gov/mailman/listinfo/trilinos-users
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20140122/62f3550a/attachment.html
More information about the Trilinos-Users
mailing list