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

Phipps, Eric T etphipp at sandia.gov
Tue Dec 10 09:33:35 MST 2013


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<mailto:romain.veltz at inria.fr>> wrote:

Hi Andrew,

I had anticipated that 'mistake' of yours. Hence, I still have my issue.

If I inherit TimeDependentMatrixFree in my interface, I have to implement in my interface applyShiftedMatrix which I cannot because I don't have access to my jacobian operator in my interface. Should I make this function  applyShiftedMatrix  return false in the interface?

I understand that I should create a second NOX::Epetra::LinearSystemAztecOO for my shifted system. Hence, I should create something like

Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> shiftedLinsys =
Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(nlPrintParams, lsParams,iReq, iJac, MFShift, *soln));

where MFShift is an Epetra_Operator. I don't understand MFShift  what should inherit. Indeed, I need the shift and the jacobian to be able to do the Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) function.

--> Basically, I don't know what kind of class (i.e. Epetra_Operator) I should give to NOX::Epetra::LinearSystemAztecOO
--> I also don't know how to access the shift value from this class which is needed to implement apply(…)


I hope I made my issues more clear,

Thank you for your help,

Romain

On Dec 9, 2013, at 5:15 PM, Andy Salinger <agsalin at sandia.gov<mailto:agsalin at sandia.gov>> wrote:


Hi Romain,

I looked a little more into it and see that you probably want to
inherit from LOCA_Epetra_Interface_TimeDependentMatrixFree
and implement applyShiftedMatrix. My previous email was not
for the matrix-free approach. This uses a different constructor
to LOCA_Epetra_Group that takes this interface. It looks like
you need to create a second NOX::Epetra::LinearSystemAztecOO
for this constructor for the shifted solve. I see some code that
constructs both linSys and shiftedLinSys on consecutive lines
with only the operator being different.
(packages/piro/src/LOCA_Epetra_LOCASolver.cpp)


Hope this helps,
Andy


On 12/08/2013 12:48 AM, Veltz Romain wrote:
Hi Andrew,

I am implementing the interface that you suggested and I am a bit stuck because I can't access the jacobian class from my interface. Maybe you can give me a hint.

I have a class FFTJac which is derived from Epetra_Operator and NOX::Interface::Jacobian which is used to define a Jacobian Operator (Matrix Free) from which I build a linear system NOX::Epetra::LinearSystemAztecOO that is passed to the group. In order to implement applyShiftedMatrix  in my interface, I should be able to access FFTJac::apply (after computeJacobian has been called….) but I don't know how to do that since I can't access the group from my interface.

Romain

On Dec 3, 2013, at 6:11 PM, "Salinger, Andrew" <agsalin at sandia.gov<mailto:agsalin at sandia.gov>> wrote:


Hi Veltz,

Yes, for Shift-Invert you need to provide an interface to apply a mass matrix.
There is the LOCA_Epetra_Interface_TimeDependent that has a computeShiftedMatrix
call that is used for both the mass matrix and the shifted matrix.

Regarding your settings: here are my recommendations
Block Size = 1 (unless you expect redundant eigenvalues due to symmetries)
Num Blocks = 200  (Or some number big enough)
Maximum Iterations = 200 (skip restart for now)
Step Size = 1 (Check tolerance every iteration)
  perhaps loosen the tolerance as well.

Andy

From: Veltz Romain <romain.veltz at inria.fr<mailto:romain.veltz at inria.fr>>
Date: Tuesday, December 3, 2013 3:34 AM
To: "Day, David" <dmday at sandia.gov<mailto:dmday at sandia.gov>>
Cc: "trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>" <trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>>
Subject: Re: [Trilinos-Users] [EXTERNAL] slow anasazi in NOX::Epetra

David,

I tried but the speed is not that great.

I also tried to use a shift-transform but I get an error with this

aList.set("Method", "Anasazi");
aList.set("Sorting Order", "LM");  // Utiliser LM pour avoir les points de bifurcation
aList.set("Block Size", 4);//1        // Size of blocks, ie Krylov subspace
aList.set("Num Blocks", 30);       // Size of Arnoldi factorization: = 3 nev
aList.set("Num Eigenvalues", 6);   // Number of eigenvalues
aList.set("Convergence Tolerance", 1.0e-12);          // Tolerance
aList.set("Step Size", 20);         // How often to check convergence
aList.set("Maximum Restarts",5);  // Maximum number of restarts
aList.set("Maximum Iterations",500);   // Maximum number of restarts
aList.set("Operator","Shift-Invert");
aList.set("Shift",0.5);

The error is the following:

================================================================
Anasazi Eigensolver starting with block size 4

************************
ERROR: LOCA::AnasaziOperator::ShiftInvert::apply()

Return Type = BadDependency
************************
[peguy:48825] *** Process received signal ***
[peguy:48825] Signal: Segmentation fault: 11 (11)
[peguy:48825] Signal code: Address not mapped (1)
[peguy:48825] Failing at address: 0x15
[peguy:48825] [ 0] 2   libsystem_c.dylib                   0x00007fff916f390a _sigtramp + 26
[peguy:48825] [ 1] 3   ???                                 0x00007ff680cfc080 0x0 + 140696699781248
[peguy:48825] [ 2] 4   NContDoG.exe                        0x0000000104c60b72 _ZN7Anasazi22BlockKrylovSchurSolMgrIdN3NOX8Abstract11MultiVectorEN4LOCA15AnasaziOperator16AbstractStrategyEE5solveEv + 22696
[peguy:48825] [ 3] 5   NContDoG.exe                        0x0000000104c3d6c6 _ZN4LOCA11Eigensolver15AnasaziStrategy18computeEigenvaluesERN3NOX8Abstract5GroupERN7Teuchos3RCPISt6vectorIdSaIdEEEESC_RNS7_INS3_11MultiVectorEEESF_ + 1936
[peguy:48825] [ 4] 6   NContDoG.exe                        0x0000000104c1e343 _ZN4LOCA7Stepper5startEv + 1391
[peguy:48825] [ 5] 7   NContDoG.exe                        0x0000000104b2293b _ZN4LOCA8Abstract8Iterator3runEv + 21
[peguy:48825] [ 6] 8   NContDoG.exe                        0x00000001048f5c8a main + 22426
[peguy:48825] [ 7] 9   libdyld.dylib                       0x00007fff9113e7e1 start + 0
[peguy:48825] *** End of error message ***

Does this mean that I need to create a class with a specific Mass matrix?


Thank you for your help,

Romain




On Nov 29, 2013, at 9:13 PM, "Day, David" <dmday at sandia.gov<mailto:dmday at sandia.gov>> wrote:

Vomain,
I too am sorry for my ignorance.  I am not sure how to diagnose solver issues through a parameter list either. Today is a US holiday.  We night need to wait a little for help.
--David Day

From: Veltz Romain <romain.veltz at inria.fr<mailto:romain.veltz at inria.fr>>
Date: Thursday, November 28, 2013 5:24 AM
To: David Day <dmday at sandia.gov<mailto:dmday at sandia.gov>>
Cc: "Thornquist, Heidi K" <hkthorn at sandia.gov<mailto:hkthorn at sandia.gov>>
Subject: Re: [EXTERNAL] [Trilinos-Users] slow anasazi in NOX::Epetra

David,

I am sorry for my ignorance but I don't know how to look at the number of times that GMRES restart.

Teuchos::ParameterList& nlPrintParams = nlParams.sublist("Printing");
nlPrintParams.set("Output Precision",3);
nlPrintParams.set("Output Processor",0);
nlPrintParams.set("Output Information",0 +
//NOX::Utils::OuterIteration +
//NOX::Utils::OuterIterationStatusTest +
//NOX::Utils::InnerIteration +
//NOX::Utils::Details +
//NOX::Utils::Warning+
NOX::Utils::StepperIteration +
//NOX::Utils::StepperDetails +
//NOX::Utils::StepperParameters +
NOX::Utils::LinearSolverDetails +
0
);

Is it the way?

Thank you again for your help to you and Heidi

Romain



On Nov 27, 2013, at 5:27 PM, "Day, David" <dmday at sandia.gov<mailto:dmday at sandia.gov>> wrote:

Romain,
There is no simple answer.   One issue that is mostly independent is simply the efficiency of the linear solve.   Here look for the number of times that GMRES must restart.   The other issues relate more directly to Anasazi.   There is a natural conflict between continuation methods and iterative linear solvers.   The continuation method drives the Jacobian matrix to be singular,  in many cases.   Loca is carefully designed to mitigate this issue, but I guess that this is why you are solving an eigenvalue problem.

You might send this to the Anasazi mailing list for more discussion.  Block Size one is a good guess, unless for some reason your problem has multiple eigenvalues.  You might try a much larger value of Num Blocks.   You might also try to experiment with the Convergence Tolerance.   There is also a trade off between the linear solver and Anasazi;  a more accurate linear solves (smaller "Tolerance") increases the cost per linear solve and can decrease the number of Anasazi iterations.

-David Day
From: Veltz Romain <romain.veltz at inria.fr<mailto:romain.veltz at inria.fr>>
Date: Wednesday, November 27, 2013 9:02 AM
To: "trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>" <trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>>
Subject: [EXTERNAL] [Trilinos-Users] slow anasazi in NOX::Epetra

Hello All,

I am doing a continuation in the NOX:Epetra interface with a Epetra_OPerator for the Jacobian. As a bi-product, I am computing the eigenvalues of the Jacobian at each step using the following piece of code.


Teuchos::ParameterList& aList = locaStepperList.sublist("Eigensolver");
aList.set("Method", "Anasazi");
aList.set("Sorting Order", "LM");  //
aList.set("Block Size", 1);//1        // Size of blocks, ie Krylov subspace
aList.set("Num Blocks", 20);//15       // Size of Arnoldi factorization: = 3 nev
aList.set("Num Eigenvalues", 5);   // Number of eigenvalues
aList.set("Convergence Tolerance", 1.0e-7);          // Tolerance
aList.set("Step Size", 1);         // How often to check convergence
aList.set("Maximum Restarts",20);  // Maximum number of restarts

Fore the linear solver, the options are here:

Teuchos::ParameterList& lsParams = newParams.sublist("Linear Solver");
lsParams.set("Aztec Solver", "GMRES");
lsParams.set("Size of Krylov Subspace",100);//default = 300
lsParams.set("Max Iterations", 100);// etait 200
lsParams.set("Tolerance", 1e-8);
lsParams.set("Output Frequency", 1);//etait 1
lsParams.set("Scaling", "None");
lsParams.set("Preconditioner Operator", "Use Jacobian");// default

I would like to point that the problem I am trying to solve is

F(V,lambda) = -V + convolution2d(J, N(lambda*V))

where N is a nonlinearity.

I have two questions:

- how do the parameters Block size and Num blocks affect the computation of the eigenvalues?

I have done the same computation in matlab with the function eigs and it is at least twice as fast.

- how can I improve the speed of the computation of my eigenvalues in Anasazi. I think this is linked to the first question.

Thank you for your help,

Best,

Romain






_______________________________________________
Trilinos-Users mailing list
Trilinos-Users at software.sandia.gov<mailto: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/20131210/0e66d41a/attachment-0001.html 


More information about the Trilinos-Users mailing list