algebra_linear 1.2.0

//! The reference LAPACK, for testing the kernels against.
//!
//! **Development only.** Enabled by the `dev_lapack` feature, which compiles
//! `c_code/reference_lapack.c` against a locally built LAPACK 3.12.1 and links
//! it; a normal build never sees this module. `ALGEBRA_LINEAR_DEV_LAPACK_ROOT`
//! points at the CMake build directory (see `build.rs`); there is deliberately no
//! default, since one would put a developer's layout into released source.
//!
//! The entry points are the same wrappers the consuming crate ships as its C shim,
//! transcribed from it — `c_code/reference_lapack.c` records the provenance and the
//! line range, and must be re-checked when that shim changes. They take `C_BOOL`
//! instead of Fortran character arguments and supply the hidden string lengths
//! themselves, which is why this module does not reproduce LAPACK's string ABI by
//! hand.
//!
//! # Column-major
//!
//! Everything here is column-major, matching the crate's storage. **LAPACKE would
//! not work**: its `LAPACKE_*` entry points are row-major and would transpose every
//! matrix, so a comparison against them would be meaningless.
//!
//! # Who calls this
//!
//! Only `crate::parity`, and only under `#[cfg(all(test, feature = "dev_lapack"))]`.
//! Nothing here is public API: the module is `pub(crate)` and the harness is a
//! private test module, so no part of it reaches `cargo doc` or docs.rs.
//!
//! # Safety
//!
//! Every function here hands raw pointers to LAPACK. The invariant each one relies
//! on is the reference's own: `n` (or `m` and `n`) describes the matrix the slices
//! hold, `ld` is that matrix's leading dimension, and the workspace is at least what
//! the `lwork = -1` query returned. The wrappers assert the slice lengths they can
//! check, and their only callers are the parity tests in `crate::parity`, which size
//! every buffer from the same `n` they pass.

#![allow(dead_code)]

use crate::Complex64;
use crate::constants::CcBool;

unsafe extern "C" {
    /// `dgesv`: solve `A X = B`. Takes `nrhs` right-hand sides
    pub(crate) fn c_dgesv(
        n: *const i32,
        nrhs: *const i32,
        a: *mut f64,
        lda: *const i32,
        ipiv: *mut i32,
        b: *mut f64,
        ldb: *const i32,
        info: *mut i32,
    );

    /// `dgetrf`: LU factorization with partial pivoting, `m x n`.
    pub(crate) fn c_dgetrf(m: *const i32, n: *const i32, a: *mut f64, lda: *const i32, ipiv: *mut i32, info: *mut i32);

    /// `dgetri`: inverse from a `dgetrf` factorization.
    pub(crate) fn c_dgetri(n: *const i32, a: *mut f64, lda: *const i32, ipiv: *const i32, work: *mut f64, lwork: *const i32, info: *mut i32);

    /// `dpotrf`: Cholesky factorization. `upper` selects the triangle to write.
    pub(crate) fn c_dpotrf(upper: CcBool, n: *const i32, a: *mut f64, lda: *const i32, info: *mut i32);

    /// `zgetrf`: LU factorization of a complex matrix.
    ///
    /// `COMPLEX64` in the C shim is `lapack_complex_double`, which is two
    /// `double`s in sequence, so `Complex64` (which is `#[repr(C)] { re, im }`)
    /// is layout-compatible and can be passed directly.
    pub(crate) fn c_zgetrf(m: *const i32, n: *const i32, a: *mut Complex64, lda: *const i32, ipiv: *mut i32, info: *mut i32);

    /// `zgetri`: complex inverse from a `zgetrf` factorization.
    pub(crate) fn c_zgetri(
        n: *const i32,
        a: *mut Complex64,
        lda: *const i32,
        ipiv: *const i32,
        work: *mut Complex64,
        lwork: *const i32,
        info: *mut i32,
    );

    /// `zgesv`: complex solve `A X = B`.
    pub(crate) fn c_zgesv(
        n: *const i32,
        nrhs: *const i32,
        a: *mut Complex64,
        lda: *const i32,
        ipiv: *mut i32,
        b: *mut Complex64,
        ldb: *const i32,
        info: *mut i32,
    );

    /// `zpotrf`: Hermitian positive-definite Cholesky.
    pub(crate) fn c_zpotrf(upper: CcBool, n: *const i32, a: *mut Complex64, lda: *const i32, info: *mut i32);

    /// `zheev`: eigenvalues (and optionally eigenvectors) of a Hermitian matrix.
    pub(crate) fn c_zheev(
        calc_v: CcBool,
        upper: CcBool,
        n: *const i32,
        a: *mut Complex64,
        lda: *const i32,
        w: *mut f64,
        work: *mut Complex64,
        lwork: *const i32,
        rwork: *mut f64,
        info: *mut i32,
    );

    /// `dsyev`: eigenvalues (and optionally eigenvectors) of a symmetric matrix.
    pub(crate) fn c_dsyev(
        calc_v: CcBool,
        upper: CcBool,
        n: *const i32,
        a: *mut f64,
        lda: *const i32,
        w: *mut f64,
        work: *mut f64,
        lwork: *const i32,
        info: *mut i32,
    );

    /// `dgeev`: eigenvalues and eigenvectors of a general matrix.
    ///
    /// Eigenvalues come back as a real/imaginary pair because a real matrix can
    /// have complex conjugate pairs. Real eigenvectors are not returned for a
    /// complex pair, and `wi` must be inspected to know which case applies.
    pub(crate) fn c_dgeev(
        calc_vl: CcBool,
        calc_vr: CcBool,
        n: *const i32,
        a: *mut f64,
        lda: *const i32,
        wr: *mut f64,
        wi: *mut f64,
        vl: *mut f64,
        ldvl: *const i32,
        vr: *mut f64,
        ldvr: *const i32,
        work: *mut f64,
        lwork: *const i32,
        info: *mut i32,
    );

    /// `dgesvd`: singular value decomposition, `m x n`.
    pub(crate) fn c_dgesvd(
        jobu_code: i32,
        jobvt_code: i32,
        m: *const i32,
        n: *const i32,
        a: *mut f64,
        lda: *const i32,
        s: *mut f64,
        u: *mut f64,
        ldu: *const i32,
        vt: *mut f64,
        ldvt: *const i32,
        work: *mut f64,
        lwork: *const i32,
        info: *mut i32,
    );
}

/// The `jobu`/`jobvt` codes the shim's `c_dgesvd` expects, mirroring the Fortran
/// characters `A`, `S`, `O` and `N`.
pub(crate) mod svd_code {
    pub const ALL: i32 = 0; // 'A'
    pub const SKINNY: i32 = 1; // 'S'
    pub const OVERWRITE: i32 = 2; // 'O'
    pub const NONE: i32 = 3; // 'N'
}

/// Convenience wrappers that handle the workspace query, so the tests do not each
/// repeat it.
///
/// LAPACK's convention is to call with `lwork = -1` to be told the workspace size
/// in `work[0]`. Getting that wrong is a silent buffer overrun, so it is done once
/// here.
pub(crate) mod query {
    use super::*;
    // `use super::*` brings in this module's own items, not its imports, so the
    // constants have to be named again here.
    use crate::constants::{C_FALSE, C_TRUE};

    /// Runs `dgeev` with a correctly sized workspace.
    ///
    /// `a` is `n x n` column-major and is destroyed. `wr`/`wi` receive the
    /// eigenvalues, `vr` the right eigenvectors (if `want_v`).
    #[allow(clippy::too_many_arguments)]
    pub(crate) fn dgeev(a: &mut [f64], n: usize, wr: &mut [f64], wi: &mut [f64], vr: Option<&mut [f64]>) -> i32 {
        assert_eq!(a.len(), n * n, "reference dgeev: a must be n*n");
        assert_eq!(wr.len(), n, "reference dgeev: wr must have n entries");
        assert_eq!(wi.len(), n, "reference dgeev: wi must have n entries");
        if let Some(v) = &vr {
            assert_eq!(v.len(), n * n, "reference dgeev: vr must be n*n when requested");
        }
        let n_i32 = n as i32;
        let ld = n_i32;
        let mut info = 0i32;
        let (want_vr, ldvr) = match &vr {
            Some(_) => (C_TRUE, n_i32),
            None => (C_FALSE, 1),
        };
        let vr_ptr = match vr {
            Some(v) => v.as_mut_ptr(),
            None => std::ptr::null_mut(),
        };

        // Workspace query.
        let mut work = vec![0.0f64; 1];
        let mut lwork = -1i32;
        // SAFETY: every pointer comes from a slice whose length the assertions at the top
        // of this function pin to the LAPACK argument it feeds, `ld` is that matrix's
        // leading dimension, and `work` is at least the size the `lwork = -1` query returned.
        unsafe {
            c_dgeev(
                C_FALSE,
                want_vr,
                &n_i32,
                a.as_mut_ptr(),
                &ld,
                wr.as_mut_ptr(),
                wi.as_mut_ptr(),
                std::ptr::null_mut(),
                &1i32,
                vr_ptr,
                &ldvr,
                work.as_mut_ptr(),
                &lwork,
                &mut info,
            );
        }
        if info != 0 {
            return info;
        }
        let need = work[0].max(1.0) as usize;
        let mut work = vec![0.0f64; need];
        lwork = need as i32;
        // SAFETY: every pointer comes from a slice whose length the assertions at the top
        // of this function pin to the LAPACK argument it feeds, `ld` is that matrix's
        // leading dimension, and `work` is at least the size the `lwork = -1` query returned.
        unsafe {
            c_dgeev(
                C_FALSE,
                want_vr,
                &n_i32,
                a.as_mut_ptr(),
                &ld,
                wr.as_mut_ptr(),
                wi.as_mut_ptr(),
                std::ptr::null_mut(),
                &1i32,
                vr_ptr,
                &ldvr,
                work.as_mut_ptr(),
                &lwork,
                &mut info,
            );
        }
        info
    }

    /// Runs `dsyev` in the lower triangle, with a correctly sized workspace.
    pub(crate) fn dsyev(a: &mut [f64], n: usize, w: &mut [f64], want_v: bool) -> i32 {
        assert_eq!(a.len(), n * n, "reference dsyev: a must be n*n");
        assert_eq!(w.len(), n, "reference dsyev: w must have n entries");
        let n_i32 = n as i32;
        let ld = n_i32;
        let mut info = 0i32;
        let jobz = if want_v { C_TRUE } else { C_FALSE };

        let mut work = vec![0.0f64; 1];
        let mut lwork = -1i32;
        // SAFETY: every pointer comes from a slice whose length the assertions at the top
        // of this function pin to the LAPACK argument it feeds, `ld` is that matrix's
        // leading dimension, and `work` is at least the size the `lwork = -1` query returned.
        unsafe {
            c_dsyev(
                jobz,
                C_FALSE,
                &n_i32,
                a.as_mut_ptr(),
                &ld,
                w.as_mut_ptr(),
                work.as_mut_ptr(),
                &lwork,
                &mut info,
            );
        }
        if info != 0 {
            return info;
        }
        let need = work[0].max(1.0) as usize;
        let mut work = vec![0.0f64; need];
        lwork = need as i32;
        // SAFETY: every pointer comes from a slice whose length the assertions at the top
        // of this function pin to the LAPACK argument it feeds, `ld` is that matrix's
        // leading dimension, and `work` is at least the size the `lwork = -1` query returned.
        unsafe {
            c_dsyev(
                jobz,
                C_FALSE,
                &n_i32,
                a.as_mut_ptr(),
                &ld,
                w.as_mut_ptr(),
                work.as_mut_ptr(),
                &lwork,
                &mut info,
            );
        }
        info
    }

    /// Runs `zheev` in the lower triangle, with a correctly sized workspace.
    pub(crate) fn zheev(a: &mut [Complex64], n: usize, w: &mut [f64], want_v: bool) -> i32 {
        assert_eq!(a.len(), n * n, "reference zheev: a must be n*n");
        assert_eq!(w.len(), n, "reference zheev: w must have n entries");
        let n_i32 = n as i32;
        let ld = n_i32;
        let mut info = 0i32;
        let jobz = if want_v { C_TRUE } else { C_FALSE };

        let mut work = vec![Complex64::new(0.0, 0.0); 1];
        let mut rwork = vec![0.0f64; (3 * n).max(1)];
        let mut lwork = -1i32;
        // SAFETY: every pointer comes from a slice whose length the assertions at the top
        // of this function pin to the LAPACK argument it feeds, `ld` is that matrix's
        // leading dimension, and `work` is at least the size the `lwork = -1` query returned.
        unsafe {
            c_zheev(
                jobz,
                C_FALSE,
                &n_i32,
                a.as_mut_ptr(),
                &ld,
                w.as_mut_ptr(),
                work.as_mut_ptr(),
                &lwork,
                rwork.as_mut_ptr(),
                &mut info,
            );
        }
        if info != 0 {
            return info;
        }
        let need = work[0].re.max(1.0) as usize;
        let mut work = vec![Complex64::new(0.0, 0.0); need];
        lwork = need as i32;
        // SAFETY: every pointer comes from a slice whose length the assertions at the top
        // of this function pin to the LAPACK argument it feeds, `ld` is that matrix's
        // leading dimension, and `work` is at least the size the `lwork = -1` query returned.
        unsafe {
            c_zheev(
                jobz,
                C_FALSE,
                &n_i32,
                a.as_mut_ptr(),
                &ld,
                w.as_mut_ptr(),
                work.as_mut_ptr(),
                &lwork,
                rwork.as_mut_ptr(),
                &mut info,
            );
        }
        info
    }

    /// Runs `dgeev`-style workspace handling for `dgesvd`, returning `s`.
    #[allow(clippy::too_many_arguments)]
    pub(crate) fn dgesvd(a: &mut [f64], m: usize, n: usize, s: &mut [f64], u: &mut [f64], vt: &mut [f64]) -> i32 {
        let (m_i32, n_i32) = (m as i32, n as i32);
        let min_mn = m.min(n);
        assert_eq!(a.len(), m * n, "reference dgesvd: a must be m*n");
        assert_eq!(u.len(), m * m, "reference dgesvd: u must be m*m");
        assert_eq!(vt.len(), n * n, "reference dgesvd: vt must be n*n");
        assert_eq!(s.len(), m.min(n), "reference dgesvd: s must be min(m,n)");
        let mut info = 0i32;

        let mut work = vec![0.0f64; 1];
        let mut lwork = -1i32;
        // SAFETY: every pointer comes from a slice whose length the assertions at the top
        // of this function pin to the LAPACK argument it feeds, `ld` is that matrix's
        // leading dimension, and `work` is at least the size the `lwork = -1` query returned.
        unsafe {
            c_dgesvd(
                svd_code::ALL,
                svd_code::ALL,
                &m_i32,
                &n_i32,
                a.as_mut_ptr(),
                &m_i32,
                s.as_mut_ptr(),
                u.as_mut_ptr(),
                &m_i32,
                vt.as_mut_ptr(),
                &n_i32,
                work.as_mut_ptr(),
                &lwork,
                &mut info,
            );
        }
        if info != 0 {
            return info;
        }
        let need = work[0].max(1.0) as usize;
        let mut work = vec![0.0f64; need];
        lwork = need as i32;
        let _ = min_mn;
        // SAFETY: every pointer comes from a slice whose length the assertions at the top
        // of this function pin to the LAPACK argument it feeds, `ld` is that matrix's
        // leading dimension, and `work` is at least the size the `lwork = -1` query returned.
        unsafe {
            c_dgesvd(
                svd_code::ALL,
                svd_code::ALL,
                &m_i32,
                &n_i32,
                a.as_mut_ptr(),
                &m_i32,
                s.as_mut_ptr(),
                u.as_mut_ptr(),
                &m_i32,
                vt.as_mut_ptr(),
                &n_i32,
                work.as_mut_ptr(),
                &lwork,
                &mut info,
            );
        }
        info
    }
}