Function russell_lab::matrix::mat_eigen

pub fn mat_eigen(
    l_real: &mut Vector,
    l_imag: &mut Vector,
    v_real: &mut Matrix,
    v_imag: &mut Matrix,
    a: &mut Matrix
) -> Result<(), StrError>

(dgeev) Performs the eigen-decomposition of a square matrix

Computes the eigenvalues l and right eigenvectors v, such that:

a ⋅ vj = lj ⋅ vj

where lj is the component j of l and vj is the column j of v.

See also: https://www.netlib.org/lapack/explore-html/d9/d28/dgeev_8f.html

Output

• l_real – (m) eigenvalues; real part
• l_imag – (m) eigenvalues; imaginary part
• v_real – (m,m) right eigenvectors (as columns); real part
• v_imag – (m,m) right eigenvectors (as columns); imaginary part

Input

• a – (m,m) general matrix (will be modified)

Note

• The matrix a will be modified

Similarity transformation

The eigen-decomposition leads to a similarity transformation like so:

a = v⋅λ⋅v⁻¹

where v is a matrix whose columns are the m linearly independent eigenvectors of a, and λ is a matrix whose diagonal are the eigenvalues of a. Thus, the following is valid:

a⋅v = v⋅λ

Let us define the error err as follows:

err := a⋅v - v⋅λ

Examples

use russell_lab::*;

fn main() -> Result<(), StrError> {
    // set matrix
    let data = [[2.0, 0.0, 0.0], [0.0, 3.0, 4.0], [0.0, 4.0, 9.0]];
    let mut a = Matrix::from(&data);

    // allocate output arrays
    let m = a.nrow();
    let mut l_real = Vector::new(m);
    let mut l_imag = Vector::new(m);
    let mut v_real = Matrix::new(m, m);
    let mut v_imag = Matrix::new(m, m);

    // perform the eigen-decomposition
    mat_eigen(&mut l_real, &mut l_imag, &mut v_real, &mut v_imag, &mut a)?;

    // check results
    assert_eq!(
        format!("{:.1}", l_real),
        "┌      ┐\n\
         │ 11.0 │\n\
         │  1.0 │\n\
         │  2.0 │\n\
         └      ┘"
    );
    assert_eq!(
        format!("{}", l_imag),
        "┌   ┐\n\
         │ 0 │\n\
         │ 0 │\n\
         │ 0 │\n\
         └   ┘"
    );

    // check eigen-decomposition (similarity transformation) of a
    // symmetric matrix with real-only eigenvalues and eigenvectors
    let a_copy = Matrix::from(&data);
    let lam = Matrix::diagonal(l_real.as_data());
    let mut a_v = Matrix::new(m, m);
    let mut v_l = Matrix::new(m, m);
    let mut err = Matrix::filled(m, m, f64::MAX);
    mat_mat_mul(&mut a_v, 1.0, &a_copy, &v_real, 0.0)?;
    mat_mat_mul(&mut v_l, 1.0, &v_real, &lam, 0.0)?;
    mat_add(&mut err, 1.0, &a_v, -1.0, &v_l)?;
    approx_eq(mat_norm(&err, Norm::Max), 0.0, 1e-15);
    Ok(())
}