[Trilinos-Users] Thyra::initializeAndReuseOp

[Heroux, Mike]

Laurent,

Yes, these capabilities are available in Trilinos in general.  I am not the lead expert on accessing them via Thyra, others can better address that, but all the underlying libraries support this feature.  In addition, most of the packages support the concept of a matrix graph as an independent object, so not only can you reuse the graph from one time step to another, but if you have multiple matrices at a given time step that have the same nonzero structure, these matrices can share the same graph.

Mike

From: Laurent Griesser <laurent.griesser at gmail.com<mailto:laurent.griesser at gmail.com>>
Date: Sunday, January 18, 2015 at 6:55 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: [Trilinos-Users] Thyra::initializeAndReuseOp

Hi,

I am new to Trilinos and still not used to the use of linear solvers. I want to solve a parabolic differential equation. Each time iteration naturally consists in solving a linear system. From one time iteration to the other, I would like to reuse my old linear system as much as I can. The matrix of my linear system is penta-diagonal and this structure doesn't change from one iteration to the other. Only the matrix coefficients do change.

I think I can do that with Thyra::initializeAndReuseOp but I am not sure what this method really does. Does it update my ILU factorization? How efficient is that compared to the uninitializeOp/initializeOp strategy?

I have partially read the documentation under

http://trilinos.org/docs/dev/packages/thyra/doc/html/classThyra_1_1LinearOpWithSolveFactoryBase.html

They distinguish between

"Reuse of factorizations for small changes in the forward operator"

and

"Updating the linear solver for major changes in the structure of the forward operator"

What does "small changes" mean? What is the efficiency penalty when I proceed with the "major changes" method?

With PETSc for example, I don't need to think about this. Usually, I can write something along those lines:

// linear system setup phase
KSPCreate(PETSC_COMM_SELF, &m_Solver);
KSPSetType(m_Solver, KSPGMRES);
KSPGMRESSetRestart(m_Solver, 1500);
KSPSetTolerances(m_Solver, 1e-5, PETSC_DEFAULT, 1e10, 1500);
KSPGMRESSetOrthogonalization(m_Solver, KSPGMRESModifiedGramSchmidtOrthogonalization);
PC pc; KSPGetPC(m_Solver, &pc); PCSetType(pc, PCILU);

/*
 * solve parabolic equation
 */
while(t < tf) {
  // update of matrix/rhs
  ...

  // solve linear system
  KSPSetOperators(m_Solver, m_A, m_A, SAME_NONZERO_PATTERN);
  KSPSolve(m_Solver, m_b, m_X);

  // post-processing
  ...
  // time update
  tf += dt;
}

And here, the fact that I mention "SAME_NON_ZERO_PATTERN" tells my solver that I don't want a complete refactorization. This is exactly what I want from trilinos. How do I achieve this?

Thanks for your advice,
Laurent