DavidsonSolver API

[JensWehner]

This is still work in progress.

So I thought a bit about the davidson solver api and I have the following ideas:

a) Eigenvalue selection/target.

We write a selection class, which can be either constructed from an enum SortRule or a std::complex containing the targets . How to access the targets internally I am not sure yet, but we could let the class handle the selection internally. So instead of if(input==SortRule::max) do this ... elsif(input==SortRule::min) do this. Just have a newritzvalues=selectorclass.select(lastiteration)

Because selector class either needs a std::complex or an enum + number of eigenvalues JDsolver would need two interfaces for both cases, or we let the user initialize a selector class themselves.

I am actually not super happy with the current spectra API, because number of eigenvalues and how to select them should go together for me. @yixuan any thoughts?

Pros:

• For the user and the simple cases nothing changes
• Can be extended for other options easily

Cons:

• another layer of abstraction
• do not know

b) The MatrixProductOperator is exactly the same as in the rest of spectra. So no changes there

c) I think we discussed already that Symmetric and nonsymmetric case should have different function names,

d) Generalized case seems not to be too different form the specific case. We could just add the Bop template parameter as in the rest of spectra and default it to a NoOp, using a default imput parameter or enable_if template stuff. Or we can just write two differently named classes.

e) Preconditioner/Correction equation. This should be another template parameter. So we can combine the JD and Davidson method in one call. The CorrectionEquationSolver should have a solve() and info() method which return the next guess and inform if the calculation was successful. info() could return the Eigen::ComputationalInfo enum.

f) Initialguess: I think we should have a method solve(Selection ) and solvewithGuess(Selection, const Eigen::MatrixXd&initialspace). In the solve(Selection ) method, we should use the random guess.

g) Setting searchspace maximum, deflation, max number iterations should be handled via individual setter functions, with sensible defaults

h) Errorhandling, keep the Spectra convention. i.e.compute returns the number of converged eigenvalues. Furthermore a info() method returns an Eigen::ComputationalInfo enum.
eigenvalues() and eigenvectors() method should return the VectorXd and MatrixXd of eigenvalues and eigenvectors as const references.

i) Logging

To not add another template parameter, I would say we implement a virtual baseclass LoggerBase. LoggerBase has a function log(parameters) which is at the end of each iteration gets an update of the iterationcout, residual, searchspace_size maybe sth else. Users can derive their own logger from LoggerBase and use a setLogger(Logger) method to inject it into the solver.

Pros:

• Easy to use
• You can save your own data in the Logger class, e.g. time between iterations

Cons:

• The data from inside the solver class you can view via the loggeris fixed and cannot be extended from outside the library
• I am not sure how to handle float double std::complex datatypes here, because the logger has to store the same types as the solver.

[JensWehner]

To put this all together for the simple case it should look like this for the simple case:

 // Construct matrix operation object using the wrapper class DenseSymMatProd
    DenseSymMatProd<double> op(M);

    // Construct eigen solver object, requesting the largest three eigenvalues
    JacobiDavidson_SymSolver<double, DenseSymMatProd<double>> eigs(op);

eigs.init();
    int nconv = eigs.compute(SortRule::LargestAlge,N);

and for the difficult case:

 // Construct matrix operation object using the wrapper class DenseSymMatProd
    DenseSymMatProd<double> op(M);

    // Construct eigen solver object, requesting the largest three eigenvalues
    JacobiDavidson_SymSolver<double, DenseSymMatProd<double>> eigs(op);

eigs.init();
    int nconv = eigs.compute(SortRule::LargestAlge,N);

[JensWehner]

@yixuan what is the reason for the <double, DenseSymMatProd<double>>, what does the template parameter do, is it always the same as the one in DenseSymMatProd

[yixuan]

I'll read your comments in more details later, but let me first give a quick response to your last question. Why two template parameters? The reason is that users may provide their own matrix operation class, say MyMatProd, and they implement it as a regular class rather than a template. The MyMatProd class uses double as the element type, but SymEigsSolver doesn't know about this. Therefore, SymEigsSolver has to know both the element type and the operator type.

[yixuan]

For question a), can you elaborate a little bit about the meaning of vector<>? Isn't target a single number?

[JensWehner]

I have not made my notes really nice yet. I thought target can be an array of number, i.e. I could target different ones at the same time? I thought I saw sth. like this in the julia implementation. Ohh I was wrong, that makes it a lot easier then. :)

[JensWehner]

Concerning the first template parameter, you can require every Matrix Product class to define a T::value_type instead of providing it as a template parameter. It is a minor issue, but it makes your API a bit cleaner for the more common case of using a simple matrix.

[yixuan]

Yes I agree with that. It is a matter of where you want to put the type definition. I can add that typedef to existing classes so it won't affect existing code but JDSolver can utilize it.

[JensWehner]

I am still very unsure of how to make a good logger