linfa_linalg/lobpcg/
mod.rs

1//!
2//! Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) is a matrix-free method for
3//! finding the large (or smallest) eigenvalues and the corresponding eigenvectors of a symmetric
4//! eigenvalue problem
5//! ```text
6//! A x = lambda x
7//! ```
8//! where A is symmetric and (x, lambda) the solution. It has the following advantages:
9//! * matrix free: does not require storing the coefficient matrix explicitely and only evaluates
10//!   matrix-vector products.
11//! * factorization-free: does not require any matrix decomposition
12//! * linear-convergence: theoretically guaranteed and practically observed
13//!
14//! See also the wikipedia article at [LOBPCG](https://en.wikipedia.org/wiki/LOBPCG)
15//!
16mod algorithm;
17mod eig;
18mod svd;
19
20use ndarray::prelude::*;
21use rand::distributions::Standard;
22use rand::prelude::*;
23
24pub use crate::{LinalgError, Order};
25pub use algorithm::lobpcg;
26pub use eig::{TruncatedEig, TruncatedEigIterator};
27pub use svd::{MagnitudeCorrection, TruncatedSvd};
28
29/// Generate random array
30pub(crate) fn random<A, Sh, D, R: Rng>(sh: Sh, mut rng: R) -> Array<A, D>
31where
32    A: NdFloat,
33    D: Dimension,
34    Sh: ShapeBuilder<Dim = D>,
35    Standard: Distribution<A>,
36{
37    ArrayBase::from_shape_fn(sh, |_| rng.gen::<A>())
38}
39
40/// The result of the eigensolver
41///
42/// In the best case the eigensolver has converged with a result better than the given threshold,
43/// then a `Ok` gives the eigenvalues, eigenvectors and norms. If an error ocurred
44/// during the process, it is returned in `Err` (together with the best result),
45/// as it could be of value. If there is no result at all, then the second field is `None`.
46/// This happens if the algorithm fails in an early stage, for example if the matrix `A` is not SPD
47pub type LobpcgResult<A> = std::result::Result<Lobpcg<A>, (LinalgError, Option<Lobpcg<A>>)>;
48pub struct Lobpcg<A> {
49    pub eigvals: Array1<A>,
50    pub eigvecs: Array2<A>,
51    pub rnorm: Vec<A>,
52}