[Trilinos-Users] bug in AnasaziBasicEigenProblem

[Thornquist, Heidi K]

Hi Simone,

The Anasazi::BasicEigenproblem is a basic implementation of the abstract Anasazi::Eigenproblem class that can be used with all three eigensolvers (block Krylov-Schur, block Davidson, and LOBPCG).  Op and AOp may be different in the case of block Krylov-Schur where AOp may be defined from the original eigenproblem definition:

AOp*V = Lambda*MOp*V

However, a spectral transformation may be required to compute the eigenvalues of interest.  Thus the real eigenproblem that block Krylov-Schur is solving is:

Op*V = Lambda*V

where Op is some spectral transformation.  Thus, the Anasazi::BasicEigenproblem allows you to define both the original eigenproblem and the one you would like to be solved.  In the case of block Davidson and LOBPCG, the problem being solved is:

Op*V = Lambda*MOp*V.

Thus, you want to make sure that Op is equivalent to the AOp if you are looking to solve:

AOp*V = Lambda*MOp*V

The easiest way to do this using Anasazi::BasicEigenproblem is to SetOperator().  I would also like to note that Anasazi::BasicEigenproblem is just an example of how one may implement the Anasazi::Eigenproblem interface.  It's not perfect, but just meant to be flexible enough to satisfy the current Anasazi eigensolvers.  Please feel free to take that implementation and make it your own.

Let me know if you have any other problems!

Thanks,  Heidi

 


 

________________________________

From: trilinos-users-bounces at software.sandia.gov on behalf of Simone Deparis
Sent: Fri 3/17/2006 1:21 PM
To: trilinos-users at software.sandia.gov
Subject: [Trilinos-Users] bug in AnasaziBasicEigenProblem



The shortest way to describe what happens is to give an example:

I have
MyProblem =
       Teuchos::rcp( new Anasazi::BasicEigenproblem<double, MV, OP>(A,
M, ivec) );
Anasazi::LOBPCG<double, MV, OP> MySolver(MyProblem, MySort, MyOM, MyPL);

MyProblem->SetProblem();

I solve my Eigenproblem once and then replace the operator A and M:
MyProblem->SetA (B);
MyProblem->SetM (N);
MyProblem->SetProblem();

Now I solve again using MySolver, but ...

   template <class ScalarType, class MV, class OP>
   ReturnType BasicEigenproblem<ScalarType, MV, OP>::SetProblem()
   {
     //----------------------------------------------------------------
     // Sanity Checks
     //----------------------------------------------------------------
(...)
     // If there is an A, but no operator, we can set them equal.
     if (_AOp.get() && !_Op.get()) { _Op = _AOp; }
(...)
    }

and in
   template <class ScalarType, class MV, class OP>
   LOBPCG<ScalarType,MV,OP>::LOBPCG(const
Teuchos::RefCountPtr<Eigenproblem<ScalarType,MV,OP> > &problem,
                                    const
Teuchos::RefCountPtr<SortManager<ScalarType,MV,OP> > &sm,
                                    const
Teuchos::RefCountPtr<OutputManager<ScalarType> > &om,
                                    Teuchos::ParameterList &pl
                                    ):
     _problem(problem),
     _sm(sm),
     _om(om),
     _pl(pl),
     _Op(_problem->GetOperator()), // <<<<<<<<<
     _MOp(_problem->GetM()),
     _Prec(_problem->GetPrec()),

Since in the constructor of MyProblem _Op, not _AOp, as been setted
equal to A, my new defined operator B will never be used!

One way around is to use SetOperator instaed of SetA

Question: Why there are two operators, _AOp and _Op?

What do you think?

Simone



--
_____________________________________________________________________
                                     Simone Deparis
   .~.                               Mech Eng Dept - MIT
   /V\   L   I   N   U   X           77 Mass Ave Room 3-264
  // \\  =================           Cambridge MA 02139
/(   )\                             USA
  ^^-^^
phone :  +1 617 452 3285  mailto:deparis at mit.edu
fax   :  +1 617 258 8559  http://www.mit.edu/~deparis
_____________________________________________________________________



_______________________________________________
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: http://software.sandia.gov/mailman/private/trilinos-users/attachments/20060320/c16168d6/attachment.html