-
Notifications
You must be signed in to change notification settings - Fork 0
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Implement DPR #21
Implement DPR #21
Conversation
Codecov Report
@@ Coverage Diff @@
## master #21 +/- ##
========================================
Coverage ? 92.0%
========================================
Files ? 42
Lines ? 1605
Branches ? 0
========================================
Hits ? 1478
Misses ? 127
Partials ? 0 Continue to review full report at Codecov.
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Great work!
I think you should add a testcase (just copy the JDSymEigsBase one) where the JDSymEigsDPR class is instantiated
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Nice work ! There are a few things I don't quite get in the correction that might need checking
include/Spectra/JDSymEigsDPR.h
Outdated
/// \return new correction vectors. | ||
Matrix CalculateCorrectionVector() const final | ||
{ | ||
Index nresidues = search_space_.size(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
the number of correction vectors is not equal to the search space size. We have a member variable called initial_search_space_size_
that also defines the number of correction vector we want to append. We could make an additional member called size_update_
or similar in the base class to be explicit though
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Lets use this size_update_
class variable
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@JensWehner can we create this size_update_
variable as a protected member variable in the base class #20
include/Spectra/JDSymEigsDPR.h
Outdated
{ | ||
Index nresidues = search_space_.size(); | ||
const Matrix& residues = Residues(); | ||
const Vector& eigenvalues = ritz_pairs_.RitzValues(); |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
there is an eigenvalues() method in the base class called like that. So at least we should rename that variable to eigvals
or similar. But we could also change the eigenvalue()
method to :
Vector eigenvalues(nev=number_eigenvalues_) const { return ritz_pairs_.RitzValues().head(nev); }
so that we could call :
const Vector& eigvals = eigenvalues(size_update_);
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks nice work, just a few things.
include/Spectra/JDSymEigsDPR.h
Outdated
for (Index k = 0; k < this->size_update_; k++) | ||
{ | ||
Vector tmp = eigvals(k) - diagonal_.array(); | ||
correction.col(k) = residues.col(k).array() / tmp.array(); | ||
} |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think you can do this without a for loop using the colwise() stuff from eigen.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I don't see how to perform that operation colwise
without passing some kind of lambda to function to operate with the residues and the diagonal
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
hmm... actually the for loop has to check if that RitzValue is already converged. So we do not append corrections for already converged vectors. @NicoRenaud correct?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Ideally we should happen only the correction vector of the unconverged eigenvalues indeed. But it would be great to have something that runs first. Then we can optimize.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
in order to avoid the loop we should create a matrix m_ij = 1. / (diag_i - lambda_j)
and do a element wise multiplication. I don't know how to create this matrix without a loop though but it might be possible
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
If you append correction vectors for already converged eigenvalues, you are endangering your orthogonalisation.
I was wrong about the for loop, it is best the way it is.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ConvergedEigenvalues()
is a function of RitzPairs, which tells you which ritzvalues are already converged.
include/Spectra/JDSymEigsDPR.h
Outdated
/// | ||
|
||
template <typename OpType> | ||
class JDSymEigsDPR : public JDSymEigsBase<OpType> |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
actually the name should be DavidsonEigsSym, I think DPR is the standard correction. @NicoRenaud ?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
DavidsonSymEig might be a better name indeed
@JensWehner I have renamed the DPR file and the example and merge the |
DPR implementation for Davidson