//! Comprehensive tests for linear algebra decompositions
//!
//! This module contains all tests for advanced linear algebra operations following strict TDD.
//! Tests are written before implementation to ensure proper behavior.

use crate::{Tensor, Device, TensorError, Result};
use crate::linalg::{SVDResult, QRResult, EigenResult, TriangularType, TriangularSolveOptions, 
                   LeastSquaresOptions, LeastSquaresResult, LeastSquaresMethod, LUResult, LUOptions};
use approx::assert_abs_diff_eq;

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

    #[test]
    fn test_svd_square_matrix() {
        let device = Device::cpu();
        
        // Create a simple 3x3 matrix for SVD
        // A = [[3, 2, 2], [2, 3, -2], [2, -2, 3]]
        let a = Tensor::from_data(
            vec![3.0, 2.0, 2.0, 2.0, 3.0, -2.0, 2.0, -2.0, 3.0],
            vec![3, 3],
            &device,
        ).unwrap();

        let svd_result = a.svd(true, true).unwrap();
        
        // Verify dimensions
        assert_eq!(svd_result.u.shape().dims(), &[3, 3]);
        assert_eq!(svd_result.s.shape().dims(), &[3]);
        assert_eq!(svd_result.vt.shape().dims(), &[3, 3]);

        // Verify reconstruction: A ≈ U * S * V^T
        let s_diag = Tensor::diag(&svd_result.s).unwrap();
        let reconstructed = svd_result.u.matmul(&s_diag).unwrap().matmul(&svd_result.vt).unwrap();
        
        let a_data = a.to_cpu().unwrap();
        let reconstructed_data = reconstructed.to_cpu().unwrap();
        
        for (&expected, &actual) in a_data.iter().zip(reconstructed_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }

    #[test]
    fn test_svd_rectangular_matrix() {
        let device = Device::cpu();
        
        // Test with 4x2 matrix
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0],
            vec![4, 2],
            &device,
        ).unwrap();

        let svd_result = a.svd(true, true).unwrap();
        
        // Verify dimensions
        assert_eq!(svd_result.u.shape().dims(), &[4, 4]);
        assert_eq!(svd_result.s.shape().dims(), &[2]);  // min(m,n)
        assert_eq!(svd_result.vt.shape().dims(), &[2, 2]);

        // Verify singular values are non-negative and ordered
        let s_data = svd_result.s.to_cpu().unwrap();
        assert!(s_data[0] >= s_data[1]);
        assert!(s_data[1] >= 0.0);
    }

    #[test]
    fn test_svd_thin_mode() {
        let device = Device::cpu();
        
        // Test thin SVD (reduced) mode
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
            vec![3, 2],
            &device,
        ).unwrap();

        let svd_result = a.svd(false, false).unwrap();
        
        // In thin mode, U should be 3x2 and V should be 2x2
        assert_eq!(svd_result.u.shape().dims(), &[3, 2]);
        assert_eq!(svd_result.s.shape().dims(), &[2]);
        assert_eq!(svd_result.vt.shape().dims(), &[2, 2]);
    }

    #[test]
    fn test_svd_orthogonality() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let svd_result = a.svd(true, true).unwrap();
        
        // Test U^T * U = I
        let u_transpose = svd_result.u.transpose(0, 1).unwrap();
        let u_t_u = u_transpose.matmul(&svd_result.u).unwrap();
        let identity = Tensor::eye(2, &device).unwrap();
        
        let u_t_u_data = u_t_u.to_cpu().unwrap();
        let identity_data = identity.to_cpu().unwrap();
        
        for (&expected, &actual) in identity_data.iter().zip(u_t_u_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }

    #[test]
    fn test_truncated_svd() {
        let device = Device::cpu();
        
        // Create a 5x4 matrix
        let a = Tensor::from_data(
            (0..20).map(|x| x as f32).collect(),
            vec![5, 4],
            &device,
        ).unwrap();

        // Get only top 2 singular values/vectors
        let svd_result = a.svd_truncated(2).unwrap();
        
        assert_eq!(svd_result.u.shape().dims(), &[5, 2]);
        assert_eq!(svd_result.s.shape().dims(), &[2]);
        assert_eq!(svd_result.vt.shape().dims(), &[2, 4]);
        
        // Singular values should be in descending order
        let s_data = svd_result.s.to_cpu().unwrap();
        assert!(s_data[0] >= s_data[1]);
        assert!(s_data[1] >= 0.0);
    }

    #[test]
    fn test_svd_zero_matrix() {
        let device = Device::cpu();
        
        let a = Tensor::zeros([3, 3], &device).unwrap();
        let svd_result = a.svd(true, true).unwrap();
        
        // All singular values should be zero
        let s_data = svd_result.s.to_cpu().unwrap();
        for &val in s_data.iter() {
            assert_abs_diff_eq!(val, 0.0, epsilon = 1e-8);
        }
    }

    #[test]
    fn test_svd_rank_deficient() {
        let device = Device::cpu();
        
        // Create rank-1 matrix: outer product of two vectors
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 3.0, 6.0, 9.0],
            vec![3, 3],
            &device,
        ).unwrap();

        let svd_result = a.svd(true, true).unwrap();
        let s_data = svd_result.s.to_cpu().unwrap();
        
        // Should have only one non-zero singular value
        assert!(s_data[0] > 1e-5);
        assert_abs_diff_eq!(s_data[1], 0.0, epsilon = 1e-5);
        assert_abs_diff_eq!(s_data[2], 0.0, epsilon = 1e-5);
    }
}

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

    #[test]
    fn test_qr_square_matrix() {
        let device = Device::cpu();
        
        // Create a 3x3 matrix
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0],
            vec![3, 3],
            &device,
        ).unwrap();

        let qr_result = a.qr().unwrap();
        
        // Verify dimensions
        assert_eq!(qr_result.q.shape().dims(), &[3, 3]);
        assert_eq!(qr_result.r.shape().dims(), &[3, 3]);

        // Verify reconstruction: A = Q * R
        let reconstructed = qr_result.q.matmul(&qr_result.r).unwrap();
        
        let a_data = a.to_cpu().unwrap();
        let reconstructed_data = reconstructed.to_cpu().unwrap();
        
        for (&expected, &actual) in a_data.iter().zip(reconstructed_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }

    #[test]
    fn test_qr_rectangular_matrix() {
        let device = Device::cpu();
        
        // Test with 4x3 matrix
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0],
            vec![4, 3],
            &device,
        ).unwrap();

        let qr_result = a.qr().unwrap();
        
        // Verify dimensions
        assert_eq!(qr_result.q.shape().dims(), &[4, 3]);  // Reduced Q
        assert_eq!(qr_result.r.shape().dims(), &[3, 3]);  // Square R

        // Verify R is upper triangular
        let r_data = qr_result.r.to_cpu().unwrap();
        let dims = qr_result.r.shape().dims();
        let rows = dims[0];
        let cols = dims[1];
        
        for i in 0..rows {
            for j in 0..cols {
                if i > j {
                    let idx = i * cols + j;
                    assert_abs_diff_eq!(r_data[idx], 0.0, epsilon = 1e-8);
                }
            }
        }
    }

    #[test]
    fn test_qr_orthogonality() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let qr_result = a.qr().unwrap();
        
        // Test Q^T * Q = I
        let q_transpose = qr_result.q.transpose(0, 1).unwrap();
        let q_t_q = q_transpose.matmul(&qr_result.q).unwrap();
        let identity = Tensor::eye(2, &device).unwrap();
        
        let q_t_q_data = q_t_q.to_cpu().unwrap();
        let identity_data = identity.to_cpu().unwrap();
        
        for (&expected, &actual) in identity_data.iter().zip(q_t_q_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }

    #[test]
    fn test_qr_tall_matrix() {
        let device = Device::cpu();
        
        // Test with 5x3 tall matrix
        let a = Tensor::from_data(
            (1..=15).map(|x| x as f32).collect(),
            vec![5, 3],
            &device,
        ).unwrap();

        let qr_result = a.qr().unwrap();
        
        assert_eq!(qr_result.q.shape().dims(), &[5, 3]);
        assert_eq!(qr_result.r.shape().dims(), &[3, 3]);
    }
}

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

    #[test]
    fn test_cholesky_simple_matrix() {
        let device = Device::cpu();
        
        // Create a simple positive definite matrix
        // A = [[4, 2], [2, 3]]
        let a = Tensor::from_data(
            vec![4.0, 2.0, 2.0, 3.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let l = a.cholesky(false).unwrap();
        
        // Verify dimensions
        assert_eq!(l.shape().dims(), &[2, 2]);

        // Verify L is lower triangular
        let l_data = l.to_cpu().unwrap();
        assert_abs_diff_eq!(l_data[1], 0.0, epsilon = 1e-8); // L[0,1] should be 0

        // Verify reconstruction: A = L * L^T
        let l_transpose = l.transpose(0, 1).unwrap();
        let reconstructed = l.matmul(&l_transpose).unwrap();
        
        let a_data = a.to_cpu().unwrap();
        let reconstructed_data = reconstructed.to_cpu().unwrap();
        
        for (&expected, &actual) in a_data.iter().zip(reconstructed_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }

    #[test]
    fn test_cholesky_3x3_matrix() {
        let device = Device::cpu();
        
        // Create a 3x3 positive definite matrix
        let a = Tensor::from_data(
            vec![4.0, 2.0, 1.0, 2.0, 3.0, 0.5, 1.0, 0.5, 2.0],
            vec![3, 3],
            &device,
        ).unwrap();

        let l = a.cholesky(false).unwrap();
        
        // Verify L is lower triangular
        let l_data = l.to_cpu().unwrap();
        assert_abs_diff_eq!(l_data[1], 0.0, epsilon = 1e-8); // L[0,1]
        assert_abs_diff_eq!(l_data[2], 0.0, epsilon = 1e-8); // L[0,2]
        assert_abs_diff_eq!(l_data[5], 0.0, epsilon = 1e-8); // L[1,2]
    }

    #[test]
    fn test_cholesky_upper_triangular() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![4.0, 2.0, 2.0, 3.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let u = a.cholesky(true).unwrap();
        
        // Verify U is upper triangular
        let u_data = u.to_cpu().unwrap();
        assert_abs_diff_eq!(u_data[2], 0.0, epsilon = 1e-8); // U[1,0] should be 0
    }

    #[test]
    fn test_cholesky_not_positive_definite() {
        let device = Device::cpu();
        
        // Create a matrix that is NOT positive definite
        let a = Tensor::from_data(
            vec![1.0, 2.0, 2.0, 1.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let result = a.cholesky(false);
        assert!(result.is_err(), "Cholesky should fail for non-positive definite matrix");
    }

    #[test]
    fn test_cholesky_identity_matrix() {
        let device = Device::cpu();
        
        let a = Tensor::eye(3, &device).unwrap();
        let l = a.cholesky(false).unwrap();
        
        // Cholesky of identity should be identity
        let l_data = l.to_cpu().unwrap();
        let identity_data = a.to_cpu().unwrap();
        
        for (&expected, &actual) in identity_data.iter().zip(l_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-8);
        }
    }
}

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

    #[test]
    fn test_eigen_symmetric_matrix() {
        let device = Device::cpu();
        
        // Create a simple symmetric matrix
        let a = Tensor::from_data(
            vec![3.0, 1.0, 1.0, 3.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let eigen_result = a.symeig(true).unwrap();
        
        // Verify dimensions
        assert_eq!(eigen_result.eigenvalues.shape().dims(), &[2]);
        assert_eq!(eigen_result.eigenvectors.shape().dims(), &[2, 2]);

        // Verify eigenvalues are in ascending order
        let vals = eigen_result.eigenvalues.to_cpu().unwrap();
        assert!(vals[0] <= vals[1]);

        // Verify eigenvector normalization
        let vecs_data = eigen_result.eigenvectors.to_cpu().unwrap();
        let v1_norm = (vecs_data[0] * vecs_data[0] + vecs_data[2] * vecs_data[2]).sqrt();
        let v2_norm = (vecs_data[1] * vecs_data[1] + vecs_data[3] * vecs_data[3]).sqrt();
        
        assert_abs_diff_eq!(v1_norm, 1.0, epsilon = 1e-6);
        assert_abs_diff_eq!(v2_norm, 1.0, epsilon = 1e-6);
    }

    #[test]
    fn test_eigen_3x3_symmetric() {
        let device = Device::cpu();
        
        // Create a 3x3 symmetric matrix
        let a = Tensor::from_data(
            vec![4.0, 1.0, 2.0, 1.0, 3.0, 0.0, 2.0, 0.0, 2.0],
            vec![3, 3],
            &device,
        ).unwrap();

        let eigen_result = a.symeig(true).unwrap();
        
        assert_eq!(eigen_result.eigenvalues.shape().dims(), &[3]);
        assert_eq!(eigen_result.eigenvectors.shape().dims(), &[3, 3]);
        
        // Check that eigenvalues are sorted
        let vals = eigen_result.eigenvalues.to_cpu().unwrap();
        assert!(vals[0] <= vals[1]);
        assert!(vals[1] <= vals[2]);
    }

    #[test]
    fn test_eigen_diagonal_matrix() {
        let device = Device::cpu();
        
        // Diagonal matrix - eigenvalues should be the diagonal elements
        let diag_vals = vec![5.0, 3.0, 7.0];
        let a = Tensor::diag(&Tensor::from_data(diag_vals.clone(), vec![3], &device).unwrap()).unwrap();
        
        let eigen_result = a.symeig(true).unwrap();
        let computed_vals = eigen_result.eigenvalues.to_cpu().unwrap();
        
        let mut sorted_diag = diag_vals.clone();
        sorted_diag.sort_by(|a, b| a.partial_cmp(b).unwrap());
        
        for (&expected, &actual) in sorted_diag.iter().zip(computed_vals.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    fn test_eigen_orthogonality() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![2.0, 1.0, 1.0, 2.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let eigen_result = a.symeig(true).unwrap();
        
        // Test V^T * V = I (eigenvectors are orthogonal)
        let v_transpose = eigen_result.eigenvectors.transpose(0, 1).unwrap();
        let v_t_v = v_transpose.matmul(&eigen_result.eigenvectors).unwrap();
        let identity = Tensor::eye(2, &device).unwrap();
        
        let v_t_v_data = v_t_v.to_cpu().unwrap();
        let identity_data = identity.to_cpu().unwrap();
        
        for (&expected, &actual) in identity_data.iter().zip(v_t_v_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }

    #[test]
    fn test_eigen_reconstruction() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![3.0, 2.0, 2.0, 3.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let eigen_result = a.symeig(true).unwrap();
        
        // Verify A = V * Λ * V^T
        let lambda_diag = Tensor::diag(&eigen_result.eigenvalues).unwrap();
        let v_transpose = eigen_result.eigenvectors.transpose(0, 1).unwrap();
        let reconstructed = eigen_result.eigenvectors
            .matmul(&lambda_diag).unwrap()
            .matmul(&v_transpose).unwrap();
        
        let a_data = a.to_cpu().unwrap();
        let reconstructed_data = reconstructed.to_cpu().unwrap();
        
        for (&expected, &actual) in a_data.iter().zip(reconstructed_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }
}

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

    #[test]
    fn test_pinv_square_invertible() {
        let device = Device::cpu();
        
        // Test with invertible 2x2 matrix
        let a = Tensor::from_data(
            vec![2.0, 1.0, 1.0, 2.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let pinv_a = a.pinv(1e-8).unwrap();
        
        // Verify dimensions
        assert_eq!(pinv_a.shape().dims(), &[2, 2]);
        
        // Verify A * A^+ = I (approximately)
        let identity_approx = a.matmul(&pinv_a).unwrap();
        let identity = Tensor::eye(2, &device).unwrap();
        
        let identity_data = identity.to_cpu().unwrap();
        let computed_data = identity_approx.to_cpu().unwrap();
        
        for (&expected, &actual) in identity_data.iter().zip(computed_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }

    #[test]
    fn test_pinv_rectangular_full_rank() {
        let device = Device::cpu();
        
        // Test with 3x2 full rank matrix
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
            vec![3, 2],
            &device,
        ).unwrap();

        let pinv_a = a.pinv(1e-8).unwrap();
        
        // Dimensions should be 2x3
        assert_eq!(pinv_a.shape().dims(), &[2, 3]);
        
        // For full rank matrices: A^+ * A = I
        let should_be_identity = pinv_a.matmul(&a).unwrap();
        let identity = Tensor::eye(2, &device).unwrap();
        
        let identity_data = identity.to_cpu().unwrap();
        let computed_data = should_be_identity.to_cpu().unwrap();
        
        for (&expected, &actual) in identity_data.iter().zip(computed_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-4);
        }
    }

    #[test]
    fn test_pinv_rank_deficient() {
        let device = Device::cpu();
        
        // Create rank-deficient matrix (rank 1)
        let a = Tensor::from_data(
            vec![1.0, 2.0, 2.0, 4.0, 3.0, 6.0],
            vec![3, 2],
            &device,
        ).unwrap();

        let pinv_a = a.pinv(1e-6).unwrap();
        
        // Should still compute pseudo-inverse
        assert_eq!(pinv_a.shape().dims(), &[2, 3]);
        
        // Test Moore-Penrose conditions: A * A^+ * A = A
        let temp = a.matmul(&pinv_a).unwrap();
        let reconstructed = temp.matmul(&a).unwrap();
        
        let original_data = a.to_cpu().unwrap();
        let reconstructed_data = reconstructed.to_cpu().unwrap();
        
        for (&expected, &actual) in original_data.iter().zip(reconstructed_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-4);
        }
    }

    #[test]
    fn test_pinv_zero_matrix() {
        let device = Device::cpu();
        
        let a = Tensor::zeros([2, 3], &device).unwrap();
        let pinv_a = a.pinv(1e-8).unwrap();
        
        // Pseudo-inverse of zero matrix should be zero
        let pinv_data = pinv_a.to_cpu().unwrap();
        for &val in pinv_data.iter() {
            assert_abs_diff_eq!(val, 0.0, epsilon = 1e-8);
        }
    }
}

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

    #[test]
    fn test_condition_number() {
        let device = Device::cpu();
        
        // Well-conditioned matrix
        let a = Tensor::eye(3, &device).unwrap();
        let cond = a.condition_number().unwrap();
        
        // Condition number of identity should be 1
        assert_abs_diff_eq!(cond, 1.0, epsilon = 1e-6);
    }

    #[test]
    fn test_matrix_rank() {
        let device = Device::cpu();
        
        // Full rank 2x2 matrix
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let rank = a.matrix_rank(1e-8).unwrap();
        assert_eq!(rank, 2);
        
        // Rank deficient matrix
        let b = Tensor::from_data(
            vec![1.0, 2.0, 2.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let rank_b = b.matrix_rank(1e-8).unwrap();
        assert_eq!(rank_b, 1);
    }

    #[test]
    fn test_matrix_norm() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![3.0, 4.0, 0.0, 5.0],
            vec![2, 2],
            &device,
        ).unwrap();

        // Frobenius norm
        let norm_fro = a.matrix_norm("fro").unwrap();
        let expected_fro = (9.0 + 16.0 + 0.0 + 25.0_f32).sqrt();
        assert_abs_diff_eq!(norm_fro, expected_fro, epsilon = 1e-6);
    }

    #[test]
    fn test_determinant() {
        let device = Device::cpu();
        
        // 2x2 matrix with known determinant
        let a = Tensor::from_data(
            vec![3.0, 2.0, 1.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let det = a.det().unwrap();
        let expected_det = 3.0 * 4.0 - 2.0 * 1.0; // 12 - 2 = 10
        assert_abs_diff_eq!(det, expected_det, epsilon = 1e-6);
    }

    #[test]
    fn test_trace() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0],
            vec![3, 3],
            &device,
        ).unwrap();

        let trace = a.trace().unwrap();
        let expected_trace = 1.0 + 5.0 + 9.0; // diagonal sum
        assert_abs_diff_eq!(trace, expected_trace, epsilon = 1e-6);
    }

    #[test]
    fn test_diag_creation() {
        let device = Device::cpu();
        
        let diag_values = Tensor::from_data(vec![1.0, 2.0, 3.0], vec![3], &device).unwrap();
        let diag_matrix = Tensor::diag(&diag_values).unwrap();
        
        assert_eq!(diag_matrix.shape().dims(), &[3, 3]);
        
        let data = diag_matrix.to_cpu().unwrap();
        // Check diagonal elements
        assert_abs_diff_eq!(data[0], 1.0, epsilon = 1e-8); // [0,0]
        assert_abs_diff_eq!(data[4], 2.0, epsilon = 1e-8); // [1,1]
        assert_abs_diff_eq!(data[8], 3.0, epsilon = 1e-8); // [2,2]
        
        // Check off-diagonal elements are zero
        assert_abs_diff_eq!(data[1], 0.0, epsilon = 1e-8); // [0,1]
        assert_abs_diff_eq!(data[2], 0.0, epsilon = 1e-8); // [0,2]
        assert_abs_diff_eq!(data[3], 0.0, epsilon = 1e-8); // [1,0]
    }

    #[test]
    fn test_eye_creation() {
        let device = Device::cpu();
        
        let identity = Tensor::eye(3, &device).unwrap();
        assert_eq!(identity.shape().dims(), &[3, 3]);
        
        let data = identity.to_cpu().unwrap();
        for i in 0..3 {
            for j in 0..3 {
                let idx = i * 3 + j;
                if i == j {
                    assert_abs_diff_eq!(data[idx], 1.0, epsilon = 1e-8);
                } else {
                    assert_abs_diff_eq!(data[idx], 0.0, epsilon = 1e-8);
                }
            }
        }
    }
}

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

    #[test]
    fn test_svd_invalid_dimensions() {
        let device = Device::cpu();
        let empty_tensor = Tensor::zeros([0, 3], &device).unwrap();
        let result = empty_tensor.svd(true, true);
        assert!(result.is_err());
    }

    #[test]
    fn test_qr_invalid_dimensions() {
        let device = Device::cpu();
        let invalid_tensor = Tensor::zeros([2, 0], &device).unwrap();
        let result = invalid_tensor.qr();
        assert!(result.is_err());
    }

    #[test]
    fn test_cholesky_non_square() {
        let device = Device::cpu();
        let non_square = Tensor::ones([3, 2], &device).unwrap();
        let result = non_square.cholesky(false);
        assert!(result.is_err());
    }

    #[test]
    fn test_symeig_non_square() {
        let device = Device::cpu();
        let non_square = Tensor::ones([2, 3], &device).unwrap();
        let result = non_square.symeig(true);
        assert!(result.is_err());
    }
}

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

    #[test]
    fn test_svd_numerical_precision() {
        let device = Device::cpu();
        
        // Create a matrix with very small and very large values
        let a = Tensor::from_data(
            vec![1e6, 1e-6, 1e-6, 1e6],
            vec![2, 2],
            &device,
        ).unwrap();

        let svd_result = a.svd(true, true).unwrap();
        
        // Should not fail and should produce reasonable results
        let s_data = svd_result.s.to_cpu().unwrap();
        assert!(s_data[0] > s_data[1]);
        assert!(s_data[1] >= 0.0);
        assert!(s_data[0].is_finite());
        assert!(s_data[1].is_finite());
    }

    #[test]
    fn test_qr_near_singular() {
        let device = Device::cpu();
        
        // Create a nearly singular matrix
        let a = Tensor::from_data(
            vec![1.0, 1.0, 1.0, 1.0 + 1e-10],
            vec![2, 2],
            &device,
        ).unwrap();

        let qr_result = a.qr().unwrap();
        
        // Should complete without numerical issues
        let q_data = qr_result.q.to_cpu().unwrap();
        let r_data = qr_result.r.to_cpu().unwrap();
        
        // Check that results are finite
        for &val in q_data.iter().chain(r_data.iter()) {
            assert!(val.is_finite(), "QR decomposition produced non-finite values");
        }
    }

    #[test]
    fn test_cholesky_barely_positive_definite() {
        let device = Device::cpu();
        
        // Create a barely positive definite matrix
        let a = Tensor::from_data(
            vec![1.0 + 1e-8, 0.0, 0.0, 1e-8],
            vec![2, 2],
            &device,
        ).unwrap();

        let l = a.cholesky(false).unwrap();
        
        // Should succeed and produce finite results
        let l_data = l.to_cpu().unwrap();
        for &val in l_data.iter() {
            assert!(val.is_finite(), "Cholesky produced non-finite values");
        }
    }
}

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

    #[test]
    fn test_svd_to_pinv_consistency() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0],
            vec![3, 2],
            &device,
        ).unwrap();

        // Compute pseudo-inverse using SVD directly
        let svd_result = a.svd(true, true).unwrap();
        let s_inv = svd_result.s.reciprocal_safe(1e-8).unwrap();
        let s_inv_diag = Tensor::diag(&s_inv).unwrap();
        let vt_transpose = svd_result.vt.transpose(0, 1).unwrap();
        let u_transpose = svd_result.u.transpose(0, 1).unwrap();
        
        let pinv_from_svd = vt_transpose
            .matmul(&s_inv_diag).unwrap()
            .matmul(&u_transpose).unwrap();
        
        // Compare with direct pinv method
        let pinv_direct = a.pinv(1e-8).unwrap();
        
        let svd_data = pinv_from_svd.to_cpu().unwrap();
        let direct_data = pinv_direct.to_cpu().unwrap();
        
        for (&expected, &actual) in svd_data.iter().zip(direct_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-4);
        }
    }

    #[test]
    fn test_qr_to_least_squares() {
        let device = Device::cpu();
        
        // Solve Ax = b using QR decomposition
        let a = Tensor::from_data(
            vec![1.0, 1.0, 1.0, 2.0, 2.0, 1.0],
            vec![3, 2],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(
            vec![6.0, 8.0, 7.0],
            vec![3, 1],
            &device,
        ).unwrap();

        let qr_result = a.qr().unwrap();
        
        // Solve Rx = Q^T b
        let qt = qr_result.q.transpose(0, 1).unwrap();
        let qtb = qt.matmul(&b).unwrap();
        
        // This would require solving upper triangular system
        // For now, just verify the QR decomposition worked
        let reconstructed = qr_result.q.matmul(&qr_result.r).unwrap();
        
        let a_data = a.to_cpu().unwrap();
        let reconstructed_data = reconstructed.to_cpu().unwrap();
        
        for (&expected, &actual) in a_data.iter().zip(reconstructed_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-5);
        }
    }

    #[test]
    fn test_eigen_to_matrix_powers() {
        let device = Device::cpu();
        
        // Test that eigendecomposition can be used for matrix powers
        let a = Tensor::from_data(
            vec![3.0, 1.0, 1.0, 3.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let eigen_result = a.symeig(true).unwrap();
        
        // Compute A^2 = V * Λ^2 * V^T
        let lambda_squared = eigen_result.eigenvalues.mul(&eigen_result.eigenvalues).unwrap();
        let lambda_squared_diag = Tensor::diag(&lambda_squared).unwrap();
        let vt = eigen_result.eigenvectors.transpose(0, 1).unwrap();
        
        let a_squared_from_eigen = eigen_result.eigenvectors
            .matmul(&lambda_squared_diag).unwrap()
            .matmul(&vt).unwrap();
        
        // Compare with direct computation
        let a_squared_direct = a.matmul(&a).unwrap();
        
        let eigen_data = a_squared_from_eigen.to_cpu().unwrap();
        let direct_data = a_squared_direct.to_cpu().unwrap();
        
        for (&expected, &actual) in direct_data.iter().zip(eigen_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-4);
        }
    }
}

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

    /// Helper function to create default lower triangular solve options
    fn default_lower_options() -> TriangularSolveOptions {
        TriangularSolveOptions {
            triangular_type: TriangularType::Lower,
            transpose: false,
            unit_diagonal: false,
            validate_triangular: true,
            strict_validation: false,
        }
    }

    /// Helper function to create default upper triangular solve options
    fn default_upper_options() -> TriangularSolveOptions {
        TriangularSolveOptions {
            triangular_type: TriangularType::Upper,
            transpose: false,
            unit_diagonal: false,
            validate_triangular: true,
            strict_validation: false,
        }
    }

    #[test]
    fn test_solve_lower_triangular_simple_2x2() {
        let device = Device::cpu();
        
        // Create a simple 2x2 lower triangular matrix:
        // L = [[2, 0], [3, 4]]
        let l = Tensor::from_data(
            vec![2.0, 0.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        // Right-hand side: b = [4, 18]
        let b = Tensor::from_data(vec![4.0, 18.0], vec![2], &device).unwrap();
        
        // Expected solution: x = [2, 3] (since L*x = b)
        // Forward substitution: 
        // x[0] = b[0] / L[0,0] = 4 / 2 = 2
        // x[1] = (b[1] - L[1,0]*x[0]) / L[1,1] = (18 - 3*2) / 4 = 12/4 = 3
        let options = default_lower_options();
        let result = l.solve_triangular(&b, &options).unwrap();
        
        assert_eq!(result.shape().dims(), &[2]);
        let result_data = result.to_cpu().unwrap();
        assert_abs_diff_eq!(result_data[0], 2.0, epsilon = 1e-6);
        assert_abs_diff_eq!(result_data[1], 3.0, epsilon = 1e-6);
    }

    #[test]
    fn test_solve_upper_triangular_simple_2x2() {
        let device = Device::cpu();
        
        // Create a simple 2x2 upper triangular matrix:
        // U = [[2, 3], [0, 4]]
        let u = Tensor::from_data(
            vec![2.0, 3.0, 0.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        // Right-hand side: b = [13, 12]
        let b = Tensor::from_data(vec![13.0, 12.0], vec![2], &device).unwrap();
        
        // Expected solution: x = [2, 3] (since U*x = b)
        // Back substitution:
        // x[1] = b[1] / U[1,1] = 12 / 4 = 3
        // x[0] = (b[0] - U[0,1]*x[1]) / U[0,0] = (13 - 3*3) / 2 = 4/2 = 2
        let options = default_upper_options();
        let result = u.solve_triangular(&b, &options).unwrap();
        
        assert_eq!(result.shape().dims(), &[2]);
        let result_data = result.to_cpu().unwrap();
        assert_abs_diff_eq!(result_data[0], 2.0, epsilon = 1e-6);
        assert_abs_diff_eq!(result_data[1], 3.0, epsilon = 1e-6);
    }

    #[test]
    fn test_solve_lower_triangular_3x3() {
        let device = Device::cpu();
        
        // Create 3x3 lower triangular matrix:
        // L = [[1, 0, 0], [2, 3, 0], [4, 5, 6]]
        let l = Tensor::from_data(
            vec![1.0, 0.0, 0.0, 2.0, 3.0, 0.0, 4.0, 5.0, 6.0],
            vec![3, 3],
            &device,
        ).unwrap();
        
        // Right-hand side: b = [1, 8, 32]
        let b = Tensor::from_data(vec![1.0, 8.0, 32.0], vec![3], &device).unwrap();
        
        // Expected solution: x = [1, 2, 3]
        // Forward substitution:
        // x[0] = 1/1 = 1
        // x[1] = (8 - 2*1)/3 = 6/3 = 2  
        // x[2] = (32 - 4*1 - 5*2)/6 = (32-4-10)/6 = 18/6 = 3
        let options = default_lower_options();
        let result = l.solve_triangular(&b, &options).unwrap();
        
        assert_eq!(result.shape().dims(), &[3]);
        let result_data = result.to_cpu().unwrap();
        assert_abs_diff_eq!(result_data[0], 1.0, epsilon = 1e-6);
        assert_abs_diff_eq!(result_data[1], 2.0, epsilon = 1e-6);
        assert_abs_diff_eq!(result_data[2], 3.0, epsilon = 1e-6);
    }

    #[test]
    fn test_solve_upper_triangular_3x3() {
        let device = Device::cpu();
        
        // Create 3x3 upper triangular matrix:
        // U = [[6, 5, 4], [0, 3, 2], [0, 0, 1]]
        let u = Tensor::from_data(
            vec![6.0, 5.0, 4.0, 0.0, 3.0, 2.0, 0.0, 0.0, 1.0],
            vec![3, 3],
            &device,
        ).unwrap();
        
        // Right-hand side: b = [32, 8, 3]
        let b = Tensor::from_data(vec![32.0, 8.0, 3.0], vec![3], &device).unwrap();
        
        // Expected solution: x = [1, 2, 3]
        // Back substitution:
        // x[2] = 3/1 = 3
        // x[1] = (8 - 2*3)/3 = 2/3 = 2
        // x[0] = (32 - 5*2 - 4*3)/6 = (32-10-12)/6 = 10/6 ≈ 1.67... wait let me recalculate
        // Actually: x = [1, 2, 3] means 6*1 + 5*2 + 4*3 = 6+10+12 = 28, not 32
        // Let me fix: if x=[1,2,3] then b should be [28, 8, 3]
        let b_correct = Tensor::from_data(vec![28.0, 8.0, 3.0], vec![3], &device).unwrap();
        
        let options = default_upper_options();
        let result = u.solve_triangular(&b_correct, &options).unwrap();
        
        assert_eq!(result.shape().dims(), &[3]);
        let result_data = result.to_cpu().unwrap();
        assert_abs_diff_eq!(result_data[0], 1.0, epsilon = 1e-6);
        assert_abs_diff_eq!(result_data[1], 2.0, epsilon = 1e-6);
        assert_abs_diff_eq!(result_data[2], 3.0, epsilon = 1e-6);
    }

    #[test]
    fn test_solve_triangular_transpose_lower() {
        let device = Device::cpu();
        
        // Solve L^T x = b where L is lower triangular
        // L = [[2, 0], [3, 4]]
        // L^T = [[2, 3], [0, 4]] (which is upper triangular)
        let l = Tensor::from_data(
            vec![2.0, 0.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        // For L^T x = b where L^T = [[2,3],[0,4]], if x=[2,3], then b=[2*2+3*3, 0*2+4*3]=[13,12]
        let b = Tensor::from_data(vec![13.0, 12.0], vec![2], &device).unwrap();
        
        let mut options = default_lower_options();
        options.transpose = true; // Solve L^T x = b
        let result = l.solve_triangular(&b, &options).unwrap();
        
        assert_eq!(result.shape().dims(), &[2]);
        let result_data = result.to_cpu().unwrap();
        assert_abs_diff_eq!(result_data[0], 2.0, epsilon = 1e-6);
        assert_abs_diff_eq!(result_data[1], 3.0, epsilon = 1e-6);
    }

    #[test]
    fn test_solve_triangular_unit_diagonal_lower() {
        let device = Device::cpu();
        
        // Unit lower triangular matrix (diagonal assumed to be 1):
        // L = [[1, 0], [3, 1]] but stored as [[_, 0], [3, _]] with unit_diagonal=true
        let l = Tensor::from_data(
            vec![999.0, 0.0, 3.0, 999.0], // Diagonal values should be ignored
            vec![2, 2],
            &device,
        ).unwrap();
        
        // If true L = [[1,0],[3,1]] and x=[2,3], then b = [1*2+0*3, 3*2+1*3] = [2, 9]
        let b = Tensor::from_data(vec![2.0, 9.0], vec![2], &device).unwrap();
        
        let mut options = default_lower_options();
        options.unit_diagonal = true;
        let result = l.solve_triangular(&b, &options).unwrap();
        
        assert_eq!(result.shape().dims(), &[2]);
        let result_data = result.to_cpu().unwrap();
        assert_abs_diff_eq!(result_data[0], 2.0, epsilon = 1e-6);
        assert_abs_diff_eq!(result_data[1], 3.0, epsilon = 1e-6);
    }

    #[test]
    fn test_solve_triangular_multiple_rhs() {
        let device = Device::cpu();
        
        // Test with multiple right-hand sides (batch solving)
        // L = [[2, 0], [3, 4]]
        let l = Tensor::from_data(
            vec![2.0, 0.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        // Multiple RHS: B = [[4, 8], [18, 24]] (2x2 matrix)
        // First column: b1=[4,18], solution x1=[2,3]
        // Second column: b2=[8,24], solution x2=[4,3] since (8-3*4)/4=3 and 8/2=4
        let b = Tensor::from_data(
            vec![4.0, 8.0, 18.0, 24.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let options = default_lower_options();
        let result = l.solve_triangular(&b, &options).unwrap();
        
        assert_eq!(result.shape().dims(), &[2, 2]);
        let result_data = result.to_cpu().unwrap();
        
        // First solution: x1 = [2, 3]
        assert_abs_diff_eq!(result_data[0], 2.0, epsilon = 1e-6); // x1[0]
        assert_abs_diff_eq!(result_data[2], 3.0, epsilon = 1e-6); // x1[1]
        
        // Second solution: x2 = [4, 3]
        assert_abs_diff_eq!(result_data[1], 4.0, epsilon = 1e-6); // x2[0] 
        assert_abs_diff_eq!(result_data[3], 3.0, epsilon = 1e-6); // x2[1]
    }

    #[test]
    fn test_solve_triangular_singular_matrix() {
        let device = Device::cpu();
        
        // Create singular matrix (zero diagonal element)
        let l = Tensor::from_data(
            vec![2.0, 0.0, 3.0, 0.0], // Zero at L[1,1]
            vec![2, 2],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(vec![4.0, 18.0], vec![2], &device).unwrap();
        
        let options = default_lower_options();
        let result = l.solve_triangular(&b, &options);
        assert!(result.is_err(), "Should fail for singular matrix");
    }

    #[test]
    fn test_solve_triangular_non_square_matrix() {
        let device = Device::cpu();
        
        // Non-square matrix should fail
        let non_square = Tensor::from_data(
            vec![1.0, 0.0, 2.0, 3.0, 0.0, 4.0],
            vec![2, 3],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(vec![1.0, 2.0], vec![2], &device).unwrap();
        
        let options = default_lower_options();
        let result = non_square.solve_triangular(&b, &options);
        assert!(result.is_err(), "Should fail for non-square matrix");
    }

    #[test]
    fn test_solve_triangular_dimension_mismatch() {
        let device = Device::cpu();
        
        // Matrix and RHS dimensions don't match
        let l = Tensor::from_data(
            vec![2.0, 0.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let b_wrong_size = Tensor::from_data(vec![4.0, 18.0, 9.0], vec![3], &device).unwrap();
        
        let options = default_lower_options();
        let result = l.solve_triangular(&b_wrong_size, &options);
        assert!(result.is_err(), "Should fail for dimension mismatch");
    }

    #[test]
    fn test_solve_triangular_non_triangular_strict() {
        let device = Device::cpu();
        
        // Matrix that's not triangular (has non-zero above diagonal)
        let not_triangular = Tensor::from_data(
            vec![2.0, 1.0, 3.0, 4.0], // 1.0 at [0,1] violates lower triangular
            vec![2, 2],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(vec![4.0, 18.0], vec![2], &device).unwrap();
        
        let mut options = default_lower_options();
        options.strict_validation = true;
        let result = not_triangular.solve_triangular(&b, &options);
        assert!(result.is_err(), "Should fail strict triangular validation");
    }

    #[test]
    fn test_solve_triangular_non_triangular_relaxed() {
        let device = Device::cpu();
        
        // Matrix with very small off-triangular elements (within tolerance)
        let nearly_triangular = Tensor::from_data(
            vec![2.0, 1e-10, 3.0, 4.0], // Very small value at [0,1]
            vec![2, 2],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(vec![4.0, 18.0], vec![2], &device).unwrap();
        
        let mut options = default_lower_options();
        options.strict_validation = false; // Relaxed validation
        let result = nearly_triangular.solve_triangular(&b, &options);
        assert!(result.is_ok(), "Should pass relaxed triangular validation");
    }

    #[test]  
    fn test_solve_triangular_no_validation() {
        let device = Device::cpu();
        
        // Matrix that's not triangular, but we skip validation
        let not_triangular = Tensor::from_data(
            vec![2.0, 1.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(vec![4.0, 18.0], vec![2], &device).unwrap();
        
        let mut options = default_lower_options();
        options.validate_triangular = false; // Skip validation entirely
        let result = not_triangular.solve_triangular(&b, &options);
        assert!(result.is_ok(), "Should succeed when validation is disabled");
    }

    #[test]
    fn test_solve_triangular_identity_matrix() {
        let device = Device::cpu();
        
        // Identity matrix (both upper and lower triangular)
        let identity = Tensor::eye(3, &device).unwrap();
        let b = Tensor::from_data(vec![1.0, 2.0, 3.0], vec![3], &device).unwrap();
        
        // Solution should be the same as b
        let options = default_lower_options();
        let result = identity.solve_triangular(&b, &options).unwrap();
        
        let result_data = result.to_cpu().unwrap();
        let b_data = b.to_cpu().unwrap();
        for (&expected, &actual) in b_data.iter().zip(result_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    fn test_solve_triangular_numerical_stability() {
        let device = Device::cpu();
        
        // Test with very small diagonal elements (but not zero)
        let l = Tensor::from_data(
            vec![1e-8, 0.0, 1.0, 1e-8],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(vec![1e-8, 1.0], vec![2], &device).unwrap();
        
        let options = default_lower_options();
        let result = l.solve_triangular(&b, &options).unwrap();
        
        // Should produce finite results despite small diagonals
        let result_data = result.to_cpu().unwrap();
        for &val in result_data.iter() {
            assert!(val.is_finite(), "Result should be finite for small diagonal elements");
        }
    }

    #[test]
    fn test_solve_triangular_zero_rhs() {
        let device = Device::cpu();
        
        let l = Tensor::from_data(
            vec![2.0, 0.0, 3.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        // Zero right-hand side should give zero solution
        let b_zero = Tensor::zeros([2], &device).unwrap();
        
        let options = default_lower_options();
        let result = l.solve_triangular(&b_zero, &options).unwrap();
        
        let result_data = result.to_cpu().unwrap();
        for &val in result_data.iter() {
            assert_abs_diff_eq!(val, 0.0, epsilon = 1e-10);
        }
    }
}

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

    // Helper function to create test matrix with known solution
    fn create_overdetermined_test_system() -> (Tensor, Tensor, Tensor) {
        let device = Device::cpu();
        
        // Create 4x3 overdetermined system
        // A x = b where x = [1, 2, 3] is exact solution for first 3 rows
        let a = Tensor::from_data(
            vec![
                1.0, 2.0, 3.0,  // Row 1: 1*1 + 2*2 + 3*3 = 14
                4.0, 5.0, 6.0,  // Row 2: 4*1 + 5*2 + 6*3 = 32
                7.0, 8.0, 9.0,  // Row 3: 7*1 + 8*2 + 9*3 = 50
                1.0, 1.0, 1.0,  // Row 4: 1*1 + 1*2 + 1*3 = 6
            ],
            vec![4, 3],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(
            vec![14.0, 32.0, 50.0, 6.0],
            vec![4],
            &device,
        ).unwrap();
        
        let x_true = Tensor::from_data(
            vec![1.0, 2.0, 3.0],
            vec![3],
            &device,
        ).unwrap();
        
        (a, b, x_true)
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_least_squares_overdetermined_qr() {
        let (a, b, x_true) = create_overdetermined_test_system();
        
        // Use QR method for overdetermined system
        let result = a.least_squares(&b, None).unwrap();
        
        // Check that solution is close to expected
        let x_data = result.solution.to_cpu().unwrap();
        let x_true_data = x_true.to_cpu().unwrap();
        
        for (&expected, &actual) in x_true_data.iter().zip(x_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-4);
        }
        
        // Check residual properties
        assert!(result.residual_norm < 1e-10);
        assert!(result.condition_number.is_finite());
        assert_eq!(result.rank, 3);
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_least_squares_underdetermined_svd() {
        let device = Device::cpu();
        
        // Create 2x3 underdetermined system (more unknowns than equations)
        // Multiple solutions exist, we want minimum norm solution
        let a = Tensor::from_data(
            vec![
                1.0, 1.0, 1.0,  // x + y + z = 6
                2.0, 1.0, 0.0,  // 2x + y = 5
            ],
            vec![2, 3],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(vec![6.0, 5.0], vec![2], &device).unwrap();
        
        // Force SVD method for underdetermined system
        let mut options = LeastSquaresOptions::default();
        options.method = LeastSquaresMethod::SVD;
        
        let result = a.least_squares(&b, Some(&options)).unwrap();
        
        // Verify solution satisfies Ax = b
        let ax = a.matmul(&result.solution).unwrap();
        let ax_data = ax.to_cpu().unwrap();
        let b_data = b.to_cpu().unwrap();
        
        for (&expected, &actual) in b_data.iter().zip(ax_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
        
        // Should find minimum norm solution
        assert!(result.residual_norm < 1e-10);
        assert_eq!(result.rank, 2);
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_least_squares_weighted() {
        let (a, b, x_true) = create_overdetermined_test_system();
        let device = Device::cpu();
        
        // Create weight matrix (diagonal)
        let weights = Tensor::from_data(
            vec![1.0, 2.0, 0.5, 1.5], // Different weights for each equation
            vec![4],
            &device,
        ).unwrap();
        
        let mut options = LeastSquaresOptions::default();
        options.weights = Some(weights);
        
        let result = a.least_squares(&b, Some(&options)).unwrap();
        
        // Solution should be influenced by weights
        let x_data = result.solution.to_cpu().unwrap();
        assert!(x_data.len() == 3);
        
        // Weighted residual should be minimized
        assert!(result.residual_norm.is_finite());
        assert!(result.condition_number.is_finite());
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_least_squares_tikhonov_regularization() {
        let (a, b, _) = create_overdetermined_test_system();
        
        // Add Tikhonov regularization
        let mut options = LeastSquaresOptions::default();
        options.regularization_lambda = Some(0.1);
        options.method = LeastSquaresMethod::NormalEquations;
        
        let result = a.least_squares(&b, Some(&options)).unwrap();
        
        // Regularized solution should be stable
        let x_data = result.solution.to_cpu().unwrap();
        assert!(x_data.len() == 3);
        
        // Should have smaller norm than unregularized solution
        let x_norm_sq: f32 = x_data.iter().map(|x| x * x).sum();
        assert!(x_norm_sq.is_finite());
        
        assert!(result.condition_number.is_finite());
        assert!(result.residual_norm.is_finite());
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_least_squares_rank_deficient() {
        let device = Device::cpu();
        
        // Create rank-deficient matrix (rank 2, but 3 columns)
        let a = Tensor::from_data(
            vec![
                1.0, 2.0, 3.0,
                2.0, 4.0, 6.0,  // Row 2 = 2 * Row 1
                1.0, 1.0, 2.0,
            ],
            vec![3, 3],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(vec![6.0, 12.0, 4.0], vec![3], &device).unwrap();
        
        let result = a.least_squares(&b, None).unwrap();
        
        // Should detect rank deficiency
        assert!(result.rank < 3);
        assert!(result.condition_number > 1e6); // High condition number indicates rank deficiency
        
        // Solution should still satisfy Ax = b approximately
        let ax = a.matmul(&result.solution).unwrap();
        let ax_data = ax.to_cpu().unwrap();
        let b_data = b.to_cpu().unwrap();
        
        for (&expected, &actual) in b_data.iter().zip(ax_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-3);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_least_squares_method_selection() {
        let (a, b, _) = create_overdetermined_test_system();
        
        // Test automatic method selection
        let result_auto = a.least_squares(&b, None).unwrap();
        
        // Test explicit QR method
        let mut options_qr = LeastSquaresOptions::default();
        options_qr.method = LeastSquaresMethod::QR;
        let result_qr = a.least_squares(&b, Some(&options_qr)).unwrap();
        
        // Test explicit SVD method
        let mut options_svd = LeastSquaresOptions::default();
        options_svd.method = LeastSquaresMethod::SVD;
        let result_svd = a.least_squares(&b, Some(&options_svd)).unwrap();
        
        // All methods should give similar results for well-conditioned problems
        let auto_data = result_auto.solution.to_cpu().unwrap();
        let qr_data = result_qr.solution.to_cpu().unwrap();
        let svd_data = result_svd.solution.to_cpu().unwrap();
        
        for i in 0..3 {
            assert_abs_diff_eq!(auto_data[i], qr_data[i], epsilon = 1e-6);
            assert_abs_diff_eq!(auto_data[i], svd_data[i], epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_least_squares_goodness_of_fit() {
        let (a, b, _) = create_overdetermined_test_system();
        
        let result = a.least_squares(&b, None).unwrap();
        
        // Check goodness of fit metrics
        assert!(result.residual_norm >= 0.0);
        assert!(result.condition_number > 0.0);
        assert!(result.rank > 0);
        assert!(result.rank <= a.shape().dims()[1]); // rank <= number of columns
        
        // R-squared should be between 0 and 1
        if let Some(r_squared) = result.r_squared {
            assert!(r_squared >= 0.0 && r_squared <= 1.0);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_least_squares_multiple_rhs() {
        let device = Device::cpu();
        
        // Test with multiple right-hand sides
        let a = Tensor::from_data(
            vec![
                1.0, 2.0,
                3.0, 4.0,
                5.0, 6.0,
            ],
            vec![3, 2],
            &device,
        ).unwrap();
        
        let b = Tensor::from_data(
            vec![
                5.0, 11.0,  // First RHS: [5, 11]
                11.0, 25.0, // Second RHS: [11, 25]  
                17.0, 39.0, // Third RHS: [17, 39]
            ],
            vec![3, 2],
            &device,
        ).unwrap();
        
        let result = a.least_squares(&b, None).unwrap();
        
        // Solution should be 2x2 matrix (2 variables, 2 RHS)
        assert_eq!(result.solution.shape().dims(), &[2, 2]);
        
        // Each column should be a valid solution
        let x_data = result.solution.to_cpu().unwrap();
        assert!(x_data.len() == 4);
    }

    #[test]
    fn test_least_squares_invalid_dimensions() {
        let device = Device::cpu();
        
        // Test dimension mismatch
        let a = Tensor::from_data(vec![1.0, 2.0, 3.0, 4.0], vec![2, 2], &device).unwrap();
        let b = Tensor::from_data(vec![1.0, 2.0, 3.0], vec![3], &device).unwrap(); // Wrong size
        
        let result = a.least_squares(&b, None);
        assert!(result.is_err());
    }

    #[test]
    fn test_least_squares_empty_matrix() {
        let device = Device::cpu();
        
        // Test empty matrix
        let a = Tensor::zeros([0, 2], &device).unwrap();
        let b = Tensor::zeros([0], &device).unwrap();
        
        let result = a.least_squares(&b, None);
        assert!(result.is_err());
    }
}

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

    /// Helper function to create default LU options
    fn default_lu_options() -> LUOptions {
        LUOptions {
            pivoting: true,
            in_place: false,
            return_permutation_matrix: false,
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_decomposition_simple_2x2() {
        let device = Device::cpu();
        
        // Create a simple 2x2 invertible matrix
        // A = [[2, 1], [1, 1]]
        let a = Tensor::from_data(
            vec![2.0, 1.0, 1.0, 1.0],
            vec![2, 2],
            &device,
        ).unwrap();

        let options = default_lu_options();
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // Verify dimensions
        assert_eq!(lu_result.l.shape().dims(), &[2, 2]);
        assert_eq!(lu_result.u.shape().dims(), &[2, 2]);
        assert_eq!(lu_result.p.shape().dims(), &[2]);

        // Verify L is lower triangular with 1s on diagonal (Doolittle)
        let l_data = lu_result.l.to_cpu().unwrap();
        assert_abs_diff_eq!(l_data[0], 1.0, epsilon = 1e-8); // L[0,0] = 1
        assert_abs_diff_eq!(l_data[1], 0.0, epsilon = 1e-8); // L[0,1] = 0
        assert_abs_diff_eq!(l_data[3], 1.0, epsilon = 1e-8); // L[1,1] = 1

        // Verify U is upper triangular
        let u_data = lu_result.u.to_cpu().unwrap();
        assert_abs_diff_eq!(u_data[2], 0.0, epsilon = 1e-8); // U[1,0] = 0

        // Verify reconstruction: P*A = L*U
        let p_permuted_a = lu_result.apply_permutation(&a).unwrap();
        let lu_product = lu_result.l.matmul(&lu_result.u).unwrap();
        
        let pa_data = p_permuted_a.to_cpu().unwrap();
        let lu_data = lu_product.to_cpu().unwrap();
        
        for (&expected, &actual) in pa_data.iter().zip(lu_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_decomposition_3x3_no_pivoting() {
        let device = Device::cpu();
        
        // Create 3x3 matrix that doesn't need pivoting
        // A = [[2, 1, 1], [4, 3, 3], [8, 7, 9]]
        let a = Tensor::from_data(
            vec![2.0, 1.0, 1.0, 4.0, 3.0, 3.0, 8.0, 7.0, 9.0],
            vec![3, 3],
            &device,
        ).unwrap();

        let mut options = default_lu_options();
        options.pivoting = false; // Test without pivoting
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // Verify dimensions
        assert_eq!(lu_result.l.shape().dims(), &[3, 3]);
        assert_eq!(lu_result.u.shape().dims(), &[3, 3]);
        
        // Without pivoting, permutation should be identity
        let p_data = lu_result.p.to_cpu().unwrap();
        assert_eq!(p_data, vec![0, 1, 2]); // Identity permutation
        
        // Verify reconstruction: A = L*U (no permutation)
        let lu_product = lu_result.l.matmul(&lu_result.u).unwrap();
        let a_data = a.to_cpu().unwrap();
        let lu_data = lu_product.to_cpu().unwrap();
        
        for (&expected, &actual) in a_data.iter().zip(lu_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_decomposition_with_pivoting_3x3() {
        let device = Device::cpu();
        
        // Create matrix that requires pivoting for numerical stability
        // A = [[1, 2, 3], [4, 5, 6], [7, 8, 10]]
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0],
            vec![3, 3],
            &device,
        ).unwrap();

        let options = default_lu_options(); // With pivoting
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // Verify L has 1s on diagonal (Doolittle algorithm)
        let l_data = lu_result.l.to_cpu().unwrap();
        assert_abs_diff_eq!(l_data[0], 1.0, epsilon = 1e-8); // L[0,0]
        assert_abs_diff_eq!(l_data[4], 1.0, epsilon = 1e-8); // L[1,1]
        assert_abs_diff_eq!(l_data[8], 1.0, epsilon = 1e-8); // L[2,2]
        
        // Verify L is lower triangular
        assert_abs_diff_eq!(l_data[1], 0.0, epsilon = 1e-8); // L[0,1]
        assert_abs_diff_eq!(l_data[2], 0.0, epsilon = 1e-8); // L[0,2]
        assert_abs_diff_eq!(l_data[5], 0.0, epsilon = 1e-8); // L[1,2]
        
        // Verify U is upper triangular
        let u_data = lu_result.u.to_cpu().unwrap();
        assert_abs_diff_eq!(u_data[3], 0.0, epsilon = 1e-8); // U[1,0]
        assert_abs_diff_eq!(u_data[6], 0.0, epsilon = 1e-8); // U[2,0]
        assert_abs_diff_eq!(u_data[7], 0.0, epsilon = 1e-8); // U[2,1]
        
        // Verify reconstruction: P*A = L*U
        let p_permuted_a = lu_result.apply_permutation(&a).unwrap();
        let lu_product = lu_result.l.matmul(&lu_result.u).unwrap();
        
        let pa_data = p_permuted_a.to_cpu().unwrap();
        let lu_data = lu_product.to_cpu().unwrap();
        
        for (&expected, &actual) in pa_data.iter().zip(lu_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_decomposition_rectangular_4x3() {
        let device = Device::cpu();
        
        // Create 4x3 matrix (more rows than columns)
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0],
            vec![4, 3],
            &device,
        ).unwrap();

        let options = default_lu_options();
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // For rectangular matrix m×n where m > n:
        // L should be m×n (4×3) and U should be n×n (3×3)
        assert_eq!(lu_result.l.shape().dims(), &[4, 3]);
        assert_eq!(lu_result.u.shape().dims(), &[3, 3]);
        assert_eq!(lu_result.p.shape().dims(), &[4]);
        
        // Verify reconstruction: P*A = L*U
        let p_permuted_a = lu_result.apply_permutation(&a).unwrap();
        let lu_product = lu_result.l.matmul(&lu_result.u).unwrap();
        
        let pa_data = p_permuted_a.to_cpu().unwrap();
        let lu_data = lu_product.to_cpu().unwrap();
        
        for (&expected, &actual) in pa_data.iter().zip(lu_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_decomposition_rectangular_3x4() {
        let device = Device::cpu();
        
        // Create 3x4 matrix (more columns than rows)
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0],
            vec![3, 4],
            &device,
        ).unwrap();

        let options = default_lu_options();
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // For rectangular matrix m×n where m < n:
        // L should be m×m (3×3) and U should be m×n (3×4)
        assert_eq!(lu_result.l.shape().dims(), &[3, 3]);
        assert_eq!(lu_result.u.shape().dims(), &[3, 4]);
        assert_eq!(lu_result.p.shape().dims(), &[3]);
        
        // Verify reconstruction: P*A = L*U
        let p_permuted_a = lu_result.apply_permutation(&a).unwrap();
        let lu_product = lu_result.l.matmul(&lu_result.u).unwrap();
        
        let pa_data = p_permuted_a.to_cpu().unwrap();
        let lu_data = lu_product.to_cpu().unwrap();
        
        for (&expected, &actual) in pa_data.iter().zip(lu_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_decomposition_singular_matrix() {
        let device = Device::cpu();
        
        // Create singular matrix (rank deficient)
        let a = Tensor::from_data(
            vec![1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 1.0, 2.0, 3.0], // Row 2 = 2*Row 1, Row 3 = Row 1
            vec![3, 3],
            &device,
        ).unwrap();

        let options = default_lu_options();
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // Should still complete but indicate singularity
        assert!(!lu_result.is_invertible);
        assert!(lu_result.determinant.abs() < 1e-10);
        
        // Verify reconstruction still works
        let p_permuted_a = lu_result.apply_permutation(&a).unwrap();
        let lu_product = lu_result.l.matmul(&lu_result.u).unwrap();
        
        let pa_data = p_permuted_a.to_cpu().unwrap();
        let lu_data = lu_product.to_cpu().unwrap();
        
        for (&expected, &actual) in pa_data.iter().zip(lu_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_solve_linear_system() {
        let device = Device::cpu();
        
        // Create invertible 3x3 system
        let a = Tensor::from_data(
            vec![2.0, 1.0, 1.0, 4.0, 3.0, 3.0, 8.0, 7.0, 9.0],
            vec![3, 3],
            &device,
        ).unwrap();
        
        // Right-hand side: b = [4, 10, 24]
        let b = Tensor::from_data(vec![4.0, 10.0, 24.0], vec![3], &device).unwrap();
        
        // Expected solution: x = [1, 1, 1] (verify: A*x = b)
        let options = default_lu_options();
        let lu_result = a.lu_decomposition(&options).unwrap();
        let solution = lu_result.solve(&b).unwrap();
        
        assert_eq!(solution.shape().dims(), &[3]);
        let sol_data = solution.to_cpu().unwrap();
        
        // Check solution
        assert_abs_diff_eq!(sol_data[0], 1.0, epsilon = 1e-6);
        assert_abs_diff_eq!(sol_data[1], 1.0, epsilon = 1e-6);
        assert_abs_diff_eq!(sol_data[2], 1.0, epsilon = 1e-6);
        
        // Verify solution by substitution: A*x = b
        let ax = a.matmul(&solution).unwrap();
        let ax_data = ax.to_cpu().unwrap();
        let b_data = b.to_cpu().unwrap();
        
        for (&expected, &actual) in b_data.iter().zip(ax_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_solve_multiple_rhs() {
        let device = Device::cpu();
        
        // Create system with multiple right-hand sides
        let a = Tensor::from_data(
            vec![2.0, 1.0, 1.0, 1.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        // Multiple RHS: B = [[3, 4], [2, 3]]
        let b = Tensor::from_data(
            vec![3.0, 4.0, 2.0, 3.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let options = default_lu_options();
        let lu_result = a.lu_decomposition(&options).unwrap();
        let solution = lu_result.solve(&b).unwrap();
        
        // Should return 2x2 solution matrix
        assert_eq!(solution.shape().dims(), &[2, 2]);
        
        // Verify each column is correct solution
        let ax = a.matmul(&solution).unwrap();
        let ax_data = ax.to_cpu().unwrap();
        let b_data = b.to_cpu().unwrap();
        
        for (&expected, &actual) in b_data.iter().zip(ax_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_determinant_computation() {
        let device = Device::cpu();
        
        // Test determinant computation from LU decomposition
        let a = Tensor::from_data(
            vec![3.0, 2.0, 1.0, 4.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let options = default_lu_options();
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // Expected determinant: 3*4 - 2*1 = 10
        let expected_det = 10.0;
        assert_abs_diff_eq!(lu_result.determinant, expected_det, epsilon = 1e-6);
        
        // Compare with direct computation
        let direct_det = a.det().unwrap();
        assert_abs_diff_eq!(lu_result.determinant, direct_det, epsilon = 1e-6);
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_in_place_decomposition() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![2.0, 1.0, 1.0, 1.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let mut options = default_lu_options();
        options.in_place = true; // In-place decomposition
        
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // Should still work correctly
        assert_eq!(lu_result.l.shape().dims(), &[2, 2]);
        assert_eq!(lu_result.u.shape().dims(), &[2, 2]);
        
        // Verify reconstruction
        let p_permuted_a = lu_result.apply_permutation(&a).unwrap();
        let lu_product = lu_result.l.matmul(&lu_result.u).unwrap();
        
        let pa_data = p_permuted_a.to_cpu().unwrap();
        let lu_data = lu_product.to_cpu().unwrap();
        
        for (&expected, &actual) in pa_data.iter().zip(lu_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_permutation_matrix_output() {
        let device = Device::cpu();
        
        let a = Tensor::from_data(
            vec![1.0, 2.0, 4.0, 5.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let mut options = default_lu_options();
        options.return_permutation_matrix = true; // Return P as matrix
        
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // P should be returned as a matrix instead of vector
        assert!(lu_result.p_matrix.is_some());
        let p_matrix = lu_result.p_matrix.unwrap();
        assert_eq!(p_matrix.shape().dims(), &[2, 2]);
        
        // P should be orthogonal: P * P^T = I
        let p_transpose = p_matrix.transpose(0, 1).unwrap();
        let pp_t = p_matrix.matmul(&p_transpose).unwrap();
        let identity = Tensor::eye(2, &device).unwrap();
        
        let pp_t_data = pp_t.to_cpu().unwrap();
        let identity_data = identity.to_cpu().unwrap();
        
        for (&expected, &actual) in identity_data.iter().zip(pp_t_data.iter()) {
            assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
        }
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_numerical_stability() {
        let device = Device::cpu();
        
        // Test with ill-conditioned matrix
        let a = Tensor::from_data(
            vec![1e-10, 1.0, 1.0, 1.0],
            vec![2, 2],
            &device,
        ).unwrap();
        
        let options = default_lu_options(); // With pivoting for stability
        let lu_result = a.lu_decomposition(&options).unwrap();
        
        // Should handle without numerical issues
        let lu_product = lu_result.l.matmul(&lu_result.u).unwrap();
        let lu_data = lu_product.to_cpu().unwrap();
        
        // Results should be finite
        for &val in lu_data.iter() {
            assert!(val.is_finite(), "LU decomposition produced non-finite values");
        }
        
        // Condition number should be computed
        assert!(lu_result.condition_number.is_finite());
        assert!(lu_result.condition_number > 1.0);
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_batch_decomposition() {
        let device = Device::cpu();
        
        // Create batch of matrices: 3x2x2 (3 matrices of size 2x2)
        let batch_a = Tensor::from_data(
            vec![
                // First matrix
                1.0, 2.0, 3.0, 4.0,
                // Second matrix  
                2.0, 1.0, 1.0, 2.0,
                // Third matrix
                3.0, 1.0, 1.0, 3.0,
            ],
            vec![3, 2, 2],
            &device,
        ).unwrap();
        
        let options = default_lu_options();
        let batch_lu_result = batch_a.lu_decomposition_batch(&options).unwrap();
        
        // Should return batched results
        assert_eq!(batch_lu_result.l.shape().dims(), &[3, 2, 2]);
        assert_eq!(batch_lu_result.u.shape().dims(), &[3, 2, 2]);
        assert_eq!(batch_lu_result.p.shape().dims(), &[3, 2]);
        
        // Verify each matrix in batch
        for i in 0..3 {
            let a_i = batch_a.select(0, i).unwrap();
            let l_i = batch_lu_result.l.select(0, i).unwrap();
            let u_i = batch_lu_result.u.select(0, i).unwrap();
            
            // Verify reconstruction for each matrix
            let lu_product = l_i.matmul(&u_i).unwrap();
            let p_permuted_a = batch_lu_result.apply_permutation_batch(&a_i, i).unwrap();
            
            let pa_data = p_permuted_a.to_cpu().unwrap();
            let lu_data = lu_product.to_cpu().unwrap();
            
            for (&expected, &actual) in pa_data.iter().zip(lu_data.iter()) {
                assert_abs_diff_eq!(expected, actual, epsilon = 1e-6);
            }
        }
    }

    #[test]
    fn test_lu_error_conditions() {
        let device = Device::cpu();
        
        // Test with non-2D tensor
        let vector = Tensor::from_data(vec![1.0, 2.0, 3.0], vec![3], &device).unwrap();
        let options = default_lu_options();
        let result = vector.lu_decomposition(&options);
        assert!(result.is_err());
        
        // Test with empty matrix
        let empty = Tensor::zeros([0, 0], &device).unwrap();
        let result = empty.lu_decomposition(&options);
        assert!(result.is_err());
        
        // Test dimension mismatch in solve
        let a = Tensor::from_data(vec![1.0, 2.0, 3.0, 4.0], vec![2, 2], &device).unwrap();
        let lu_result = a.lu_decomposition(&options);
        // This will fail due to unimplemented, but the error checking structure is set up
        assert!(lu_result.is_err()); // Will fail because implementation doesn't exist yet
    }

    #[test]
    #[should_panic] // This test should fail until implementation is complete
    fn test_lu_vs_other_decompositions() {
        let device = Device::cpu();
        
        // Test consistency with other decomposition methods
        let a = Tensor::from_data(
            vec![4.0, 2.0, 1.0, 2.0, 3.0, 0.5, 1.0, 0.5, 2.0],
            vec![3, 3],
            &device,
        ).unwrap();
        
        // Solve system using LU
        let b = Tensor::from_data(vec![7.0, 5.5, 3.5], vec![3], &device).unwrap();
        
        let lu_options = default_lu_options();
        let lu_result = a.lu_decomposition(&lu_options).unwrap();
        let lu_solution = lu_result.solve(&b).unwrap();
        
        // Solve same system using QR
        let qr_result = a.qr().unwrap();
        let qt = qr_result.q.transpose(0, 1).unwrap();
        let qtb = qt.matmul(&b).unwrap();
        let triangular_options = TriangularSolveOptions {
            triangular_type: TriangularType::Upper,
            transpose: false,
            unit_diagonal: false,
            validate_triangular: false,
            strict_validation: false,
        };
        let qr_solution = qr_result.r.solve_triangular(&qtb, &triangular_options).unwrap();
        
        // Solutions should be similar
        let lu_data = lu_solution.to_cpu().unwrap();
        let qr_data = qr_solution.to_cpu().unwrap();
        
        for (&lu_val, &qr_val) in lu_data.iter().zip(qr_data.iter()) {
            assert_abs_diff_eq!(lu_val, qr_val, epsilon = 1e-5);
        }
    }
}