[Trilinos-Users] Anasazi Block Krylov Schur Trouble With Preconditioners

[Thornquist, Heidi K]

Hi Davood,

If I am understanding you correctly, all you are trying to do is change the preconditioner for the linear solver used in applying the shift-invert operator in BKS.  That shift-invert operator is applying K^{-1}*M using the Anasazi::EpetraGenOp class, where K^{-1} is being approximated by the AztecOO GMRES solver, which you want to precondition by an ILU because the eigenproblem is non-Hermitian.

Given those assumptions, this error indicates that the matrix K, or the preconditioner, is ill-conditioned and the first linear solve was not able to be performed.  This kind of failure is unrecoverable for BKS.  If you don't expect your operator K to be ill-conditioned, or nearly singular, I would suggest making sure they are correct.  This could be double checked by solving K*x = M*b, where b is a random vector, using AztecOO.  If you are looking for very small eigenvalues, then I would suggest using a non-zero shift (sigma) away from the origin, so that you apply the operator (K-sigma*M)^{-1}*M instead of K^{-1}*M.  Chris described this process in detail in a previous email, so I've copied that below just in case you need it.

Either way, the operator K seems to be the problem and the fact that any change of preconditioner results in the same behavior supports this suggestion.

Regards,
Heidi

--
Heidi K. Thornquist
Electrical & Microsystems Modeling
Sandia National Laboratories
P.O. Box 5800, MS 0316
Albuquerque, NM  87185-0316
Office:  505-284-8426
Fax:     505-284-5451





>From Chris 2/18/2008:

When looking for interior eigenvalues with a Arnoldi solver (like Anasazi's BKS), the spectral transformation is typically a shift-invert transformation. Consider the following generalized eigenvalue problem:
   A*x = B*x*lambda,
where A and B are linear operators, x is a non-zero vector, and lambda is a scalar.
Remove sigma*B*x (for scalar sigma, your pre-specified value of interest) from both sides:
  (A - sigma*B)*x = B*x*(lambda - sigma)
Premultiply both sides by the inverse of (A - sigma*B) and 1/(lambda - sigma):
  inv(A - sigma*B)*B x = x 1/(lambda - sigma)
This is a new (standard) eigenvalue problem, where the eigenvalues closest to sigma have become the eigenvalues of largest magnitude, i.e., the eigenvalues which an Arnoldi solve will naturally find:
  Op = inv(A-sigma*B)*B, (the spectral transformation)
 Op*x = x theta,
where
   theta = 1/(lambda-sigma), or
   lambda = 1/theta+sigma   (the back-transformed eigenvalue)
Note also that if A and B are symmetric and B is positive semi-definite, then the operator Op is symmetric with respect to the inner product defined by B. In this scenario, you want to pass both Op and B to the Anasazi::Eigenproblem and specify that the eigenproblem is symmetric/Hermitian, as in the examples cited below.
Note, in the current version of Anasazi, because Anasazi is not aware of the spectral transformation present in Op, it cannot back-transform the computed eigenvalues, i.e., turn the thetas corresponding to the transformed eigenproblem into the lambdas corresponding to the original eigenproblem. You must do this after the solver returns. See the examples cited below.


On 2/25/08 8:51 PM, "Davood Ansari" <david.ansari at gmail.com> wrote:

Hi everyone

Well, I am trying to modify the example with Anasazi BKS to fit into my problem.
The major change I am making is that I use ifpack preconditioner other than that
of the example ( BlockKrylovSchurEpetraExGenAztecOO.cpp ). The issue is
that the code numerically fails when any preconditioner other than that of the original
code ( BlockKrylovSchurEpetraExGenAztecOO.cpp ) is used. The error message is
as:

Aztec status AZ_ill_cond: GMRES hessenberg ill-conditioned
Error! Caught exception in BlockKrylovSchur::iterate() at iteration 1
std::exception
terminate called after throwing an instance of 'Anasazi::OperatorError'
  what():  /opt/Trilinos/Parallel/include/AnasaziEpetraAdapter.hpp:792:

Throw number = 1

Throw test that evaluated to true: Op.Apply( x, y ) != 0

Error in Epetra_Operator::Apply()!

For the above case I had used the following chunk of code to build my preconditioner: (I have also tried
preconditioners form ifpack's factory class )

    Teuchos::RCP<Ifpack_IlukGraph> IlukGraph;
    Teuchos::RCP<Ifpack_CrsRiluk> ILUK;
    // This block examaplifies the use of Ifpack (without using the Factory class) derived from example progrma from Michael A. Heroux
    if (prec_ilu_fill>-1)
    {
        Epetra_Time timer(Comm);
        double elapsed_time = timer.ElapsedTime();
        IlukGraph = Teuchos::rcp( new Ifpack_IlukGraph((*A).Graph(), prec_ilu_fill, prec_ptn_overlap) );
        assert(IlukGraph->ConstructFilledGraph()==0);
        ILUK = Teuchos::rcp( new Ifpack_CrsRiluk(*IlukGraph) );
        ILUK->SetAbsoluteThreshold(prec_abs_tresh);
        ILUK->SetRelativeThreshold(prec_rel_tresh);
        int initerr = ILUK->InitValues(*A);
        if (initerr!=0) cout << Comm << "InitValues error = " << initerr;
        assert(ILUK->Factor()==0);

        elapsed_time = timer.ElapsedTime() - elapsed_time;
        if(rank == 0)
        {
            cout << " done." << std::endl;
            cout << "\tTime to construct " << prec_type << " preconditioner: " << elapsed_time << "sec."<<endl;
        }
    }


The original code, however, uses the following chunk to build the preconditioner:

   Teuchos::RCP<Ifpack_CrsIct> ICT;
    ICT = Teuchos::rcp( new Ifpack_CrsIct(*A, prec_drop_tol, prec_ilut_fill) );
    ICT->SetAbsoluteThreshold(prec_abs_tresh);
    ICT->SetRelativeThreshold(prec_rel_tresh);
    int initerr = ICT->InitValues(*A);
    if (initerr != 0) cout << "\tInitValues error = " << initerr ;
    int info = ICT->Factor();
    assert( info==0 );

    if (rank==0)
    {
        bool transA = false;
        double Cond_Est;
        ICT->Condest(transA, Cond_Est);
        cout << "\tCondition number estimate for this preconditoner = " << Cond_Est << endl;
    }

Any idea why should other preconditioners (or even ICT made through the ifpack factory) should lead to this failure ?

Yours
Davood


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://software.sandia.gov/mailman/private/trilinos-users/attachments/20080226/3272f480/attachment.html