oxicuda_solver/dense/

ldlt.rs

//! Symmetric Indefinite Factorization (LDL^T / Bunch-Kaufman).
//!
//! Computes `P * A * P^T = L * D * L^T` where:
//! - P is a permutation matrix (encoded in pivot_info)
//! - L is unit lower triangular
//! - D is block diagonal with 1x1 and 2x2 blocks
//!
//! # Algorithm — Bunch-Kaufman Pivoting
//!
//! At each step, the algorithm decides whether to use a 1x1 or 2x2 pivot
//! by examining the magnitudes of diagonal and off-diagonal elements:
//!
//! 1. Let `alpha = (1 + sqrt(17)) / 8` (~0.6404).
//! 2. Find `lambda = max |A[i, k]|` for i != k (largest off-diagonal in column k).
//! 3. If `|A[k,k]| >= alpha * lambda`: use 1x1 pivot at position k.
//! 4. Otherwise, find `sigma = max |A[i, r]|` for i != r, where r is the row
//!    achieving lambda.
//! 5. If `|A[k,k]| * sigma >= alpha * lambda^2`: use 1x1 pivot.
//! 6. If `|A[r,r]| >= alpha * sigma`: use 1x1 pivot at position r (with swap).
//! 7. Otherwise: use 2x2 pivot at positions (k, r).
//!
//! This guarantees bounded element growth in L.

#![allow(dead_code)]

use oxicuda_blas::types::{FillMode, GpuFloat};
use oxicuda_memory::DeviceBuffer;

use crate::error::{SolverError, SolverResult};
use crate::handle::SolverHandle;

// ---------------------------------------------------------------------------
// GpuFloat <-> f64 conversion helpers
// ---------------------------------------------------------------------------

fn to_f64<T: GpuFloat>(val: T) -> f64 {
    if T::SIZE == 4 {
        f32::from_bits(val.to_bits_u64() as u32) as f64
    } else {
        f64::from_bits(val.to_bits_u64())
    }
}

fn from_f64<T: GpuFloat>(val: f64) -> T {
    if T::SIZE == 4 {
        T::from_bits_u64(u64::from((val as f32).to_bits()))
    } else {
        T::from_bits_u64(val.to_bits())
    }
}

// ---------------------------------------------------------------------------
// Constants
// ---------------------------------------------------------------------------

/// Bunch-Kaufman pivot threshold: alpha = (1 + sqrt(17)) / 8.
const BUNCH_KAUFMAN_ALPHA: f64 = 0.6403882032022076;

// ---------------------------------------------------------------------------
// Result type
// ---------------------------------------------------------------------------

/// LDL^T factorization result.
///
/// Contains pivot information describing the block diagonal structure of D:
/// - `pivot_info[k] > 0`: 1x1 pivot block at position k; rows/cols were swapped
///   with row/col `pivot_info[k] - 1`.
/// - `pivot_info[k] < 0`: start of a 2x2 pivot block at positions (k, k+1);
///   rows/cols were swapped with row/col `(-pivot_info[k]) - 1`.
pub struct LdltResult {
    /// Pivot block sizes: positive for 1x1, negative for start of 2x2.
    pub pivot_info: DeviceBuffer<i32>,
}

impl std::fmt::Debug for LdltResult {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        f.debug_struct("LdltResult")
            .field("pivot_info_len", &self.pivot_info.len())
            .finish()
    }
}

// ---------------------------------------------------------------------------
// Public API
// ---------------------------------------------------------------------------

/// Computes the Bunch-Kaufman LDL^T factorization of a symmetric indefinite matrix.
///
/// On exit, the specified triangle of `a` is overwritten with L and D:
/// - The unit lower triangular factor L (with unit diagonal implicit).
/// - The block diagonal factor D (1x1 and 2x2 blocks on the diagonal/super-diagonal).
///
/// # Arguments
///
/// * `handle` — solver handle.
/// * `a` — symmetric matrix (n x n, column-major), overwritten with L and D.
/// * `n` — matrix dimension.
/// * `uplo` — which triangle to read/write (Lower or Upper).
///
/// # Returns
///
/// An [`LdltResult`] containing the pivot information.
///
/// # Errors
///
/// Returns [`SolverError::DimensionMismatch`] for invalid dimensions.
/// Returns [`SolverError::SingularMatrix`] if the matrix is singular.
pub fn ldlt<T: GpuFloat>(
    handle: &mut SolverHandle,
    a: &mut DeviceBuffer<T>,
    n: usize,
    uplo: FillMode,
) -> SolverResult<LdltResult> {
    if n == 0 {
        let pivot_info = DeviceBuffer::<i32>::zeroed(0)?;
        return Ok(LdltResult { pivot_info });
    }
    if a.len() < n * n {
        return Err(SolverError::DimensionMismatch(format!(
            "ldlt: buffer too small ({} < {})",
            a.len(),
            n * n
        )));
    }
    if uplo == FillMode::Full {
        return Err(SolverError::DimensionMismatch(
            "ldlt: uplo must be Upper or Lower, not Full".into(),
        ));
    }

    // Workspace for the host-side factorization.
    let ws = n * n * std::mem::size_of::<f64>();
    handle.ensure_workspace(ws)?;

    // Read the matrix into host memory for the factorization.
    let mut a_host = vec![0.0_f64; n * n];
    read_device_to_host(a, &mut a_host, n * n)?;

    // Perform the Bunch-Kaufman factorization on the host.
    let mut ipiv = vec![0_i32; n];
    bunch_kaufman_factorize(&mut a_host, n, uplo, &mut ipiv)?;

    // Write back to device.
    let a_device: Vec<T> = a_host.iter().map(|&v| from_f64(v)).collect();
    write_host_to_device(a, &a_device, n * n)?;

    let mut pivot_info = DeviceBuffer::<i32>::zeroed(n)?;
    write_host_to_device_i32(&mut pivot_info, &ipiv, n)?;

    Ok(LdltResult { pivot_info })
}

/// Solves `A * x = b` using the LDL^T factorization.
///
/// The LDL^T factors must have been computed by [`ldlt`].
///
/// Algorithm:
/// 1. Apply row permutations P to b.
/// 2. Forward substitution: solve `L * y = P * b`.
/// 3. Block diagonal solve: solve `D * z = y`.
/// 4. Backward substitution: solve `L^T * w = z`.
/// 5. Apply column permutations P^T to w to get x.
///
/// # Arguments
///
/// * `handle` — solver handle.
/// * `a` — LDL^T-factored matrix (output of `ldlt`).
/// * `pivot_info` — pivot information from `ldlt`.
/// * `b` — right-hand side matrix (n x nrhs), overwritten with solution.
/// * `n` — system dimension.
/// * `nrhs` — number of right-hand side columns.
/// * `uplo` — which triangle contains the factor.
///
/// # Errors
///
/// Returns [`SolverError::DimensionMismatch`] for invalid dimensions.
pub fn ldlt_solve<T: GpuFloat>(
    handle: &mut SolverHandle,
    a: &DeviceBuffer<T>,
    pivot_info: &DeviceBuffer<i32>,
    b: &mut DeviceBuffer<T>,
    n: usize,
    nrhs: usize,
    uplo: FillMode,
) -> SolverResult<()> {
    if n == 0 || nrhs == 0 {
        return Ok(());
    }
    if a.len() < n * n {
        return Err(SolverError::DimensionMismatch(
            "ldlt_solve: factor buffer too small".into(),
        ));
    }
    if pivot_info.len() < n {
        return Err(SolverError::DimensionMismatch(
            "ldlt_solve: pivot_info buffer too small".into(),
        ));
    }
    if b.len() < n * nrhs {
        return Err(SolverError::DimensionMismatch(
            "ldlt_solve: B buffer too small".into(),
        ));
    }

    // Workspace.
    let ws = (n * n + n * nrhs) * std::mem::size_of::<f64>();
    handle.ensure_workspace(ws)?;

    // Read to host.
    let mut a_host = vec![0.0_f64; n * n];
    read_device_to_host(a, &mut a_host, n * n)?;

    let mut ipiv = vec![0_i32; n];
    read_device_to_host_i32(pivot_info, &mut ipiv, n)?;

    let mut b_host = vec![0.0_f64; n * nrhs];
    read_device_to_host(b, &mut b_host, n * nrhs)?;

    // Solve on host.
    bunch_kaufman_solve(&a_host, &ipiv, &mut b_host, n, nrhs, uplo)?;

    // Write back.
    let b_device: Vec<T> = b_host.iter().map(|&v| from_f64(v)).collect();
    write_host_to_device(b, &b_device, n * nrhs)?;

    Ok(())
}

// ---------------------------------------------------------------------------
// Bunch-Kaufman factorization (host-side)
// ---------------------------------------------------------------------------

/// Bunch-Kaufman factorization on a host-side column-major matrix.
fn bunch_kaufman_factorize(
    a: &mut [f64],
    n: usize,
    uplo: FillMode,
    ipiv: &mut [i32],
) -> SolverResult<()> {
    match uplo {
        FillMode::Lower => bunch_kaufman_lower(a, n, ipiv),
        FillMode::Upper => {
            // Normalize to a lower-storage representation and reuse the same
            // decomposition path to avoid divergent structural implementations.
            mirror_upper_to_lower(a, n);
            bunch_kaufman_lower(a, n, ipiv)
        }
        FillMode::Full => Err(SolverError::DimensionMismatch(
            "ldlt: uplo must be Lower or Upper".into(),
        )),
    }
}

/// Lower-triangular Bunch-Kaufman: P*A*P^T = L*D*L^T.
fn bunch_kaufman_lower(a: &mut [f64], n: usize, ipiv: &mut [i32]) -> SolverResult<()> {
    let mut k = 0;

    while k < n {
        // Find the largest off-diagonal in column k (below diagonal).
        let (lambda, r_idx) = column_max_offdiag(a, n, k, true);

        let akk = a[k * n + k].abs();

        if akk < 1e-300 && lambda < 1e-300 {
            // Entire column is zero — matrix is singular at this point.
            return Err(SolverError::SingularMatrix);
        }

        if akk >= BUNCH_KAUFMAN_ALPHA * lambda {
            // Case 1: Use 1x1 pivot at position k.
            perform_1x1_pivot_lower(a, n, k);
            ipiv[k] = (k + 1) as i32; // 1-based, positive => 1x1
            k += 1;
        } else {
            // Find sigma = max |A[i, r]| for i != r.
            let (sigma, _) = column_max_offdiag(a, n, r_idx, true);

            if akk * sigma >= BUNCH_KAUFMAN_ALPHA * lambda * lambda {
                // Case 2: 1x1 pivot at k is fine.
                perform_1x1_pivot_lower(a, n, k);
                ipiv[k] = (k + 1) as i32;
                k += 1;
            } else if a[r_idx * n + r_idx].abs() >= BUNCH_KAUFMAN_ALPHA * sigma {
                // Case 3: 1x1 pivot at r (swap k <-> r first).
                if r_idx != k {
                    swap_rows_and_cols(a, n, k, r_idx);
                }
                perform_1x1_pivot_lower(a, n, k);
                ipiv[k] = (r_idx + 1) as i32;
                k += 1;
            } else {
                // Case 4: 2x2 pivot at (k, k+1).
                if k + 1 >= n {
                    // Edge case: can't form 2x2 block at the last row.
                    perform_1x1_pivot_lower(a, n, k);
                    ipiv[k] = (k + 1) as i32;
                    k += 1;
                } else {
                    if r_idx != k + 1 {
                        swap_rows_and_cols(a, n, k + 1, r_idx);
                    }
                    perform_2x2_pivot_lower(a, n, k)?;
                    ipiv[k] = -((r_idx + 1) as i32); // Negative => start of 2x2
                    ipiv[k + 1] = ipiv[k];
                    k += 2;
                }
            }
        }
    }

    Ok(())
}

/// Upper-triangular Bunch-Kaufman: P*A*P^T = U^T*D*U.
fn bunch_kaufman_upper(a: &mut [f64], n: usize, ipiv: &mut [i32]) -> SolverResult<()> {
    if n == 0 {
        return Ok(());
    }

    let mut k = n;

    while k > 0 {
        let col = k - 1;
        let (lambda, r_idx) = column_max_offdiag(a, n, col, false);
        let akk = a[col * n + col].abs();

        if akk < 1e-300 && lambda < 1e-300 {
            return Err(SolverError::SingularMatrix);
        }

        if akk >= BUNCH_KAUFMAN_ALPHA * lambda {
            ipiv[col] = (col + 1) as i32;
            k -= 1;
        } else {
            let (sigma, _) = column_max_offdiag(a, n, r_idx, false);

            if akk * sigma >= BUNCH_KAUFMAN_ALPHA * lambda * lambda {
                ipiv[col] = (col + 1) as i32;
                k -= 1;
            } else if a[r_idx * n + r_idx].abs() >= BUNCH_KAUFMAN_ALPHA * sigma {
                if r_idx != col {
                    swap_rows_and_cols(a, n, col, r_idx);
                }
                ipiv[col] = (r_idx + 1) as i32;
                k -= 1;
            } else {
                if col == 0 {
                    ipiv[col] = (col + 1) as i32;
                    k -= 1;
                } else {
                    let col2 = col - 1;
                    if r_idx != col2 {
                        swap_rows_and_cols(a, n, col2, r_idx);
                    }
                    ipiv[col] = -((r_idx + 1) as i32);
                    ipiv[col2] = ipiv[col];
                    k -= 2;
                }
            }
        }
    }

    Ok(())
}

// ---------------------------------------------------------------------------
// Pivot operations
// ---------------------------------------------------------------------------

/// Finds the maximum absolute off-diagonal element in a column.
/// Returns (max_value, row_index).
fn column_max_offdiag(a: &[f64], n: usize, col: usize, lower: bool) -> (f64, usize) {
    let mut max_val = 0.0_f64;
    let mut max_idx = col;

    if lower {
        for i in (col + 1)..n {
            let val = a[col * n + i].abs();
            if val > max_val {
                max_val = val;
                max_idx = i;
            }
        }
    } else {
        for i in 0..col {
            let val = a[col * n + i].abs();
            if val > max_val {
                max_val = val;
                max_idx = i;
            }
        }
    }

    (max_val, max_idx)
}

/// Swaps rows and columns i and j in a symmetric matrix (column-major).
fn swap_rows_and_cols(a: &mut [f64], n: usize, i: usize, j: usize) {
    if i == j {
        return;
    }
    // Swap row i and row j.
    for col in 0..n {
        a.swap(col * n + i, col * n + j);
    }
    // Swap col i and col j.
    for row in 0..n {
        a.swap(i * n + row, j * n + row);
    }
}

/// Performs a 1x1 pivot step at position k for lower-triangular factorization.
fn perform_1x1_pivot_lower(a: &mut [f64], n: usize, k: usize) {
    let akk = a[k * n + k];
    if akk.abs() < 1e-300 {
        return; // Cannot pivot on zero.
    }
    let inv_akk = 1.0 / akk;

    // Scale column k below diagonal: L[i, k] = A[i, k] / A[k, k].
    for i in (k + 1)..n {
        a[k * n + i] *= inv_akk;
    }

    // Update trailing submatrix: A[i, j] -= L[i, k] * D[k, k] * L[j, k].
    for j in (k + 1)..n {
        let ljk = a[k * n + j];
        for i in j..n {
            let lik = a[k * n + i];
            a[j * n + i] -= lik * akk * ljk;
        }
    }
}

/// Performs a 2x2 pivot step at positions (k, k+1) for lower-triangular factorization.
fn perform_2x2_pivot_lower(a: &mut [f64], n: usize, k: usize) -> SolverResult<()> {
    if k + 1 >= n {
        return Err(SolverError::InternalError(
            "ldlt: 2x2 pivot at boundary".into(),
        ));
    }

    // Extract 2x2 block D.
    let d11 = a[k * n + k];
    let d21 = a[k * n + (k + 1)];
    let d22 = a[(k + 1) * n + (k + 1)];

    // Invert the 2x2 block: D^{-1} = adj(D) / det(D).
    let det = d11 * d22 - d21 * d21;
    if det.abs() < 1e-300 {
        return Err(SolverError::SingularMatrix);
    }
    let inv_det = 1.0 / det;

    // Compute L columns below the 2x2 block.
    // [L[i,k], L[i,k+1]] = [A[i,k], A[i,k+1]] * D^{-1}
    for i in (k + 2)..n {
        let aik = a[k * n + i];
        let aik1 = a[(k + 1) * n + i];

        a[k * n + i] = (d22 * aik - d21 * aik1) * inv_det;
        a[(k + 1) * n + i] = (-d21 * aik + d11 * aik1) * inv_det;
    }

    // Update trailing submatrix.
    for j in (k + 2)..n {
        let ljk = a[k * n + j];
        let ljk1 = a[(k + 1) * n + j];

        for i in j..n {
            let lik = a[k * n + i];
            let lik1 = a[(k + 1) * n + i];

            // A[i,j] -= L[i,k]*D[k,k]*L[j,k] + L[i,k]*D[k,k+1]*L[j,k+1]
            //         + L[i,k+1]*D[k+1,k]*L[j,k] + L[i,k+1]*D[k+1,k+1]*L[j,k+1]
            a[j * n + i] -=
                lik * d11 * ljk + lik * d21 * ljk1 + lik1 * d21 * ljk + lik1 * d22 * ljk1;
        }
    }

    Ok(())
}

// ---------------------------------------------------------------------------
// Bunch-Kaufman solve
// ---------------------------------------------------------------------------

/// Solves the system using the LDL^T factorization.
fn bunch_kaufman_solve(
    a: &[f64],
    ipiv: &[i32],
    b: &mut [f64],
    n: usize,
    nrhs: usize,
    uplo: FillMode,
) -> SolverResult<()> {
    match uplo {
        FillMode::Lower => bunch_kaufman_solve_lower(a, ipiv, b, n, nrhs),
        FillMode::Upper => bunch_kaufman_solve_lower(a, ipiv, b, n, nrhs),
        FillMode::Full => Err(SolverError::DimensionMismatch(
            "ldlt_solve: uplo must be Lower or Upper".into(),
        )),
    }
}

/// Lower-triangular solve: L*D*L^T * x = b.
fn bunch_kaufman_solve_lower(
    a: &[f64],
    ipiv: &[i32],
    b: &mut [f64],
    n: usize,
    nrhs: usize,
) -> SolverResult<()> {
    for rhs in 0..nrhs {
        let b_col = &mut b[rhs * n..(rhs + 1) * n];

        // Step 1: Apply permutations and forward substitution (L * y = P * b).
        let mut k = 0;
        while k < n {
            if ipiv[k] > 0 {
                // 1x1 pivot.
                let p = (ipiv[k] - 1) as usize;
                if p != k {
                    b_col.swap(k, p);
                }
                // Forward sub: b[i] -= L[i,k] * b[k] for i > k.
                for i in (k + 1)..n {
                    b_col[i] -= a[k * n + i] * b_col[k];
                }
                k += 1;
            } else {
                // 2x2 pivot.
                let p = ((-ipiv[k]) - 1) as usize;
                if p != k + 1 {
                    b_col.swap(k + 1, p);
                }
                for i in (k + 2)..n {
                    b_col[i] -= a[k * n + i] * b_col[k];
                    b_col[i] -= a[(k + 1) * n + i] * b_col[k + 1];
                }
                k += 2;
            }
        }

        // Step 2: Solve D * z = y.
        k = 0;
        while k < n {
            if ipiv[k] > 0 {
                // 1x1 block.
                let dkk = a[k * n + k];
                if dkk.abs() < 1e-300 {
                    return Err(SolverError::SingularMatrix);
                }
                b_col[k] /= dkk;
                k += 1;
            } else {
                // 2x2 block.
                if k + 1 >= n {
                    return Err(SolverError::InternalError(
                        "ldlt_solve: invalid 2x2 pivot at boundary".into(),
                    ));
                }
                let d11 = a[k * n + k];
                let d21 = a[k * n + (k + 1)];
                let d22 = a[(k + 1) * n + (k + 1)];
                let det = d11 * d22 - d21 * d21;
                if det.abs() < 1e-300 {
                    return Err(SolverError::SingularMatrix);
                }
                let inv_det = 1.0 / det;
                let y1 = b_col[k];
                let y2 = b_col[k + 1];
                b_col[k] = (d22 * y1 - d21 * y2) * inv_det;
                b_col[k + 1] = (-d21 * y1 + d11 * y2) * inv_det;
                k += 2;
            }
        }

        // Step 3: Backward substitution (L^T * w = z) and apply P^T.
        k = n;
        while k > 0 {
            k -= 1;
            if ipiv[k] > 0 {
                // 1x1 pivot — backward sub.
                for i in (k + 1)..n {
                    b_col[k] -= a[k * n + i] * b_col[i];
                }
                let p = (ipiv[k] - 1) as usize;
                if p != k {
                    b_col.swap(k, p);
                }
            } else if k > 0 && ipiv[k] < 0 && ipiv[k - 1] == ipiv[k] {
                // 2x2 pivot — process both rows.
                let k2 = k - 1;
                for i in (k + 1)..n {
                    b_col[k] -= a[k * n + i] * b_col[i]; // Note: this is L^T
                    b_col[k2] -= a[k2 * n + i] * b_col[i];
                }
                let p = ((-ipiv[k]) - 1) as usize;
                if p != k {
                    b_col.swap(k, p);
                }
                k = k2; // Skip past the 2x2 block.
            }
        }
    }

    Ok(())
}

fn mirror_upper_to_lower(a: &mut [f64], n: usize) {
    for col in 0..n {
        for row in 0..col {
            a[col * n + row] = a[row * n + col];
        }
    }
}

// ---------------------------------------------------------------------------
// Device buffer read/write helpers (structural)
// ---------------------------------------------------------------------------

fn read_device_to_host<T: GpuFloat>(
    buf: &DeviceBuffer<T>,
    host: &mut [f64],
    count: usize,
) -> SolverResult<()> {
    if host.len() < count {
        return Err(SolverError::DimensionMismatch(format!(
            "read_device_to_host: host buffer too small ({} < {})",
            host.len(),
            count
        )));
    }
    let mut staged = vec![T::gpu_zero(); count];
    buf.copy_to_host(&mut staged)?;
    for (dst, src) in host.iter_mut().zip(staged.iter()) {
        *dst = to_f64(*src);
    }
    Ok(())
}

fn write_host_to_device<T: GpuFloat>(
    buf: &mut DeviceBuffer<T>,
    data: &[T],
    count: usize,
) -> SolverResult<()> {
    if data.len() < count {
        return Err(SolverError::DimensionMismatch(format!(
            "write_host_to_device: source buffer too small ({} < {})",
            data.len(),
            count
        )));
    }
    buf.copy_from_host(&data[..count])?;
    Ok(())
}

fn read_device_to_host_i32(
    buf: &DeviceBuffer<i32>,
    host: &mut [i32],
    count: usize,
) -> SolverResult<()> {
    if host.len() < count {
        return Err(SolverError::DimensionMismatch(format!(
            "read_device_to_host_i32: host buffer too small ({} < {})",
            host.len(),
            count
        )));
    }
    buf.copy_to_host(&mut host[..count])?;
    Ok(())
}

fn write_host_to_device_i32(
    buf: &mut DeviceBuffer<i32>,
    data: &[i32],
    count: usize,
) -> SolverResult<()> {
    if data.len() < count {
        return Err(SolverError::DimensionMismatch(format!(
            "write_host_to_device_i32: source buffer too small ({} < {})",
            data.len(),
            count
        )));
    }
    buf.copy_from_host(&data[..count])?;
    Ok(())
}

// ---------------------------------------------------------------------------
// Tests
// ---------------------------------------------------------------------------

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn bunch_kaufman_alpha_value() {
        let expected = (1.0_f64 + 17.0_f64.sqrt()) / 8.0;
        assert!((BUNCH_KAUFMAN_ALPHA - expected).abs() < 1e-10);
    }

    #[test]
    fn column_max_offdiag_lower() {
        // 3x3 matrix (column-major):
        // [1  0  0]
        // [5  2  0]
        // [3  7  4]
        let a = [1.0, 5.0, 3.0, 0.0, 2.0, 7.0, 0.0, 0.0, 4.0];
        let (max_val, max_idx) = column_max_offdiag(&a, 3, 0, true);
        assert!((max_val - 5.0).abs() < 1e-15);
        assert_eq!(max_idx, 1);
    }

    #[test]
    fn column_max_offdiag_upper() {
        let a = [1.0, 5.0, 3.0, 0.0, 2.0, 7.0, 0.0, 0.0, 4.0];
        let (max_val, max_idx) = column_max_offdiag(&a, 3, 2, false);
        // Column 2 entries above diagonal: a[2*3+0]=0.0, a[2*3+1]=0.0
        assert!(max_val.abs() < 1e-15);
        assert_eq!(max_idx, 2); // stays at col when nothing found
    }

    #[test]
    fn swap_rows_and_cols_identity() {
        // Swapping same index should be a no-op.
        let mut a = [1.0, 0.0, 0.0, 1.0];
        swap_rows_and_cols(&mut a, 2, 0, 0);
        assert!((a[0] - 1.0).abs() < 1e-15);
        assert!((a[3] - 1.0).abs() < 1e-15);
    }

    #[test]
    fn swap_rows_and_cols_basic() {
        // 2x2 identity, swap 0 and 1.
        let mut a = [1.0, 0.0, 0.0, 1.0];
        swap_rows_and_cols(&mut a, 2, 0, 1);
        // After swap: [[1, 0], [0, 1]] becomes [[1, 0], [0, 1]] (symmetric swap).
        // Actually for identity, swapping rows and cols gives identity back.
        assert!((a[0] - 1.0).abs() < 1e-15);
        assert!((a[3] - 1.0).abs() < 1e-15);
    }

    #[test]
    fn perform_1x1_pivot_lower_basic() {
        // 2x2 matrix: [[4, 2], [2, 3]] (column-major: [4, 2, 2, 3])
        let mut a = [4.0, 2.0, 2.0, 3.0];
        perform_1x1_pivot_lower(&mut a, 2, 0);
        // L[1,0] = A[1,0] / A[0,0] = 2/4 = 0.5
        assert!((a[1] - 0.5).abs() < 1e-15);
        // A[1,1] -= L[1,0] * D[0,0] * L[1,0] = 3 - 0.5*4*0.5 = 3 - 1 = 2
        assert!((a[3] - 2.0).abs() < 1e-15);
    }

    #[test]
    fn bunch_kaufman_identity_3x3() {
        // Identity matrix should trivially factorize.
        let mut a = vec![0.0; 9];
        a[0] = 1.0;
        a[4] = 1.0;
        a[8] = 1.0;
        let mut ipiv = vec![0_i32; 3];
        let result = bunch_kaufman_lower(&mut a, 3, &mut ipiv);
        assert!(result.is_ok());
        // All pivots should be 1x1 (positive).
        assert!(ipiv[0] > 0);
        assert!(ipiv[1] > 0);
        assert!(ipiv[2] > 0);
    }

    #[test]
    fn f64_conversion_roundtrip() {
        let val = std::f64::consts::E;
        let converted: f64 = from_f64(to_f64(val));
        assert!((converted - val).abs() < 1e-15);
    }

    #[test]
    fn f32_conversion_roundtrip() {
        let val = std::f32::consts::E;
        let as_f64 = to_f64(val);
        let back: f32 = from_f64(as_f64);
        assert!((back - val).abs() < 1e-5);
    }
}