[Trilinos-Users] Smallest Eigenvalues
Alicia Klinvex
aklinvex at purdue.edu
Wed Apr 1 20:25:16 MDT 2015
Hello Pate,
Unfortunately, I don't know that much about Epetra and Ifpack. The great
thing about the Trilinos mailing list though is that somebody who DOES know
about those packages is bound to see your message.
My initial guess would be that the issue is with your &*K. Maybe try
K.get() instead (
http://trilinos.org/docs/dev/packages/teuchos/doc/html/classTeuchos_1_1RCP.html
)?
Best wishes,
Alicia
On Mon, Mar 30, 2015 at 10:29 PM, Pate Motter <pate.motter at colorado.edu>
wrote:
> Hi Alicia,
>
> I have been playing with Epetra some more, but am getting stuck with a
> segfault. I took the example at (
> http://trilinos.org/docs/r11.12/packages/anasazi/doc/html/BlockKrylovSchur_2BlockKrylovSchurEpetraExGenBelos_8cpp-example.html
> ) and tried to change it as little as possible to what I want.
>
> The only thing I have changed so far is the default K and M CrsMatrix
> objects to the following:
> //Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp(
> const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
> //Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp(
> const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );
> Epetra_CrsMatrix *A;
> EpetraExt::MatrixMarketFileToCrsMatrix("LFAT5.mtx", Comm, A);
> Teuchos::RCP<Epetra_CrsMatrix> K(A, false), M(A, false);
>
> The code segfaults when creating the Ifpack factory just a few lines later:
> Teuchos::RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp(
> Factory.Create(PrecType, &*K, OverlapLevel) );
>
> It seems to be an issue with how I am assigning the matrices via raw
> pointers and the RCPs. I wasn't sure if there was a clear way to avoid this
> since reading in the MatrixMarket file has to be passed a raw pointer to
> the Epetra_CrsMatrix.
>
> Any help would be greatly appreciated.
>
> Thanks again,
> Pate
>
>
> On 3/26/2015 9:30:27 PM, Alicia Klinvex <aklinvex at purdue.edu> wrote:
> My advice is to stick to Epetra for now.
>
> You tried to create a new Tpetra::Operator, but Tpetra::Operator is an
> abstract base class...so you can't do that. (I'm a big fan of the
> explanation of abstract base classes here:
> http://www.cplusplus.com/doc/tutorial/polymorphism/) That's what the
> error message is trying to tell you.
>
> - Alicia
>
> On Thu, Mar 26, 2015 at 12:04 AM, Pate Motter <pate.motter at colorado.edu>
> wrote:
>
>> 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/20150401/442b165d/attachment.html>
More information about the Trilinos-Users
mailing list