[Trilinos-Users] Smallest Eigenvalues
Pate Motter
pate.motter at colorado.edu
Wed Mar 25 22:04:53 MDT 2015
I have no preference to using Tpetra, just all of the documentation
indicates it to be the successor to Epetra and so it seemed logical to
choose it if starting a new project.
The error was:
"error: object of abstract class type "Tpetra::Operator<ST={double},
LO={int}, GO={int64_t={long}}, NT>" is not allowed"
for the line: "RCP<op_type> belosPrec = rcp(new op_type(prec.get() ) );"
I was trying to decipher which type of object it was expecting within the
template argument.
Thanks,
Pate
On Wed, Mar 25, 2015 at 8:53 PM, Alicia Klinvex <aklinvex at purdue.edu> wrote:
> Responses are below
>
> - Alicia
>
> On Wed, Mar 25, 2015 at 10:18 PM, Pate Motter <pate.motter at colorado.edu>
> wrote:
>
>> I doubt that the problem is you being unclear, just me being unfamiliar
>> with the inner workings. Especially when porting the documented examples
>> from Epetra to Tpetra.
>>
>> I attempted to make the Tpetra version of the code you provided, and am
>> having troubles with creating the Belos operator. Specifically at the
>> portion documented as: "Create the Belos preconditioned operator from
>> the Ifpack preconditioner NOTE: This is necessary because Belos expects
>> an operator to apply the preconditioner with Apply() NOT ApplyInverse()."
>>
>
> >>> Oh that. Ifpack preconditioners are applied via a function called
> ApplyInverse, so you have to do some special things to get Belos to apply
> them correctly. (Since you seem to be using Ifpack2, I won't get into what
> those things are. If I misunderstood, and you wanted to use Ifpack, let me
> know.) Ifpack2 preconditioners are applied via a function called apply, so
> you don't have to create that strange preconditioned operator.
>
>>
>> My attempt is added below, but it won't compile when attempting to create
>> belosPrec or BelosOp. I can't figure out what types it is wanting so that
>> it plays nicely.
>>
>
> >>> Please send the error the compiler gives you, and I'll tell you how to
> read it.
>
>>
>> Should I just go ahead and follow the trend of the documentation/examples
>> and just revert back to only using Epetra?
>>
>
> >>> It will probably be easier to just take the Epetra example I sent you
> and modify it to solve your problem, but if you desperately want to use
> Tpetra, it's definitely possible to do so. Is there a Tpetra functionality
> you need that is not available in Epetra?
>
>>
>> -Pate
>>
>> RCP<MV> calcSmallestEigenValues(const RCP<const MAT> &A) {
>> typedef Tpetra::Operator<ST, LO, GO, NT> op_type;
>> typedef Ifpack2::Preconditioner <ST, LO, GO, NT> prec_type;
>> Ifpack2::Factory factory;
>> Teuchos::ParameterList plist;
>>
>> // Create preconditioner
>> std::string precondType = "ICT";
>> RCP<prec_type> prec = factory.create(precondType, A);
>> prec->setParameters(plist);
>> prec->initialize();
>> prec->compute();
>>
>> // Set up Belos operator
>> RCP<Belos::LinearProblem<ST, MV, OP> > My_LP =
>> rcp( new Belos::LinearProblem<ST, MV, OP>() );
>> My_LP->setOperator(A);
>> RCP<Belos::Operator> belosPrec = rcp(new Belos::Operator( prec.get() ) );
>> // Compiler error here
>> My_LP->setLeftPrec(belosPrec);
>> RCP<Teuchos::ParameterList> My_List = rcp(new
>> Teuchos::ParameterList());
>> RCP<Belos::Operator> BelosOp = rcp (new Belos::Operator(My_LP,
>> My_List));
>> }
>>
>>
>>
>> On Wed, Mar 25, 2015 at 11:43 AM, Alicia Klinvex <aklinvex at purdue.edu>
>> wrote:
>>
>>> Hello Pate,
>>>
>>> Have you gotten it to work on any matrices? My advice is to try it with
>>> a small matrix first to see if your driver works at all.
>>>
>>> I'm afraid I can't tell you whether your driver is correct based on the
>>> code snippet you provided. I don't see the part where you create the
>>> Anasazi eigensolver. As far as I can tell, what you are doing is creating
>>> a GMRES solver and solving a linear system. Is the linear solve failing?
>>> When you say it fails, do you mean it says "Belos: # iterations, but did
>>> not finish"?
>>>
>>> This is an example which does what you seem to want:
>>> http://trilinos.org/docs/dev/packages/anasazi/doc/html/BlockKrylovSchur_2BlockKrylovSchurEpetraExGenBelos_8cpp-example.html.
>>> That example shows how to use the Anasazi eigensolver BKS with a Belos
>>> linear solver and Ifpack preconditioner.
>>>
>>> As always, let me know if any of this is unclear, and I'd be happy to
>>> elaborate :-)
>>>
>>> Best wishes,
>>> Alicia
>>>
>>> On Wed, Mar 25, 2015 at 12:58 PM, Pate Motter <pate.motter at colorado.edu>
>>> wrote:
>>>
>>>> Hi again Alicia,
>>>>
>>>> I am having some issues with getting my eigenvalue solver working. It
>>>> compiles, but fails to converge.
>>>>
>>>> I am currently using the default parameters which could be affecting it
>>>> somewhat. Since this will be run on a variety of matrices with no common
>>>> properties, I am not sure how to best code it so that it works as
>>>> generically as possible. This was pieced together from a few different
>>>> Trilinos example codes since there was not one doing the exact mixture of
>>>> Belos/Anasazi/Tpetra. So there may be some slight issues when I combined
>>>> them.
>>>>
>>>> Thanks for the help,
>>>> Pate
>>>>
>>>> RCP<MV> calcEigenValues(const RCP<const MAT> &A, std::string eigenType)
>>>> {
>>>>
>>>> // Preconditioner
>>>> const std::string precType = "ILUT";
>>>> Teuchos::ParameterList plist; // leave empty to use default params
>>>>
>>>> // Setup preconditioner using Ifpack2
>>>> typedef Ifpack2::Preconditioner<ST, LO, GO, NT> prec_type;
>>>> RCP<prec_type> prec;
>>>> Ifpack2::Factory factory;
>>>> prec = factory.create(precType, A);
>>>> prec->setParameters(plist);
>>>> prec->initialize();
>>>> prec->compute();
>>>> // Create vectors
>>>> RCP<OP> M = prec;
>>>> RCP<MV> X = rcp (new MV (A->getDomainMap (), 1) );
>>>> RCP<MV> B = rcp (new MV (A->getRangeMap(), 1) );
>>>> B->randomize();
>>>>
>>>> // Solve using Belos
>>>> RCP<Teuchos::ParameterList> solverParams;
>>>> Belos::SolverFactory<ST, MV, OP> solverFactory;
>>>> RCP<Belos::SolverManager<ST, MV, OP> > solver =
>>>> solverFactory.create("GMRES", solverParams);
>>>> typedef Belos::LinearProblem<ST, MV, OP> problem_type;
>>>> RCP<problem_type> problem =
>>>> rcp (new problem_type (A, X, B) );
>>>> problem->setRightPrec(M);
>>>> problem->setProblem();
>>>> solver->setProblem(problem);
>>>> Belos::ReturnType result = solver->solve();
>>>> const int numIters = solver->getNumIters();
>>>> if (result == Belos::Converged) {
>>>> *fos << "Belos: " << numIters << " iteration(s) to finish" <<
>>>> std::endl;
>>>> } else {
>>>> *fos << "Belos: " << numIters << " iteration(s), but did not
>>>> finish" << std::endl;
>>>> }
>>>>
>>>> On 3/17/2015 5:33:01 PM, Alicia Klinvex <aklinvex at purdue.edu> wrote:
>>>> If your matrix has the potential to be indefinite, here are some
>>>> additional options (although BKS will certainly work too):
>>>>
>>>> 1. Use TraceMin-Davidson, and ask it to compute the harmonic Ritz
>>>> values. This is brand new, so if you try this and see any wonky behavior,
>>>> let me know.
>>>> 2. Use spectrum folding. This method has the potential to be very
>>>> slow (since you are in essence squaring the condition number of your
>>>> matrix), but it will allow you to use any eigensolver you wish. Example:
>>>> http://trilinos.org/docs/dev/packages/anasazi/doc/html/LOBPCGCustomStatusTest_8cpp-example.html
>>>>
>>>> If you decide to stick with BKS and have trouble getting your driver to
>>>> work, let us know :-)
>>>>
>>>> - Alicia
>>>>
>>>> On Tue, Mar 17, 2015 at 6:51 PM, Pate Motter <pate.motter at colorado.edu>
>>>> wrote:
>>>>
>>>>> Hi Alicia,
>>>>>
>>>>> Thank you for the information. This is not guaranteed to be an SPD
>>>>> matrix so the first option is probably the best. I remember seeing an
>>>>> example going through that process, but could not get it to work either. I
>>>>> will retry that though as it could have been an error elsewhere.
>>>>>
>>>>> -Pate
>>>>>
>>>>> On 3/17/2015 4:47:03 PM, Alicia Klinvex <aklinvex at purdue.edu> wrote:
>>>>> Hello Pate,
>>>>>
>>>>> If you want to find the smallest eigenvalues, I'd recommend one of the
>>>>> following options:
>>>>> 1. Using BKS in shift-and-invert mode. There are several examples
>>>>> demonstrating how to do this on the Anasazi page. If you want to go this
>>>>> route, you will have to create a linear solver (either direct or iterative)
>>>>> and pass it to the BKS solver manager.
>>>>> 2. Pick a different eigensolver. Is your matrix symmetric positive
>>>>> definite? If so, you can use LOBPCG, RTR, or TraceMin-Davidson. If your
>>>>> matrix is not SPD, your choices are more limited, and we can discuss your
>>>>> options further.
>>>>>
>>>>> Let me know which route you want to go, and I'll try to help you with
>>>>> it :-)
>>>>>
>>>>> Best wishes,
>>>>> Alicia
>>>>>
>>>>> On Tue, Mar 17, 2015 at 6:16 PM, Pate Motter <pate.motter at colorado.edu
>>>>> > wrote:
>>>>>
>>>>>> Hi,
>>>>>>
>>>>>> I am having difficulties converging when finding the smallest
>>>>>> magnitude eigenvalues in a sparse, square, real matrix. I am trying to use
>>>>>> a combination of Tpetra and Anasazi to accomplish this based on a few of
>>>>>> the examples I was able to find online. It currently works for finding
>>>>>> largest magnitude.
>>>>>>
>>>>>> Thanks,
>>>>>> Pate Motter
>>>>>>
>>>>>> RCP<MV> calcEigenValues(const RCP<MAT> &A, std::string eigenType) {
>>>>>> // Get norm
>>>>>> ST mat_norm = A->getFrobeniusNorm();
>>>>>>
>>>>>> // Eigensolver parameters
>>>>>> int nev = 1;
>>>>>> int blockSize = 1;
>>>>>> int numBlocks = 10*nev / blockSize;
>>>>>> ST tol = 1e-6;
>>>>>>
>>>>>> // Create parameters to pass to the solver
>>>>>> Teuchos::ParameterList MyPL;
>>>>>> MyPL.set("Block Size", blockSize );
>>>>>> MyPL.set("Convergence Tolerance", tol);
>>>>>> MyPL.set("Num Blocks", numBlocks);
>>>>>>
>>>>>> // Default to largest magnitude
>>>>>> if (eigenType.compare("SM") == 0) {
>>>>>> MyPL.set("Which", "SM");
>>>>>> } else if (eigenType.compare("SR") == 0) {
>>>>>> MyPL.set("Which", "SR");
>>>>>> } else if (eigenType.compare("LR") == 0) {
>>>>>> MyPL.set("Which", "LR");
>>>>>> } else {
>>>>>> MyPL.set("Which", "LM");
>>>>>> }
>>>>>>
>>>>>> // Create multivector for a initial vector to start the solver
>>>>>> RCP<MV> ivec = rcp (new MV(A->getRowMap(), blockSize));
>>>>>> MVT::MvRandom(*ivec);
>>>>>>
>>>>>> // Create eigenproblem
>>>>>> RCP<Anasazi::BasicEigenproblem<ST, MV, OP> > MyProblem =
>>>>>> Teuchos::rcp(new Anasazi::BasicEigenproblem<ST, MV, OP>(A, ivec));
>>>>>>
>>>>>> MyProblem->setHermitian(false);
>>>>>> MyProblem->setNEV(nev);
>>>>>> MyProblem->setProblem();
>>>>>>
>>>>>> Anasazi::BlockKrylovSchurSolMgr<ST, MV, OP> MySolverMgr(MyProblem,
>>>>>> MyPL);
>>>>>>
>>>>>> // Solve the problem
>>>>>> Anasazi::ReturnType returnCode = MySolverMgr.solve();
>>>>>> if (returnCode != Anasazi::Converged) {
>>>>>> *fos << "unconverged" << ", ";
>>>>>> }
>>>>>>
>>>>>> // Get the results
>>>>>> Anasazi::Eigensolution<ST, MV> sol = MyProblem->getSolution();
>>>>>> std::vector<Anasazi::Value<ST> > evals = sol.Evals;
>>>>>> RCP<MV> evecs = sol.Evecs;
>>>>>> int numev = sol.numVecs;
>>>>>>
>>>>>> // Compute residual
>>>>>> if (numev > 0) {
>>>>>> Teuchos::SerialDenseMatrix<int, ST> T(numev,numev);
>>>>>> for (int i = 0; i < numev; i++) {
>>>>>> T(i,i) = evals[i].realpart;
>>>>>> }
>>>>>> std::vector<ST> normR(sol.numVecs);
>>>>>> MV Kvec(A->getRowMap(), MVT::GetNumberVecs(*evecs));
>>>>>> OPT::Apply(*A, *evecs, Kvec);
>>>>>> MVT::MvTimesMatAddMv(-1.0, *evecs, T, 1.0, Kvec);
>>>>>> MVT::MvNorm(Kvec, normR);
>>>>>> for (int i=0; i<numev; i++) {
>>>>>> *fos << evals[i].realpart << ", " << normR[i]/mat_norm << ", ";
>>>>>> }
>>>>>> }
>>>>>> return evecs;
>>>>>> }
>>>>>>
>>>>>>
>>>>>> _______________________________________________
>>>>>> Trilinos-Users mailing list
>>>>>> Trilinos-Users at software.sandia.gov
>>>>>> https://software.sandia.gov/mailman/listinfo/trilinos-users
>>>>>>
>>>>>>
>>>>>
>>>>
>>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://software.sandia.gov/pipermail/trilinos-users/attachments/20150325/19d7a50d/attachment.html>
More information about the Trilinos-Users
mailing list