scirs2-linalg 0.4.2

//! Sparse eigenvalue decomposition for large sparse matrices
//!
//! This module provides efficient algorithms for computing eigenvalues and eigenvectors
//! of large sparse matrices. These algorithms are particularly useful when:
//! - Only a few eigenvalues/eigenvectors are needed
//! - The matrix is too large to fit in memory as a dense matrix
//! - The matrix has a high sparsity ratio
//!
//! ## Planned Algorithms
//!
//! - **Lanczos Algorithm**: For symmetric sparse matrices, finding extreme eigenvalues
//! - **Arnoldi Method**: For non-symmetric sparse matrices, finding eigenvalues near a target
//! - **Shift-and-Invert**: For finding interior eigenvalues efficiently
//! - **Jacobi-Davidson**: For generalized sparse eigenvalue problems
//!
//! ## Future Implementation
//!
//! This module currently provides placeholder implementations and will be fully
//! implemented in future versions to support:
//! - CSR (Compressed Sparse Row) matrix format
//! - Integration with external sparse linear algebra libraries
//! - Memory-efficient iterative solvers
//! - Parallel sparse matrix operations

use scirs2_core::ndarray::{s, Array1, Array2, ArrayView1, ArrayView2, ScalarOperand};
use scirs2_core::numeric::Complex;
use scirs2_core::numeric::{Float, NumAssign};
use scirs2_core::random::prelude::*;
use std::iter::Sum;

use crate::error::{LinalgError, LinalgResult};

// Type alias for sparse eigenvalue results
pub type SparseEigenResult<F> = LinalgResult<(Array1<Complex<F>>, Array2<Complex<F>>)>;

// Type alias for QR decomposition results
pub type QrResult<F> = LinalgResult<(Array2<Complex<F>>, Array2<Complex<F>>)>;

/// Sparse matrix trait for eigenvalue computations
///
/// This trait defines the interface that sparse matrix types should implement
/// to be compatible with sparse eigenvalue algorithms.
pub trait SparseMatrix<F> {
    /// Get the number of rows
    fn nrows(&self) -> usize;

    /// Get the number of columns  
    fn ncols(&self) -> usize;

    /// Matrix-vector multiplication: y = A * x
    fn matvec(&self, x: &ArrayView1<F>, y: &mut Array1<F>) -> LinalgResult<()>;

    /// Check if the matrix is symmetric
    fn is_symmetric(&self) -> bool;

    /// Get the sparsity ratio (number of non-zeros / total elements)
    fn sparsity(&self) -> f64;
}

/// Compute a few eigenvalues and eigenvectors of a large sparse matrix using Lanczos algorithm.
///
/// The Lanczos algorithm is an iterative method that is particularly effective for
/// symmetric sparse matrices when only a few eigenvalues are needed.
///
/// # Arguments
///
/// * `matrix` - Sparse matrix implementing the SparseMatrix trait
/// * `k` - Number of eigenvalues to compute
/// * `which` - Which eigenvalues to find ("largest", "smallest", "target")
/// * `target` - Target value for "target" mode (ignored for other modes)
/// * `max_iter` - Maximum number of iterations
/// * `tol` - Convergence tolerance
///
/// # Returns
///
/// * Tuple (eigenvalues, eigenvectors) with k eigenvalues and eigenvectors
///
/// # Examples
///
/// ```rust,ignore
/// use scirs2_linalg::eigen::sparse::{lanczos, SparseMatrix};
///
/// // This is a placeholder example - actual implementation pending
/// // let sparsematrix = create_sparsematrix();
/// // let (w, v) = lanczos(&sparsematrix, 5, "largest", 0.0, 100, 1e-6).expect("Operation failed");
/// ```
///
/// # Note
///
/// This function implements a parallel Lanczos algorithm for symmetric sparse matrices.
#[allow(dead_code)]
pub fn lanczos<F, M>(
    matrix: &M,
    k: usize,
    which: &str,
    target: F,
    max_iter: usize,
    tol: F,
) -> SparseEigenResult<F>
where
    F: Float + NumAssign + Sum + Send + Sync + 'static + Default,
    M: SparseMatrix<F> + Sync,
{
    let n = matrix.nrows();
    if n != matrix.ncols() {
        return Err(LinalgError::ShapeError(
            "Matrix must be square for eigenvalue decomposition".to_string(),
        ));
    }

    if k >= n {
        return Err(LinalgError::InvalidInputError(
            "Number of eigenvalues requested must be less than matrix size".to_string(),
        ));
    }

    // Store all Lanczos basis vectors for full reorthogonalization
    // (prevents loss of orthogonality in finite-precision arithmetic)
    let max_steps = max_iter.min(n);
    let mut v_basis: Vec<Array1<F>> = Vec::with_capacity(max_steps + 1);

    // Random initial vector
    let mut rng = scirs2_core::random::rng();
    let mut v_init = Array1::<F>::zeros(n);
    for i in 0..n {
        v_init[i] = F::from(rng.random::<f64>()).unwrap_or(F::one());
    }

    // Normalize initial vector
    let norm = v_init.iter().map(|x| (*x) * (*x)).sum::<F>().sqrt();
    if norm < F::from(1e-15).unwrap_or(F::epsilon()) {
        return Err(LinalgError::InvalidInputError(
            "Initial Lanczos vector has near-zero norm".to_string(),
        ));
    }
    v_init.mapv_inplace(|x| x / norm);
    v_basis.push(v_init);

    // Tridiagonal matrix elements
    let mut alpha = Vec::with_capacity(max_steps);
    let mut beta = Vec::with_capacity(max_steps);

    // Main Lanczos iteration with full reorthogonalization (Paige's implementation)
    // Storing all basis vectors prevents the numerical breakdown that occurs
    // with the three-term recurrence alone in finite-precision arithmetic.
    let mut v_next = Array1::<F>::zeros(n);
    for iter in 0..max_steps {
        let v_curr = &v_basis[iter];

        // Matrix-vector multiplication: w = A * v_curr
        matrix.matvec(&v_curr.view(), &mut v_next)?;

        // Three-term recurrence: subtract beta_{j-1} * v_{j-1}
        if iter > 0 {
            let beta_prev = beta[iter - 1];
            let v_prev = &v_basis[iter - 1];
            for j in 0..n {
                v_next[j] -= beta_prev * v_prev[j];
            }
        }

        // Compute alpha_j = <v_j, w>
        let alpha_curr = v_curr
            .iter()
            .zip(v_next.iter())
            .map(|(v, w)| (*v) * (*w))
            .sum::<F>();
        alpha.push(alpha_curr);

        // Subtract alpha_j * v_j
        for j in 0..n {
            v_next[j] -= alpha_curr * v_curr[j];
        }

        // Full reorthogonalization (modified Gram-Schmidt) against all previous vectors.
        // This is critical for numerical stability: without it, the three-term recurrence
        // accumulates floating-point errors that cause spurious (ghost) eigenvalues.
        // Two passes of MGS provide near-full double orthogonality.
        for _pass in 0..2 {
            for prev_v in v_basis.iter() {
                let proj = prev_v
                    .iter()
                    .zip(v_next.iter())
                    .map(|(p, w)| (*p) * (*w))
                    .sum::<F>();
                for j in 0..n {
                    v_next[j] -= proj * prev_v[j];
                }
            }
        }

        // Compute beta_{j+1} = ||w||
        let beta_curr = v_next.iter().map(|x| (*x) * (*x)).sum::<F>().sqrt();

        // Check for invariant subspace (exact breakdown)
        if beta_curr < tol {
            break;
        }

        beta.push(beta_curr);

        // Normalize for next iteration
        let v_new = v_next.mapv(|x| x / beta_curr);
        v_basis.push(v_new);
        v_next = Array1::<F>::zeros(n);

        // Check convergence of eigenvalues every few iterations
        if iter >= k && iter % 5 == 0 && check_lanczos_convergence(&alpha, &beta, k, tol) {
            break;
        }
    }

    // Solve tridiagonal eigenvalue problem
    let (eigenvals, eigenvecs) = solve_tridiagonal_eigenproblem(&alpha, &beta, which, target, k)?;

    // Convert to complex format
    let complex_eigenvals = eigenvals.mapv(|x| Complex::new(x, F::zero()));
    let complex_eigenvecs = eigenvecs.mapv(|x| Complex::new(x, F::zero()));

    Ok((complex_eigenvals, complex_eigenvecs))
}

// Helper function to check Lanczos convergence
#[allow(dead_code)]
fn check_lanczos_convergence<F: Float>(_alpha: &[F], beta: &[F], k: usize, tol: F) -> bool {
    // Simple convergence check based on beta values
    if beta.len() < k {
        return false;
    }

    let recent_betas = &beta[beta.len().saturating_sub(k)..];
    recent_betas
        .iter()
        .all(|&b| b < tol * F::from(10.0).expect("Operation failed"))
}

// Helper function to solve tridiagonal eigenvalue problem
#[allow(dead_code)]
fn solve_tridiagonal_eigenproblem<F: Float + NumAssign + Sum + Send + Sync + 'static>(
    alpha: &[F],
    beta: &[F],
    which: &str,
    target: F,
    k: usize,
) -> LinalgResult<(Array1<F>, Array2<F>)> {
    let n = alpha.len();
    if n == 0 {
        return Err(LinalgError::InvalidInputError(
            "Empty tridiagonal matrix".to_string(),
        ));
    }

    // Create tridiagonal matrix
    let mut trimatrix = Array2::<F>::zeros((n, n));

    // Fill diagonal
    for i in 0..n {
        trimatrix[[i, i]] = alpha[i];
    }

    // Fill off-diagonals
    for i in 0..n.saturating_sub(1) {
        if i < beta.len() {
            trimatrix[[i, i + 1]] = beta[i];
            trimatrix[[i + 1, i]] = beta[i];
        }
    }

    // Use QR algorithm for small tridiagonal matrices
    let (eigenvals, eigenvecs) = qr_algorithm_tridiagonal(&trimatrix)?;

    // Select requested eigenvalues based on 'which' parameter
    let selected_indices = select_eigenvalues(&eigenvals, which, target, k);

    let mut result_eigenvals = Array1::<F>::zeros(k);
    let mut result_eigenvecs = Array2::<F>::zeros((n, k));

    for (i, &idx) in selected_indices.iter().enumerate() {
        result_eigenvals[i] = eigenvals[idx];
        for j in 0..n {
            result_eigenvecs[[j, i]] = eigenvecs[[j, idx]];
        }
    }

    Ok((result_eigenvals, result_eigenvecs))
}

// Helper function for QR algorithm on tridiagonal matrices with Wilkinson shift.
//
// Bug fix: previously `qr_decomposition_tridiagonal` returned Q^T (the product of
// the Givens rotation matrices G_k), not Q = G_k^T. This caused `r.dot(&q)` to
// compute R*Q^T instead of R*Q, yielding incorrect eigenvalues.
// Fixed by calling `qr_decomposition_tridiagonal` which now returns the correct Q.
#[allow(dead_code)]
fn qr_algorithm_tridiagonal<F: Float + NumAssign + Sum + 'static>(
    matrix: &Array2<F>,
) -> LinalgResult<(Array1<F>, Array2<F>)> {
    let n = matrix.nrows();
    let mut a = matrix.clone();
    let mut q_total = Array2::<F>::eye(n);

    let max_iterations = 1000;
    let tolerance = F::from(1e-12).unwrap_or(F::epsilon());

    for _iter in 0..max_iterations {
        // Check for convergence: all subdiagonals near zero
        let mut converged = true;
        for i in 0..n - 1 {
            if a[[i + 1, i]].abs() > tolerance {
                converged = false;
                break;
            }
        }

        if converged {
            break;
        }

        // Wilkinson shift: use eigenvalue of the bottom-right 2x2 submatrix
        // closest to a[[n-1, n-1]] to accelerate convergence.
        let shift = if n >= 2 {
            let d = (a[[n - 2, n - 2]] - a[[n - 1, n - 1]]) / F::from(2.0).unwrap_or(F::one());
            let b_sq = a[[n - 1, n - 2]] * a[[n - 1, n - 2]];
            let sign_d = if d >= F::zero() { F::one() } else { -F::one() };
            a[[n - 1, n - 1]] - sign_d * b_sq / (d.abs() + (d * d + b_sq).sqrt())
        } else {
            a[[0, 0]]
        };

        // Shift: A - mu*I
        for i in 0..n {
            a[[i, i]] -= shift;
        }

        // QR decomposition step — returns the correct orthogonal Q
        let (q, r) = qr_decomposition_tridiagonal(&a)?;

        // A_new = R * Q + mu*I  (equivalent to Q^T * A * Q unshifted, then restore shift)
        a = r.dot(&q);

        // Restore shift
        for i in 0..n {
            a[[i, i]] += shift;
        }

        // Accumulate eigenvectors: V = V * Q
        q_total = q_total.dot(&q);
    }

    // Extract eigenvalues from diagonal
    let eigenvals = (0..n).map(|i| a[[i, i]]).collect::<Array1<F>>();

    Ok((eigenvals, q_total))
}

// QR decomposition for tridiagonal (and general) matrices via Givens rotations.
//
// Returns (Q, R) such that matrix = Q * R, where Q is orthogonal and R is upper triangular.
//
// Bug fix: the previous implementation accumulated `q = G_1 * G_2 * ...` (product of
// left-multiplication Givens matrices), which equals Q^T, not Q. This was because
// `apply_givens_rotation_transpose` applied column operations equivalent to right-multiplying
// by G_k. Since G_k * A = R implies A = G_k^T * R = Q * R with Q = G_k^T, the correct Q
// is the transpose of the accumulated product. We now return q.t() to get the correct Q.
#[allow(dead_code)]
fn qr_decomposition_tridiagonal<F: Float + NumAssign + Sum>(
    matrix: &Array2<F>,
) -> LinalgResult<(Array2<F>, Array2<F>)> {
    let n = matrix.nrows();
    // g_product accumulates G_1 * G_2 * ... (product of row-Givens matrices)
    // which equals Q^T. We return its transpose to get Q.
    let mut g_product = Array2::<F>::eye(n);
    let mut r = matrix.clone();

    let eps = F::from(1e-15).unwrap_or(F::epsilon());

    // Use Givens rotations for tridiagonal matrices
    for i in 0..n - 1 {
        let a = r[[i, i]];
        let b = r[[i + 1, i]];

        if b.abs() > eps {
            let (c, s) = givens_rotation(a, b);

            // Left-multiply R by G_k (zeroes out r[i+1, i])
            apply_givens_rotation(&mut r, i, i + 1, c, s);

            // Accumulate G_k into g_product (right-multiply by G_k)
            apply_givens_rotation_transpose(&mut g_product, i, i + 1, c, s);
        }
    }

    // g_product = G_1 * G_2 * ... = Q^T, so Q = g_product^T
    let q = g_product.t().to_owned();

    Ok((q, r))
}

// Helper function for Givens rotation
#[allow(dead_code)]
fn givens_rotation<F: Float>(a: F, b: F) -> (F, F) {
    if b.abs() < F::from(1e-15).expect("Operation failed") {
        (F::one(), F::zero())
    } else {
        let r = (a * a + b * b).sqrt();
        (a / r, -b / r)
    }
}

// Apply Givens rotation to matrix
#[allow(dead_code)]
fn apply_givens_rotation<F: Float + NumAssign>(
    matrix: &mut Array2<F>,
    i: usize,
    j: usize,
    c: F,
    s: F,
) {
    let n = matrix.ncols();
    for k in 0..n {
        let temp1 = matrix[[i, k]];
        let temp2 = matrix[[j, k]];
        matrix[[i, k]] = c * temp1 - s * temp2;
        matrix[[j, k]] = s * temp1 + c * temp2;
    }
}

// Apply Givens rotation transpose to matrix
#[allow(dead_code)]
fn apply_givens_rotation_transpose<F: Float + NumAssign>(
    matrix: &mut Array2<F>,
    i: usize,
    j: usize,
    c: F,
    s: F,
) {
    let n = matrix.nrows();
    for k in 0..n {
        let temp1 = matrix[[k, i]];
        let temp2 = matrix[[k, j]];
        matrix[[k, i]] = c * temp1 + s * temp2;
        matrix[[k, j]] = -s * temp1 + c * temp2;
    }
}

// Helper function to select eigenvalues based on criteria
#[allow(dead_code)]
fn select_eigenvalues<F: Float>(
    eigenvals: &Array1<F>,
    which: &str,
    target: F,
    k: usize,
) -> Vec<usize> {
    let mut indices_and_values: Vec<(usize, F)> = eigenvals
        .iter()
        .enumerate()
        .map(|(i, &val)| (i, val))
        .collect();

    match which {
        "largest" | "LM" => {
            indices_and_values.sort_by(|a, b| b.1.partial_cmp(&a.1).expect("Operation failed"));
        }
        "smallest" | "SM" => {
            indices_and_values.sort_by(|a, b| a.1.partial_cmp(&b.1).expect("Operation failed"));
        }
        "target" | "nearest" => {
            indices_and_values.sort_by(|a, b| {
                let dist_a = (a.1 - target).abs();
                let dist_b = (b.1 - target).abs();
                dist_a.partial_cmp(&dist_b).expect("Operation failed")
            });
        }
        _ => {
            // Default to largest
            indices_and_values.sort_by(|a, b| b.1.partial_cmp(&a.1).expect("Operation failed"));
        }
    }

    indices_and_values
        .into_iter()
        .take(k)
        .map(|(idx, _)| idx)
        .collect()
}

/// Compute eigenvalues near a target value using the Arnoldi method.
///
/// The Arnoldi method is a generalization of the Lanczos algorithm that works
/// for non-symmetric matrices. It's particularly effective when combined with
/// shift-and-invert to find eigenvalues near a specific target value.
///
/// # Arguments
///
/// * `matrix` - Sparse matrix implementing the SparseMatrix trait
/// * `k` - Number of eigenvalues to compute
/// * `target` - Target eigenvalue around which to search
/// * `max_iter` - Maximum number of iterations
/// * `tol` - Convergence tolerance
///
/// # Returns
///
/// * Tuple (eigenvalues, eigenvectors) with k eigenvalues closest to target
///
/// # Examples
///
/// ```rust,ignore
/// use scirs2_linalg::eigen::sparse::{arnoldi, SparseMatrix};
///
/// // This is a placeholder example - actual implementation pending
/// // let sparsematrix = create_sparsematrix();
/// // let (w, v) = arnoldi(&sparsematrix, 3, 1.5, 100, 1e-6).expect("Operation failed");
/// ```
///
/// # Note
///
/// This function implements a parallel Arnoldi method for non-symmetric sparse matrices.
#[allow(dead_code)]
pub fn arnoldi<F, M>(
    matrix: &M,
    k: usize,
    target: Complex<F>,
    max_iter: usize,
    tol: F,
) -> SparseEigenResult<F>
where
    F: Float + NumAssign + Sum + Send + Sync + 'static + Default,
    M: SparseMatrix<F> + Sync,
{
    let n = matrix.nrows();
    if n != matrix.ncols() {
        return Err(LinalgError::ShapeError(
            "Matrix must be square for eigenvalue decomposition".to_string(),
        ));
    }

    if k >= n {
        return Err(LinalgError::InvalidInputError(
            "Number of eigenvalues requested must be less than matrix size".to_string(),
        ));
    }

    let m = (max_iter + 1).min(n);

    // Arnoldi vectors (Krylov basis)
    let mut v_vectors = vec![Array1::<F>::zeros(n); m + 1];

    // Hessenberg matrix
    let mut hmatrix = Array2::<F>::zeros((m + 1, m));

    // Random initial vector
    let mut rng = scirs2_core::random::rng();
    for v_elem in &mut v_vectors[0] {
        *v_elem = F::from(rng.random::<f64>()).expect("Operation failed");
    }

    // Normalize initial vector
    let norm = v_vectors[0].iter().map(|x| (*x) * (*x)).sum::<F>().sqrt();
    for v_elem in &mut v_vectors[0] {
        *v_elem /= norm;
    }

    // Main Arnoldi iteration
    let mut actual_m = 0;
    for j in 0..m {
        actual_m = j + 1;

        // Matrix-vector multiplication: w = A * v_j
        let mut w = Array1::<F>::zeros(n);
        matrix.matvec(&v_vectors[j].view(), &mut w)?;

        // Modified Gram-Schmidt orthogonalization
        for i in 0..=j {
            // h[i][j] = <w, v_i>
            let h_ij = w
                .iter()
                .zip(v_vectors[i].iter())
                .map(|(w_val, v_val)| (*w_val) * (*v_val))
                .sum::<F>();
            hmatrix[[i, j]] = h_ij;

            // w = w - h[i][j] * v_i
            for l in 0..n {
                w[l] -= h_ij * v_vectors[i][l];
            }
        }

        // h[j+1][j] = ||w||
        let h_j1_j = w.iter().map(|x| (*x) * (*x)).sum::<F>().sqrt();

        // Check for breakdown or convergence
        if h_j1_j < tol {
            break;
        }

        if j + 1 < m {
            hmatrix[[j + 1, j]] = h_j1_j;

            // v_{j+1} = w / h[j+1][j]
            for l in 0..n {
                v_vectors[j + 1][l] = w[l] / h_j1_j;
            }
        }

        // Check convergence of Ritz values every few iterations
        if j >= k && j % 5 == 0 && check_arnoldi_convergence(&hmatrix, j + 1, k, tol) {
            break;
        }
    }

    // Extract the m x m upper Hessenberg matrix
    let h_reduced = hmatrix.slice(s![..actual_m, ..actual_m]).to_owned();

    // Solve eigenvalue problem for Hessenberg matrix
    let (ritz_values, ritz_vectors) = solve_hessenberg_eigenproblem(&h_reduced)?;

    // Convert Ritz values to eigenvalue estimates
    let eigenvals = if target.im == F::zero() {
        // Real target - select closest real eigenvalues
        select_closest_real_eigenvalues(&ritz_values, target.re, k)
    } else {
        // Complex target - select closest eigenvalues
        select_closest_complex_eigenvalues(&ritz_values, target, k)
    };

    // Compute eigenvectors by combining Ritz vectors with Arnoldi basis
    let mut eigenvecs = Array2::<Complex<F>>::zeros((n, k));
    let v_basis = v_vectors[..actual_m]
        .iter()
        .map(|v| v.mapv(|x| Complex::new(x, F::zero())))
        .collect::<Vec<_>>();

    for (i, &ritz_idx) in eigenvals.iter().enumerate() {
        for j in 0..n {
            let mut eigenvec_j = Complex::new(F::zero(), F::zero());
            for l in 0..actual_m {
                eigenvec_j += ritz_vectors[[l, ritz_idx]] * v_basis[l][j];
            }
            eigenvecs[[j, i]] = eigenvec_j;
        }
    }

    let final_eigenvals = eigenvals
        .iter()
        .map(|&idx| ritz_values[idx])
        .collect::<Array1<_>>();

    Ok((final_eigenvals, eigenvecs))
}

// Helper function to check Arnoldi convergence
#[allow(dead_code)]
fn check_arnoldi_convergence<F: Float>(hmatrix: &Array2<F>, m: usize, k: usize, tol: F) -> bool {
    // Simple convergence check based on subdiagonal elements
    if m < k + 1 {
        return false;
    }

    // Check if the last k subdiagonal elements are small
    (0..k).all(|i| {
        let row = m - 1 - i;
        let col = m - 2 - i;
        if row < hmatrix.nrows() && col < hmatrix.ncols() {
            hmatrix[[row, col]].abs() < tol * F::from(10.0).expect("Operation failed")
        } else {
            true
        }
    })
}

// Helper function to solve Hessenberg eigenvalue problem
#[allow(dead_code)]
fn solve_hessenberg_eigenproblem<F: Float + NumAssign + Sum + 'static>(
    hmatrix: &Array2<F>,
) -> SparseEigenResult<F> {
    let n = hmatrix.nrows();

    // For simplicity, convert to general eigenvalue problem
    // In practice, specialized Hessenberg QR algorithm would be better
    let mut matrix_complex = Array2::<Complex<F>>::zeros((n, n));
    for i in 0..n {
        for j in 0..n {
            matrix_complex[[i, j]] = Complex::new(hmatrix[[i, j]], F::zero());
        }
    }

    // Use QR algorithm for complex matrices
    qr_algorithm_complex(&matrix_complex)
}

// Simplified QR algorithm for complex matrices
#[allow(dead_code)]
fn qr_algorithm_complex<F: Float + NumAssign + Sum + 'static>(
    matrix: &Array2<Complex<F>>,
) -> SparseEigenResult<F> {
    let n = matrix.nrows();
    let mut a = matrix.clone();
    let mut q_total = Array2::<Complex<F>>::eye(n);

    let max_iterations = 1000;
    let tolerance = F::from(1e-12).expect("Operation failed");

    for _iter in 0..max_iterations {
        // Check for convergence (simplified)
        let mut converged = true;
        for i in 0..n - 1 {
            if a[[i + 1, i]].norm() > tolerance {
                converged = false;
                break;
            }
        }

        if converged {
            break;
        }

        // Simplified QR step (this should use specialized Hessenberg QR)
        let (q, r) = householder_qr_complex(&a)?;
        a = r.dot(&q);
        q_total = q_total.dot(&q);
    }

    // Extract eigenvalues from diagonal
    let eigenvals = (0..n).map(|i| a[[i, i]]).collect::<Array1<_>>();

    Ok((eigenvals, q_total))
}

// Simplified Householder QR for complex matrices
#[allow(dead_code)]
fn householder_qr_complex<F: Float + NumAssign + Sum>(matrix: &Array2<Complex<F>>) -> QrResult<F> {
    let (m, n) = matrix.dim();
    let mut q = Array2::<Complex<F>>::eye(m);
    let mut r = matrix.clone();

    let min_dim = m.min(n);

    for k in 0..min_dim {
        // Extract column for Householder reflection
        let x = r.slice(s![k.., k]).to_owned();
        let (house_vec, tau) = householder_vector_complex(&x);

        // Apply Householder reflection to R
        apply_householder_left_complex(&mut r, &house_vec, tau, k);

        // Apply to Q (accumulate transformations)
        apply_householder_right_complex(&mut q, &house_vec, tau.conj(), k);
    }

    Ok((q, r))
}

// Helper function for complex Householder vector
#[allow(dead_code)]
fn householder_vector_complex<F: Float + NumAssign + Sum>(
    x: &Array1<Complex<F>>,
) -> (Array1<Complex<F>>, Complex<F>) {
    let n = x.len();
    if n == 0 {
        return (Array1::zeros(0), Complex::new(F::zero(), F::zero()));
    }

    let norm_x = x.iter().map(|z| z.norm_sqr()).sum::<F>().sqrt();

    if norm_x == F::zero() {
        return (Array1::zeros(n), Complex::new(F::zero(), F::zero()));
    }

    let mut v = x.clone();
    let sign = if x[0].re >= F::zero() {
        F::one()
    } else {
        -F::one()
    };
    v[0] += Complex::new(sign * norm_x, F::zero());

    let norm_v = v.iter().map(|z| z.norm_sqr()).sum::<F>().sqrt();
    if norm_v > F::zero() {
        v.mapv_inplace(|z| z / norm_v);
    }

    let tau = Complex::new(F::from(2.0).expect("Operation failed"), F::zero());

    (v, tau)
}

// Apply Householder reflection from left
#[allow(dead_code)]
fn apply_householder_left_complex<F: Float + NumAssign>(
    matrix: &mut Array2<Complex<F>>,
    house_vec: &Array1<Complex<F>>,
    tau: Complex<F>,
    k: usize,
) {
    let (m, n) = matrix.dim();
    let house_len = house_vec.len();

    for j in k..n {
        let mut sum = Complex::new(F::zero(), F::zero());
        for i in 0..house_len {
            if k + i < m {
                sum += house_vec[i].conj() * matrix[[k + i, j]];
            }
        }

        for i in 0..house_len {
            if k + i < m {
                matrix[[k + i, j]] -= tau * house_vec[i] * sum;
            }
        }
    }
}

// Apply Householder reflection from right
#[allow(dead_code)]
fn apply_householder_right_complex<F: Float + NumAssign>(
    matrix: &mut Array2<Complex<F>>,
    house_vec: &Array1<Complex<F>>,
    tau: Complex<F>,
    k: usize,
) {
    let (m, _n) = matrix.dim();
    let house_len = house_vec.len();

    for i in 0..m {
        let mut sum = Complex::new(F::zero(), F::zero());
        for j in 0..house_len {
            if k + j < matrix.ncols() {
                sum += matrix[[i, k + j]] * house_vec[j];
            }
        }

        for j in 0..house_len {
            if k + j < matrix.ncols() {
                matrix[[i, k + j]] -= sum * tau.conj() * house_vec[j].conj();
            }
        }
    }
}

// Helper functions for eigenvalue selection
#[allow(dead_code)]
fn select_closest_real_eigenvalues<F: Float>(
    eigenvals: &Array1<Complex<F>>,
    target: F,
    k: usize,
) -> Vec<usize> {
    let mut real_eigenvals: Vec<(usize, F)> = eigenvals
        .iter()
        .enumerate()
        .filter(|(_, z)| z.im.abs() < F::from(1e-10).expect("Operation failed"))
        .map(|(i, z)| (i, z.re))
        .collect();

    real_eigenvals.sort_by(|a, b| {
        let dist_a = (a.1 - target).abs();
        let dist_b = (b.1 - target).abs();
        dist_a.partial_cmp(&dist_b).expect("Operation failed")
    });

    real_eigenvals
        .into_iter()
        .take(k)
        .map(|(idx, _)| idx)
        .collect()
}

#[allow(dead_code)]
fn select_closest_complex_eigenvalues<F: Float>(
    eigenvals: &Array1<Complex<F>>,
    target: Complex<F>,
    k: usize,
) -> Vec<usize> {
    let mut eigenvals_with_dist: Vec<(usize, F)> = eigenvals
        .iter()
        .enumerate()
        .map(|(i, z)| {
            let diff = *z - target;
            (i, diff.norm())
        })
        .collect();

    eigenvals_with_dist.sort_by(|a, b| a.1.partial_cmp(&b.1).expect("Operation failed"));

    eigenvals_with_dist
        .into_iter()
        .take(k)
        .map(|(idx, _)| idx)
        .collect()
}

/// Solve sparse generalized eigenvalue problem Ax = λBx using iterative methods.
///
/// This function solves the generalized eigenvalue problem for sparse matrices
/// using specialized algorithms that avoid forming dense factorizations.
///
/// # Arguments
///
/// * `a` - Sparse matrix A
/// * `b` - Sparse matrix B (should be positive definite for symmetric case)
/// * `k` - Number of eigenvalues to compute
/// * `which` - Which eigenvalues to find ("largest", "smallest", "target")
/// * `target` - Target value for "target" mode
/// * `max_iter` - Maximum number of iterations
/// * `tol` - Convergence tolerance
///
/// # Returns
///
/// * Tuple (eigenvalues, eigenvectors) with k generalized eigenvalues and eigenvectors
///
/// # Examples
///
/// ```rust,ignore
/// use scirs2_linalg::eigen::sparse::{eigs_gen, SparseMatrix};
///
/// // This is a placeholder example - actual implementation pending
/// // let (w, v) = eigs_gen(&sparse_a, &sparse_b, 4, "smallest", 0.0, 100, 1e-6).expect("Operation failed");
/// ```
///
/// # Note
///
/// This function is currently a placeholder and will be implemented in a future version.
#[allow(dead_code)]
pub fn eigs_gen<F, M1, M2>(
    _a: &M1,
    _b: &M2,
    _k: usize,
    _which: &str,
    _target: F,
    _max_iter: usize,
    _tol: F,
) -> SparseEigenResult<F>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
    M1: SparseMatrix<F>,
    M2: SparseMatrix<F>,
{
    Err(LinalgError::NotImplementedError(
        "Sparse generalized eigenvalue solver not yet implemented".to_string(),
    ))
}

/// Compute singular values and vectors of a sparse matrix using iterative methods.
///
/// This function computes the largest or smallest singular values of a sparse matrix
/// without forming the normal equations, which can be numerically unstable for
/// ill-conditioned matrices.
///
/// # Arguments
///
/// * `matrix` - Sparse matrix
/// * `k` - Number of singular values to compute
/// * `which` - Which singular values to find ("largest" or "smallest")
/// * `max_iter` - Maximum number of iterations
/// * `tol` - Convergence tolerance
///
/// # Returns
///
/// * Tuple (singular_values, left_vectors, right_vectors)
///
/// # Examples
///
/// ```rust,ignore
/// use scirs2_linalg::eigen::sparse::{svds, SparseMatrix};
///
/// // This is a placeholder example - actual implementation pending
/// // let (s, u, vt) = svds(&sparsematrix, 6, "largest", 100, 1e-6).expect("Operation failed");
/// ```
///
/// # Note
///
/// This function is currently a placeholder and will be implemented in a future version.
#[allow(dead_code)]
pub fn svds<F, M>(
    matrix: &M,
    _k: usize,
    _which: &str,
    _max_iter: usize,
    _tol: F,
) -> LinalgResult<(Array1<F>, Array2<F>, Array2<F>)>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
    M: SparseMatrix<F>,
{
    Err(LinalgError::NotImplementedError(
        "Sparse SVD solver not yet implemented".to_string(),
    ))
}

/// Convert a dense matrix to sparse format for eigenvalue computations.
///
/// This is a utility function that can detect sparsity in dense matrices and
/// convert them to an appropriate sparse format for more efficient eigenvalue
/// computations when the matrix is sufficiently sparse.
///
/// # Arguments
///
/// * `densematrix` - Dense matrix to convert
/// * `threshold` - Sparsity threshold (elements with absolute value below this are considered zero)
///
/// # Returns
///
/// * A sparse matrix representation suitable for sparse eigenvalue algorithms
///
/// # Examples
///
/// ```rust,ignore
/// use scirs2_core::ndarray::Array2;
/// use scirs2_linalg::eigen::sparse::dense_to_sparse;
///
/// // This is a placeholder example - actual implementation pending
/// // let dense = Array2::eye(1000);
/// // let sparse = dense_to_sparse(&dense.view(), 1e-12).expect("Operation failed");
/// ```
///
/// # Note
///
/// This function is currently a placeholder and will be implemented in a future version.
/// CSR (Compressed Sparse Row) matrix for eigenvalue computations.
///
/// Stores non-zero entries in three arrays:
/// - `data`: the non-zero values, stored row by row
/// - `indices`: the column index for each non-zero
/// - `indptr`: row pointers -- `indptr[i]..indptr[i+1]` gives the range
///   into `data`/`indices` for row `i`
pub struct CsrMatrix<F> {
    nrows: usize,
    ncols: usize,
    data: Vec<F>,
    indices: Vec<usize>,
    indptr: Vec<usize>,
}

impl<F> CsrMatrix<F>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    /// Create a new CSR matrix.
    ///
    /// # Arguments
    /// * `nrows` - Number of rows
    /// * `ncols` - Number of columns
    /// * `data` - Non-zero values
    /// * `indices` - Column indices for each non-zero
    /// * `indptr` - Row pointers (length should be nrows + 1 for well-formed CSR)
    pub fn new(
        nrows: usize,
        ncols: usize,
        data: Vec<F>,
        indices: Vec<usize>,
        indptr: Vec<usize>,
    ) -> Self {
        Self {
            nrows,
            ncols,
            data,
            indices,
            indptr,
        }
    }

    /// Number of stored non-zero entries.
    pub fn nnz(&self) -> usize {
        self.data.len()
    }

    /// Create a CSR matrix from a dense array, dropping entries below `threshold`.
    pub fn from_dense(dense: &ArrayView2<F>, threshold: F) -> Self {
        let (m, n) = dense.dim();
        let mut data = Vec::new();
        let mut indices = Vec::new();
        let mut indptr = Vec::with_capacity(m + 1);

        indptr.push(0);
        for i in 0..m {
            for j in 0..n {
                let val = dense[[i, j]];
                if val.abs() >= threshold {
                    data.push(val);
                    indices.push(j);
                }
            }
            indptr.push(data.len());
        }

        Self {
            nrows: m,
            ncols: n,
            data,
            indices,
            indptr,
        }
    }
}

/// Convert a dense matrix to CSR sparse format for eigenvalue computations.
///
/// Scans the dense matrix and extracts entries with absolute value >= `threshold`.
///
/// # Arguments
/// * `densematrix` - Dense matrix to convert
/// * `threshold` - Sparsity threshold
///
/// # Returns
/// A `CsrMatrix` wrapped in `Box<dyn SparseMatrix<F>>`.
#[allow(dead_code)]
pub fn dense_to_sparse<F>(
    densematrix: &ArrayView2<F>,
    threshold: F,
) -> LinalgResult<Box<dyn SparseMatrix<F>>>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    Ok(Box::new(CsrMatrix::from_dense(densematrix, threshold)))
}

impl<F> SparseMatrix<F> for CsrMatrix<F>
where
    F: Float + NumAssign + Sum + Send + Sync + ScalarOperand + 'static,
{
    fn nrows(&self) -> usize {
        self.nrows
    }

    fn ncols(&self) -> usize {
        self.ncols
    }

    fn matvec(&self, x: &ArrayView1<F>, y: &mut Array1<F>) -> LinalgResult<()> {
        if x.len() != self.ncols {
            return Err(LinalgError::DimensionError(format!(
                "CsrMatrix::matvec: x has {} elements but matrix has {} columns",
                x.len(),
                self.ncols
            )));
        }
        if y.len() != self.nrows {
            return Err(LinalgError::DimensionError(format!(
                "CsrMatrix::matvec: y has {} elements but matrix has {} rows",
                y.len(),
                self.nrows
            )));
        }

        // CSR matvec: y[i] = sum_{k in row_range(i)} data[k] * x[indices[k]]
        for i in 0..self.nrows {
            let row_start = if i < self.indptr.len() {
                self.indptr[i]
            } else {
                self.data.len()
            };
            let row_end = if i + 1 < self.indptr.len() {
                self.indptr[i + 1]
            } else {
                self.data.len()
            };
            let mut acc = F::zero();
            for k in row_start..row_end {
                if k < self.data.len() && k < self.indices.len() {
                    let col = self.indices[k];
                    if col < x.len() {
                        acc += self.data[k] * x[col];
                    }
                }
            }
            y[i] = acc;
        }

        Ok(())
    }

    fn is_symmetric(&self) -> bool {
        if self.nrows != self.ncols {
            return false;
        }
        let n = self.nrows;
        // Check: for every (i, j, val), verify (j, i) has the same value
        for i in 0..n {
            let row_start = if i < self.indptr.len() {
                self.indptr[i]
            } else {
                return false;
            };
            let row_end = if i + 1 < self.indptr.len() {
                self.indptr[i + 1]
            } else {
                self.data.len()
            };
            for k in row_start..row_end {
                if k >= self.data.len() || k >= self.indices.len() {
                    continue;
                }
                let j = self.indices[k];
                if j >= n {
                    return false;
                }
                let val = self.data[k];
                // Search for (j, i) entry
                let ji_start = if j < self.indptr.len() {
                    self.indptr[j]
                } else {
                    return false;
                };
                let ji_end = if j + 1 < self.indptr.len() {
                    self.indptr[j + 1]
                } else {
                    self.data.len()
                };
                let mut found = false;
                for kk in ji_start..ji_end {
                    if kk < self.indices.len()
                        && self.indices[kk] == i
                        && kk < self.data.len()
                        && (self.data[kk] - val).abs() < F::from(1e-14).unwrap_or(F::epsilon())
                    {
                        found = true;
                        break;
                    }
                }
                if !found {
                    return false;
                }
            }
        }
        true
    }

    fn sparsity(&self) -> f64 {
        let total = (self.nrows as f64) * (self.ncols as f64);
        if total == 0.0 {
            return 0.0;
        }
        self.data.len() as f64 / total
    }
}

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

    /// Build a 5x5 symmetric tridiagonal SPD matrix in CSR format (1D Laplacian).
    fn tridiag_csr_5x5() -> CsrMatrix<f64> {
        let n = 5;
        let mut data = Vec::new();
        let mut indices = Vec::new();
        let mut indptr = vec![0usize];

        for i in 0..n {
            if i > 0 {
                data.push(-1.0);
                indices.push(i - 1);
            }
            data.push(2.0);
            indices.push(i);
            if i < n - 1 {
                data.push(-1.0);
                indices.push(i + 1);
            }
            indptr.push(data.len());
        }

        CsrMatrix::new(n, n, data, indices, indptr)
    }

    #[test]
    fn test_csr_matvec() {
        let csr = tridiag_csr_5x5();
        let x = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0, 5.0]);
        let mut y = Array1::zeros(5);
        csr.matvec(&x.view(), &mut y).expect("matvec failed");

        // Row 0: 2*1 - 1*2 = 0
        assert!((y[0] - 0.0).abs() < 1e-12);
        // Row 1: -1*1 + 2*2 - 1*3 = 0
        assert!((y[1] - 0.0).abs() < 1e-12);
        // Row 2: -1*2 + 2*3 - 1*4 = 0
        assert!((y[2] - 0.0).abs() < 1e-12);
        // Row 3: -1*3 + 2*4 - 1*5 = 0
        assert!((y[3] - 0.0).abs() < 1e-12);
        // Row 4: -1*4 + 2*5 = 6
        assert!((y[4] - 6.0).abs() < 1e-12);
    }

    #[test]
    fn test_csr_is_symmetric() {
        let csr = tridiag_csr_5x5();
        assert!(
            csr.is_symmetric(),
            "Tridiagonal Laplacian should be symmetric"
        );
    }

    #[test]
    fn test_csr_sparsity() {
        let csr = tridiag_csr_5x5();
        // 5x5 with 13 nonzeros: 13/25 = 0.52
        let sp = csr.sparsity();
        assert!(
            (sp - 0.52).abs() < 0.01,
            "Sparsity should be ~0.52, got {sp}"
        );
    }

    #[test]
    fn test_csrmatrix_empty() {
        let csr = CsrMatrix::<f64>::new(5, 5, vec![], vec![], vec![0, 0, 0, 0, 0, 0]);
        assert_eq!(csr.nrows(), 5);
        assert_eq!(csr.ncols(), 5);
        assert_eq!(csr.nnz(), 0);
        assert_eq!(csr.sparsity(), 0.0);
    }

    #[test]
    fn test_csr_from_dense() {
        let dense = Array2::from_shape_fn((3, 3), |(i, j)| {
            if i == j {
                2.0_f64
            } else if (i as isize - j as isize).abs() == 1 {
                -1.0
            } else {
                0.0
            }
        });
        let csr = CsrMatrix::from_dense(&dense.view(), 1e-14);
        assert_eq!(csr.nrows(), 3);
        assert_eq!(csr.ncols(), 3);
        assert_eq!(csr.nnz(), 7); // 3 diagonal + 4 off-diagonal
        assert!(csr.is_symmetric());
    }

    #[test]
    fn test_dense_to_sparse_fn() {
        let dense = Array2::<f64>::eye(4);
        let sparse = dense_to_sparse(&dense.view(), 1e-12).expect("dense_to_sparse failed");
        assert_eq!(sparse.nrows(), 4);
        assert_eq!(sparse.ncols(), 4);

        // Should be able to do matvec
        let x = Array1::from_vec(vec![1.0, 2.0, 3.0, 4.0]);
        let mut y = Array1::zeros(4);
        sparse.matvec(&x.view(), &mut y).expect("matvec failed");
        for i in 0..4 {
            assert!((y[i] - x[i]).abs() < 1e-12, "Identity matvec failed at {i}");
        }
    }

    #[test]
    #[ignore = "flaky: random initial vector causes rare convergence failure in parallel test runs"]
    fn test_lanczos_with_csr() {
        let csr = tridiag_csr_5x5();
        // The 1D Laplacian eigenvalues are in (0, 4).
        // Lanczos should be able to find the 2 largest eigenvalues.
        let result = lanczos(&csr, 2, "largest", 0.0_f64, 100, 1e-6);
        assert!(
            result.is_ok(),
            "Lanczos should succeed on tridiagonal: {:?}",
            result.as_ref().err()
        );
        let (eigenvals, _) = result.expect("already checked");
        assert_eq!(eigenvals.len(), 2);
        // Eigenvalues should be real and in reasonable range
        let max_eig = eigenvals
            .iter()
            .map(|z| z.re)
            .fold(f64::NEG_INFINITY, f64::max);
        assert!(
            max_eig > 0.0 && max_eig < 5.0,
            "Eigenvalue should be in range (0, 5), got {max_eig}"
        );
    }
}