[Trilinos-Users] Smallest Eigenvalues
Pate Motter
pate.motter at colorado.edu
Thu Apr 2 14:19:39 MDT 2015
Mike,
I made the change you recommended and it is still giving me the same
results.
As of right now I am trying to tweak the example to just run on any of the
matrices found in the Florida Sparse Matrix Collection and determine their
smallest eigenvalues. For testing purposes right now I am just working on
square, general, positive-definite matrices.
I have run the example code on them and can never get them to converge even
though they converge when using different non-Trilinos applications. I
don't if I am not setting it up correctly or what, but even when following
the one example that solves this particular problem it's still giving me a
lot of issues.
Is the example code attempting to factor in something else that I am not
thinking of? I am currently setting both the K and M matrices as the same
matrix-market file so that the amount of code in the example is changed as
little as possible. Is using the example a bad way to go about this?
Epetra_CrsMatrix *A, *B;
EpetraExt::MatrixMarketFileToCrsMatrix("../goodwin.mtx", Comm, A, false,
true);
EpetraExt::MatrixMarketFileToCrsMatrix("../goodwin.mtx", Comm, B, false,
true);
Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp(
const_cast<Epetra_CrsMatrix *>(A), false );
Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp(
const_cast<Epetra_CrsMatrix *>(B), false );
and I get the following output:
Reading MatrixMarket file ../goodwin.mtx
Sorting local nonzeros
Constructing the matrix
Inserting global values
Completing matrix fill
File Read time (secs): 0.131942
K dims:7320, 7320
M dims:7320, 7320
Anasazi::BlockKrylovSchurSolMgr::solve() caught unexpected exception from
Anasazi::BlockKrylovSchur::iterate() at iteration 0
/home/pamo8800/project/TRILINOS_INTEL_MPI/lib/cmake/Trilinos/../../../include/AnasaziSVQBOrthoManager.hpp:701:
Throw number = 1
Throw test that evaluated to true: info != 0
Anasazi::SVQBOrthoManager::findBasis(): Error code 2 from HEEV
Anasazi::BlockKrylovSchurSolMgr::solve() returning Unconverged with no
solutions
-Pate
On 4/2/2015 6:44:47 AM, Heroux, Michael A <maherou at sandia.gov> wrote:
Pate,
See if the following lines work to instantiate your K and M:
Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp(
const_cast<Epetra_CrsMatrix *>(A), false );
Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp(
const_cast<Epetra_CrsMatrix *>(M), false );
If that doesn’t work, check the integrity of A, K and M. See if the
attributes are sane by checking the dimensions and norms.
Mike
From: Pate Motter <pate.motter at colorado.edu>
Date: Monday, March 30, 2015 at 9:29 PM
To: Alicia M Klinvex <aklinvex at purdue.edu>
Cc: "trilinos-users at software.sandia.gov" <trilinos-users at software.sandia.gov
>
Subject: [EXTERNAL] Re: [Trilinos-Users] Smallest Eigenvalues
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/20150402/b40c2a40/attachment.html>
More information about the Trilinos-Users
mailing list