Function russell_lab::matrix::mat_inverse

pub fn mat_inverse(ai: &mut Matrix, a: &Matrix) -> Result<f64, StrError>

(dgetrf, dgetri) Computes the inverse of a square matrix and returns its determinant

ai := a⁻¹

See: https://www.netlib.org/lapack/explore-html/d3/d6a/dgetrf_8f.html

And: https://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f.html

Output

• ai – (m,m) inverse matrix
• Returns the matrix determinant

Input

• a – (m,m) matrix, symmetric or not

Examples

First – 2 x 2 square matrix

use russell_lab::{mat_inverse, Matrix, StrError};

fn main() -> Result<(), StrError> {
    // set matrix
    let mut a = Matrix::from(&[
        [-1.0,  1.5],
        [ 1.0, -1.0],
    ]);
    let a_copy = a.clone();

    // compute inverse matrix
    let mut ai = Matrix::new(2, 2);
    mat_inverse(&mut ai, &mut a)?;

    // compare with solution
    let ai_correct = "┌     ┐\n\
                      │ 2 3 │\n\
                      │ 2 2 │\n\
                      └     ┘";
    assert_eq!(format!("{}", ai), ai_correct);

    // check if a⋅ai == identity
    let (m, n) = a.dims();
    let mut a_ai = Matrix::new(m, m);
    for i in 0..m {
        for j in 0..m {
            for k in 0..n {
                a_ai.add(i, j, a_copy.get(i, k) * ai.get(k, j));
            }
        }
    }
    let identity = "┌     ┐\n\
                    │ 1 0 │\n\
                    │ 0 1 │\n\
                    └     ┘";
    assert_eq!(format!("{}", a_ai), identity);
    Ok(())
}

Second – 3 x 3 square matrix

use russell_lab::{mat_inverse, Matrix, StrError};

fn main() -> Result<(), StrError> {
    // set matrix
    let mut a = Matrix::from(&[
        [1.0, 2.0, 3.0],
        [0.0, 1.0, 4.0],
        [5.0, 6.0, 0.0],
    ]);
    let a_copy = a.clone();

    // compute inverse matrix
    let mut ai = Matrix::new(3, 3);
    mat_inverse(&mut ai, &mut a)?;

    // compare with solution
    let ai_correct = "┌             ┐\n\
                      │ -24  18   5 │\n\
                      │  20 -15  -4 │\n\
                      │  -5   4   1 │\n\
                      └             ┘";
    assert_eq!(format!("{}", ai), ai_correct);

    // check if a⋅ai == identity
    let (m, n) = a.dims();
    let mut a_ai = Matrix::new(m, m);
    for i in 0..m {
        for j in 0..m {
            for k in 0..n {
                a_ai.add(i, j, a_copy.get(i, k) * ai.get(k, j));
            }
        }
    }
    let identity = "┌       ┐\n\
                    │ 1 0 0 │\n\
                    │ 0 1 0 │\n\
                    │ 0 0 1 │\n\
                    └       ┘";
    assert_eq!(format!("{}", a_ai), identity);
    Ok(())
}