runmat_runtime/builtins/math/linalg/factor/

lu.rs

//! MATLAB-compatible `lu` builtin with CPU-backed semantics.

use crate::builtins::common::spec::{
    BroadcastSemantics, BuiltinFusionSpec, BuiltinGpuSpec, ConstantStrategy, GpuOpKind,
    ProviderHook, ReductionNaN, ResidencyPolicy, ScalarType, ShapeRequirements,
};
use crate::builtins::common::{gpu_helpers, tensor};
use crate::builtins::math::linalg::type_resolvers::matrix_unary_type;
use crate::{build_runtime_error, BuiltinResult, RuntimeError};

use num_complex::Complex64;
use runmat_accelerate_api::{GpuTensorHandle, ProviderLuResult};
use runmat_builtins::{
    BuiltinCompletionPolicy, BuiltinDescriptor, BuiltinErrorDescriptor, BuiltinExtensionDescriptor,
    BuiltinExtensionMode, BuiltinIntegerBackendRule, BuiltinIntegerCapabilityDescriptor,
    BuiltinIntegerComputationDomain, BuiltinIntegerInputAvailability,
    BuiltinIntegerInputCapability, BuiltinIntegerOutputClassRule, BuiltinIntegerOverflowRule,
    BuiltinIntegerOverloadKind, BuiltinIntegerScalarDoubleRule, BuiltinOutputMode,
    BuiltinParamArity, BuiltinParamDescriptor, BuiltinParamType, BuiltinSignatureDescriptor,
};
use runmat_macros::runtime_builtin;
use runmat_value::{ComplexTensor, Tensor, Value};

const BUILTIN_NAME: &str = "lu";

const LU_INTEGER_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
    id: "lu-integer-input",
    mode: BuiltinExtensionMode::RunMatOnly,
    description: "lu with integer input is a RunMat extension",
    error_identifier: Some("RunMat:compatibility:LuIntegerInputExtension"),
};
const LU_LOGICAL_EXTENSION: BuiltinExtensionDescriptor = BuiltinExtensionDescriptor {
    id: "lu-logical-input",
    mode: BuiltinExtensionMode::RunMatOnly,
    description: "lu with logical input is a RunMat extension",
    error_identifier: Some("RunMat:compatibility:LuLogicalInputExtension"),
};
pub const LU_EXTENSIONS: [BuiltinExtensionDescriptor; 2] =
    [LU_INTEGER_EXTENSION, LU_LOGICAL_EXTENSION];
const LU_INTEGER_INPUTS: [BuiltinIntegerInputCapability; 1] = [BuiltinIntegerInputCapability {
    name: "A",
    classes: &crate::builtins::common::integer_capability::ALL_INTEGER_CLASSES,
    availability: BuiltinIntegerInputAvailability::RunMatOnly,
    scalar_double: BuiltinIntegerScalarDoubleRule::NotApplicable,
    notes: "RunMat mode admits integer matrices at an explicit binary64 factorization boundary.",
}];
pub const LU_INTEGER_CAPABILITIES: [BuiltinIntegerCapabilityDescriptor; 1] =
    [BuiltinIntegerCapabilityDescriptor {
        form: "[L,U,P] = lu(integer_A)",
        inputs: &LU_INTEGER_INPUTS,
        computation_domain: BuiltinIntegerComputationDomain::FloatingPoint,
        output_class: BuiltinIntegerOutputClassRule::Double,
        overflow: BuiltinIntegerOverflowRule::Error,
        backend: BuiltinIntegerBackendRule::GatherFallback,
        overload: BuiltinIntegerOverloadKind::FunctionSpecific,
        notes: "Gated RunMat extension; documented MATLAB input classes remain single and double.",
    }];

const LU_OUTPUT_COMBINED: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
    name: "LU",
    ty: BuiltinParamType::NumericArray,
    arity: BuiltinParamArity::Required,
    default: None,
    description: "Combined LU factors.",
}];

const LU_OUTPUT_LU: [BuiltinParamDescriptor; 2] = [
    BuiltinParamDescriptor {
        name: "L",
        ty: BuiltinParamType::NumericArray,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Lower-triangular factor.",
    },
    BuiltinParamDescriptor {
        name: "U",
        ty: BuiltinParamType::NumericArray,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Upper-triangular factor.",
    },
];

const LU_OUTPUT_LUP: [BuiltinParamDescriptor; 3] = [
    BuiltinParamDescriptor {
        name: "L",
        ty: BuiltinParamType::NumericArray,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Lower-triangular factor.",
    },
    BuiltinParamDescriptor {
        name: "U",
        ty: BuiltinParamType::NumericArray,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Upper-triangular factor.",
    },
    BuiltinParamDescriptor {
        name: "P",
        ty: BuiltinParamType::NumericArray,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Permutation matrix or vector based on pivot mode.",
    },
];

const LU_INPUTS_A: [BuiltinParamDescriptor; 1] = [BuiltinParamDescriptor {
    name: "A",
    ty: BuiltinParamType::NumericArray,
    arity: BuiltinParamArity::Required,
    default: None,
    description: "Input matrix to factorize.",
}];

const LU_INPUTS_A_MODE: [BuiltinParamDescriptor; 2] = [
    BuiltinParamDescriptor {
        name: "A",
        ty: BuiltinParamType::NumericArray,
        arity: BuiltinParamArity::Required,
        default: None,
        description: "Input matrix to factorize.",
    },
    BuiltinParamDescriptor {
        name: "pivotMode",
        ty: BuiltinParamType::StringScalar,
        arity: BuiltinParamArity::Required,
        default: Some("\"matrix\""),
        description: "Permutation mode (`\"matrix\"` or `\"vector\"`).",
    },
];

const LU_SIGNATURES: [BuiltinSignatureDescriptor; 6] = [
    BuiltinSignatureDescriptor {
        label: "LU = lu(A)",
        inputs: &LU_INPUTS_A,
        outputs: &LU_OUTPUT_COMBINED,
    },
    BuiltinSignatureDescriptor {
        label: "LU = lu(A, pivotMode)",
        inputs: &LU_INPUTS_A_MODE,
        outputs: &LU_OUTPUT_COMBINED,
    },
    BuiltinSignatureDescriptor {
        label: "[L, U] = lu(A)",
        inputs: &LU_INPUTS_A,
        outputs: &LU_OUTPUT_LU,
    },
    BuiltinSignatureDescriptor {
        label: "[L, U] = lu(A, pivotMode)",
        inputs: &LU_INPUTS_A_MODE,
        outputs: &LU_OUTPUT_LU,
    },
    BuiltinSignatureDescriptor {
        label: "[L, U, P] = lu(A)",
        inputs: &LU_INPUTS_A,
        outputs: &LU_OUTPUT_LUP,
    },
    BuiltinSignatureDescriptor {
        label: "[L, U, P] = lu(A, pivotMode)",
        inputs: &LU_INPUTS_A_MODE,
        outputs: &LU_OUTPUT_LUP,
    },
];

const LU_ERROR_INVALID_ARGUMENT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
    code: "RM.LU.INVALID_ARGUMENT",
    identifier: Some("RunMat:lu:InvalidArgument"),
    when: "Option arguments or requested output count are invalid.",
    message: "lu currently supports at most three outputs",
};

const LU_ERROR_INVALID_INPUT: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
    code: "RM.LU.INVALID_INPUT",
    identifier: Some("RunMat:lu:InvalidInput"),
    when: "Input is unsupported for LU factorization.",
    message: "lu: expected numeric or logical input values",
};

const LU_ERROR_INTERNAL: BuiltinErrorDescriptor = BuiltinErrorDescriptor {
    code: "RM.LU.INTERNAL",
    identifier: Some("RunMat:lu:Internal"),
    when: "Runtime cannot materialize LU outputs.",
    message: "lu: internal runtime failure",
};

const LU_ERRORS: [BuiltinErrorDescriptor; 3] = [
    LU_ERROR_INVALID_ARGUMENT,
    LU_ERROR_INVALID_INPUT,
    LU_ERROR_INTERNAL,
];

pub const LU_DESCRIPTOR: BuiltinDescriptor = BuiltinDescriptor {
    signatures: &LU_SIGNATURES,
    output_mode: BuiltinOutputMode::ByRequestedOutputCount,
    completion_policy: BuiltinCompletionPolicy::Public,
    errors: &LU_ERRORS,
};

#[runmat_macros::register_gpu_spec(builtin_path = "crate::builtins::math::linalg::factor::lu")]
pub const GPU_SPEC: BuiltinGpuSpec = BuiltinGpuSpec {
    name: "lu",
    op_kind: GpuOpKind::Custom("lu-factor"),
    supported_precisions: &[ScalarType::F32, ScalarType::F64],
    broadcast: BroadcastSemantics::None,
    provider_hooks: &[ProviderHook::Custom("lu")],
    constant_strategy: ConstantStrategy::InlineLiteral,
    residency: ResidencyPolicy::NewHandle,
    nan_mode: ReductionNaN::Include,
    two_pass_threshold: None,
    workgroup_size: None,
    accepts_nan_mode: false,
    notes: "Prefers the provider `lu` hook; automatically gathers and falls back to the CPU implementation when no provider support is registered.",
};

fn lu_error(error: &'static BuiltinErrorDescriptor) -> RuntimeError {
    lu_error_with_message(error.message, error)
}

fn lu_error_with_message(
    message: impl Into<String>,
    error: &'static BuiltinErrorDescriptor,
) -> RuntimeError {
    let mut builder = build_runtime_error(message).with_builtin(BUILTIN_NAME);
    if let Some(identifier) = error.identifier {
        builder = builder.with_identifier(identifier);
    }
    builder.build()
}

fn lu_invalid_argument(message: impl Into<String>) -> RuntimeError {
    lu_error_with_message(message, &LU_ERROR_INVALID_ARGUMENT)
}

fn lu_invalid_input(message: impl Into<String>) -> RuntimeError {
    lu_error_with_message(message, &LU_ERROR_INVALID_INPUT)
}

fn lu_internal_error(message: impl Into<String>) -> RuntimeError {
    lu_error_with_message(message, &LU_ERROR_INTERNAL)
}

fn with_lu_context(mut error: RuntimeError) -> RuntimeError {
    if error.message() == "interaction pending..." {
        return build_runtime_error("interaction pending...")
            .with_builtin(BUILTIN_NAME)
            .build();
    }
    if error.context.builtin.is_none() {
        error.context = error.context.with_builtin(BUILTIN_NAME);
    }
    error
}

#[runmat_macros::register_fusion_spec(builtin_path = "crate::builtins::math::linalg::factor::lu")]
pub const FUSION_SPEC: BuiltinFusionSpec = BuiltinFusionSpec {
    name: "lu",
    shape: ShapeRequirements::Any,
    constant_strategy: ConstantStrategy::InlineLiteral,
    elementwise: None,
    reduction: None,
    emits_nan: false,
    notes: "LU decomposition is not part of expression fusion; calls execute eagerly on the CPU.",
};

#[runtime_builtin(
    name = "lu",
    category = "math/linalg/factor",
    summary = "Compute LU decompositions with partial pivoting.",
    keywords = "lu,factorization,decomposition,permutation",
    accel = "sink",
    sink = true,
    type_resolver(matrix_unary_type),
    descriptor(crate::builtins::math::linalg::factor::lu::LU_DESCRIPTOR),
    extensions(crate::builtins::math::linalg::factor::lu::LU_EXTENSIONS),
    integer_capabilities(crate::builtins::math::linalg::factor::lu::LU_INTEGER_CAPABILITIES),
    builtin_path = "crate::builtins::math::linalg::factor::lu"
)]
async fn lu_builtin(value: Value, rest: Vec<Value>) -> crate::BuiltinResult<Value> {
    let eval = evaluate(value, &rest).await?;
    if let Some(out_count) = crate::output_count::current_output_count() {
        if out_count == 0 {
            return Ok(Value::OutputList(Vec::new()));
        }
        if out_count == 1 {
            return Ok(Value::OutputList(vec![eval.combined()]));
        }
        if out_count == 2 {
            return Ok(Value::OutputList(vec![eval.lower(), eval.upper()]));
        }
        if out_count == 3 {
            return Ok(Value::OutputList(vec![
                eval.lower(),
                eval.upper(),
                eval.permutation(),
            ]));
        }
        return Err(lu_error(&LU_ERROR_INVALID_ARGUMENT));
    }
    Ok(eval.combined())
}

/// Output form for `lu`, reused by both the builtin wrapper and the VM multi-output path.
#[derive(Clone)]
pub struct LuEval {
    combined: Value,
    lower: Value,
    upper: Value,
    perm_matrix: Value,
    perm_vector: Value,
    pivot_mode: PivotMode,
}

impl LuEval {
    /// Combined LU factor (single-output form).
    pub fn combined(&self) -> Value {
        self.combined.clone()
    }

    /// Lower-triangular factor.
    pub fn lower(&self) -> Value {
        self.lower.clone()
    }

    /// Upper-triangular factor.
    pub fn upper(&self) -> Value {
        self.upper.clone()
    }

    /// Permutation value respecting the selected pivot mode.
    pub fn permutation(&self) -> Value {
        match self.pivot_mode {
            PivotMode::Matrix => self.perm_matrix.clone(),
            PivotMode::Vector => self.perm_vector.clone(),
        }
    }

    /// Permutation matrix (always available, useful for tests).
    pub fn permutation_matrix(&self) -> Value {
        self.perm_matrix.clone()
    }

    /// Pivot vector (always available, useful for tests).
    pub fn pivot_vector(&self) -> Value {
        self.perm_vector.clone()
    }

    /// The pivot mode that was requested.
    pub fn pivot_mode(&self) -> PivotMode {
        self.pivot_mode
    }

    fn from_components(components: LuComponents, pivot_mode: PivotMode) -> BuiltinResult<Self> {
        let combined = matrix_to_value(&components.combined)?;
        let lower = matrix_to_value(&components.lower)?;
        let upper = matrix_to_value(&components.upper)?;
        let perm_matrix = matrix_to_value(&components.permutation)?;
        let perm_vector = pivot_vector_to_value(&components.pivot_vector)?;
        Ok(Self {
            combined,
            lower,
            upper,
            perm_matrix,
            perm_vector,
            pivot_mode,
        })
    }

    fn from_provider(
        mut result: ProviderLuResult,
        pivot_mode: PivotMode,
        provenance: runmat_accelerate_api::GpuHandleProvenance,
    ) -> Self {
        for handle in [
            &mut result.combined,
            &mut result.lower,
            &mut result.upper,
            &mut result.perm_matrix,
            &mut result.perm_vector,
        ] {
            runmat_accelerate_api::set_handle_provenance(handle, provenance);
            runmat_accelerate_api::mark_residency(handle);
        }
        Self {
            combined: Value::GpuTensor(result.combined),
            lower: Value::GpuTensor(result.lower),
            upper: Value::GpuTensor(result.upper),
            perm_matrix: Value::GpuTensor(result.perm_matrix),
            perm_vector: Value::GpuTensor(result.perm_vector),
            pivot_mode,
        }
    }
}

/// Permutation output mode.
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum PivotMode {
    Matrix,
    Vector,
}

impl Default for PivotMode {
    fn default() -> Self {
        Self::Matrix
    }
}

/// Evaluate `lu` while preserving all output forms for later extraction.
pub async fn evaluate(value: Value, args: &[Value]) -> BuiltinResult<LuEval> {
    let pivot_mode = parse_pivot_mode(args)?;
    ensure_lu_extensions(&value).await?;
    crate::builtins::common::validation::reject_typed_complex_integer(&value, BUILTIN_NAME)?;
    match value {
        Value::GpuTensor(handle) => {
            if let Some(eval) = evaluate_gpu(&handle, pivot_mode).await? {
                return Ok(eval);
            }
            let owner = gpu_helpers::exact_provider_for_handle(&handle);
            let explicit = runmat_accelerate_api::handle_is_explicit(&handle);
            let tensor = gpu_helpers::gather_tensor_async(&handle)
                .await
                .map_err(with_lu_context)?;
            let eval = evaluate_host_value(Value::Tensor(tensor), pivot_mode).await?;
            if explicit {
                let owner = owner.ok_or_else(|| {
                    lu_invalid_input("lu: no exact owner for explicit gpuArray input")
                })?;
                restore_lu_eval_to_provider(eval, owner)
            } else {
                Ok(eval)
            }
        }
        other => evaluate_host_value(other, pivot_mode).await,
    }
}

fn restore_lu_eval_to_provider(
    eval: LuEval,
    owner: &'static dyn runmat_accelerate_api::AccelProvider,
) -> BuiltinResult<LuEval> {
    fn upload(
        owner: &'static dyn runmat_accelerate_api::AccelProvider,
        value: &Value,
    ) -> BuiltinResult<GpuTensorHandle> {
        let handle = match value {
            Value::Tensor(tensor) => gpu_helpers::upload_tensor(owner, tensor)
                .map_err(|error| lu_internal_error(format!("lu: GPU upload failed: {error}")))?,
            Value::ComplexTensor(tensor) => gpu_helpers::upload_complex_tensor(owner, tensor)?,
            _ => return Err(lu_internal_error("lu: unexpected host factor value")),
        };
        Ok(handle)
    }
    let mut uploaded = Vec::with_capacity(5);
    for value in [
        &eval.combined,
        &eval.lower,
        &eval.upper,
        &eval.perm_matrix,
        &eval.perm_vector,
    ] {
        match upload(owner, value) {
            Ok(handle) => uploaded.push(handle),
            Err(error) => {
                for handle in &uploaded {
                    gpu_helpers::free_unprotected_exact_owner(handle, &[]);
                }
                return Err(error);
            }
        }
    }
    let mut uploaded = uploaded.into_iter();
    let mut combined = uploaded.next().expect("combined upload");
    let mut lower = uploaded.next().expect("lower upload");
    let mut upper = uploaded.next().expect("upper upload");
    let mut perm_matrix = uploaded.next().expect("permutation upload");
    let mut perm_vector = uploaded.next().expect("pivot upload");
    for handle in [
        &mut combined,
        &mut lower,
        &mut upper,
        &mut perm_matrix,
        &mut perm_vector,
    ] {
        runmat_accelerate_api::mark_handle_explicit(handle);
        runmat_accelerate_api::mark_residency(handle);
    }
    Ok(LuEval {
        combined: Value::GpuTensor(combined),
        lower: Value::GpuTensor(lower),
        upper: Value::GpuTensor(upper),
        perm_matrix: Value::GpuTensor(perm_matrix),
        perm_vector: Value::GpuTensor(perm_vector),
        pivot_mode: eval.pivot_mode,
    })
}

async fn ensure_lu_extensions(value: &Value) -> BuiltinResult<()> {
    let extension = match value {
        Value::Int(_) => Some(&LU_INTEGER_EXTENSION),
        Value::Tensor(tensor) if tensor.integer_storage().is_some() => Some(&LU_INTEGER_EXTENSION),
        Value::GpuTensor(handle)
            if runmat_accelerate_api::handle_integer_type(handle).is_some() =>
        {
            Some(&LU_INTEGER_EXTENSION)
        }
        Value::Bool(_) | Value::LogicalArray(_) => Some(&LU_LOGICAL_EXTENSION),
        Value::GpuTensor(handle) if runmat_accelerate_api::handle_is_logical(handle) => {
            Some(&LU_LOGICAL_EXTENSION)
        }
        _ => None,
    };
    if let Some(extension) = extension {
        crate::compatibility::ensure_builtin_extension_enabled(extension, BUILTIN_NAME)?;
    }
    if crate::builtins::common::validation::value_has_native_integer_class(value)
        && !crate::builtins::common::validation::native_integer_value_is_exact_f64_async(value)
            .await?
    {
        return Err(lu_invalid_input(
            "lu: integer input lies outside the exact binary64 interval",
        ));
    }
    Ok(())
}

async fn evaluate_host_value(value: Value, pivot_mode: PivotMode) -> BuiltinResult<LuEval> {
    let matrix = extract_matrix(value).await?;
    let components = lu_factor(matrix)?;
    LuEval::from_components(components, pivot_mode)
}

async fn evaluate_gpu(
    handle: &GpuTensorHandle,
    pivot_mode: PivotMode,
) -> BuiltinResult<Option<LuEval>> {
    if let Some(provider) = gpu_helpers::exact_provider_for_handle(handle) {
        if let Ok(result) = provider.lu(handle).await {
            if valid_provider_lu_result(&result, handle, provider) {
                let provenance = runmat_accelerate_api::handle_provenance(handle)
                    .unwrap_or(runmat_accelerate_api::GpuHandleProvenance::Automatic);
                return Ok(Some(LuEval::from_provider(result, pivot_mode, provenance)));
            }
            free_invalid_provider_lu_result(&result, handle);
        }
    }
    Ok(None)
}

fn valid_provider_lu_result(
    result: &ProviderLuResult,
    input: &GpuTensorHandle,
    owner: &'static dyn runmat_accelerate_api::AccelProvider,
) -> bool {
    let rows = input.shape.first().copied().unwrap_or(1);
    let outputs = [
        &result.combined,
        &result.lower,
        &result.upper,
        &result.perm_matrix,
        &result.perm_vector,
    ];
    let expected = [
        input.shape.clone(),
        vec![rows, rows],
        input.shape.clone(),
        vec![rows, rows],
        vec![rows, 1],
    ];
    outputs.iter().zip(expected.iter()).all(|(output, shape)| {
        output.shape == *shape
            && output.device_id == input.device_id
            && !gpu_helpers::same_gpu_handle(output, input)
            && runmat_accelerate_api::handle_storage(output)
                == runmat_accelerate_api::GpuTensorStorage::Real
            && runmat_accelerate_api::handle_integer_type(output).is_none()
            && !runmat_accelerate_api::handle_is_logical(output)
            && runmat_accelerate_api::handle_precision(output)
                == runmat_accelerate_api::handle_precision(input)
            && gpu_helpers::exact_provider_for_handle(output)
                .is_some_and(|candidate| std::ptr::eq(candidate, owner))
    }) && outputs.iter().enumerate().all(|(index, output)| {
        outputs
            .iter()
            .skip(index + 1)
            .all(|other| !gpu_helpers::same_gpu_handle(output, other))
    })
}

fn free_invalid_provider_lu_result(result: &ProviderLuResult, input: &GpuTensorHandle) {
    let outputs = [
        &result.combined,
        &result.lower,
        &result.upper,
        &result.perm_matrix,
        &result.perm_vector,
    ];
    for (index, output) in outputs.iter().enumerate() {
        if outputs[..index]
            .iter()
            .any(|prior| gpu_helpers::same_gpu_handle(output, prior))
        {
            continue;
        }
        gpu_helpers::free_unprotected_exact_owner(output, &[input]);
    }
}

fn parse_pivot_mode(args: &[Value]) -> BuiltinResult<PivotMode> {
    if args.is_empty() {
        return Ok(PivotMode::Matrix);
    }
    if args.len() > 1 {
        return Err(lu_invalid_argument("lu: too many option arguments"));
    }
    let Some(option) = tensor::value_to_string(&args[0]) else {
        return Err(lu_invalid_argument(
            "lu: option must be a string or character vector",
        ));
    };
    match option.trim().to_ascii_lowercase().as_str() {
        "matrix" => Ok(PivotMode::Matrix),
        "vector" => Ok(PivotMode::Vector),
        other => Err(lu_invalid_argument(format!("lu: unknown option '{other}'"))),
    }
}

async fn extract_matrix(value: Value) -> BuiltinResult<RowMajorMatrix> {
    match value {
        Value::Tensor(t) => RowMajorMatrix::from_tensor(&t),
        Value::ComplexTensor(ct) => RowMajorMatrix::from_complex_tensor(&ct),
        Value::GpuTensor(handle) => {
            let tensor = gpu_helpers::gather_tensor_async(&handle)
                .await
                .map_err(with_lu_context)?;
            RowMajorMatrix::from_tensor(&tensor)
        }
        Value::LogicalArray(logical) => {
            let tensor = tensor::logical_to_tensor(&logical)
                .map_err(|err| lu_invalid_input(format!("lu: {err}")))?;
            RowMajorMatrix::from_tensor(&tensor)
        }
        Value::Num(n) => Ok(RowMajorMatrix::from_scalar(Complex64::new(n, 0.0))),
        Value::Int(i) => Ok(RowMajorMatrix::from_scalar(Complex64::new(i.to_f64(), 0.0))),
        Value::Bool(b) => Ok(RowMajorMatrix::from_scalar(Complex64::new(
            if b { 1.0 } else { 0.0 },
            0.0,
        ))),
        Value::Complex(re, im) => Ok(RowMajorMatrix::from_scalar(Complex64::new(re, im))),
        Value::CharArray(_) | Value::String(_) | Value::StringArray(_) => Err(lu_invalid_input(
            "lu: character data is not supported; convert to numeric values first",
        )),
        other => Err(lu_invalid_input(format!(
            "lu: unsupported input type {:?}",
            other
        ))),
    }
}

struct LuComponents {
    combined: RowMajorMatrix,
    lower: RowMajorMatrix,
    upper: RowMajorMatrix,
    permutation: RowMajorMatrix,
    pivot_vector: Vec<f64>,
}

fn lu_factor(mut matrix: RowMajorMatrix) -> BuiltinResult<LuComponents> {
    let rows = matrix.rows;
    let cols = matrix.cols;
    let min_dim = rows.min(cols);
    let mut perm: Vec<usize> = (0..rows).collect();

    for k in 0..min_dim {
        // Select pivot row with maximal absolute value in column k.
        let mut pivot_row = k;
        let mut pivot_abs = 0.0;
        for r in k..rows {
            let val = matrix.get(r, k);
            let abs = val.norm();
            if abs > pivot_abs {
                pivot_abs = abs;
                pivot_row = r;
            }
        }

        if pivot_row != k {
            matrix.swap_rows(pivot_row, k);
            perm.swap(pivot_row, k);
        }

        if pivot_abs == 0.0 {
            // Entire column is effectively zero; set multipliers to zero and continue.
            for r in (k + 1)..rows {
                matrix.set(r, k, Complex64::new(0.0, 0.0));
            }
            continue;
        }

        let pivot_value = matrix.get(k, k);
        for r in (k + 1)..rows {
            let factor = matrix.get(r, k) / pivot_value;
            matrix.set(r, k, factor);
            for c in (k + 1)..cols {
                let updated = matrix.get(r, c) - factor * matrix.get(k, c);
                matrix.set(r, c, updated);
            }
        }
    }

    let combined = matrix.clone();
    let lower = build_lower(&matrix);
    let upper = build_upper(&matrix);
    let mut permutation = build_permutation(rows, &perm);
    permutation.single = matrix.single;
    let pivot_vector: Vec<f64> = perm.iter().map(|idx| (*idx + 1) as f64).collect();

    Ok(LuComponents {
        combined,
        lower,
        upper,
        permutation,
        pivot_vector,
    })
}

fn build_lower(matrix: &RowMajorMatrix) -> RowMajorMatrix {
    let rows = matrix.rows;
    let cols = matrix.cols;
    let min_dim = rows.min(cols);
    let mut lower = RowMajorMatrix::identity(rows);
    lower.single = matrix.single;
    for i in 0..rows {
        for j in 0..min_dim {
            if i > j {
                lower.set(i, j, matrix.get(i, j));
            }
        }
    }
    lower
}

fn build_upper(matrix: &RowMajorMatrix) -> RowMajorMatrix {
    let rows = matrix.rows;
    let cols = matrix.cols;
    let mut upper = RowMajorMatrix::zeros(rows, cols);
    upper.single = matrix.single;
    for i in 0..rows {
        for j in 0..cols {
            if i <= j {
                upper.set(i, j, matrix.get(i, j));
            }
        }
    }
    upper
}

fn build_permutation(rows: usize, perm: &[usize]) -> RowMajorMatrix {
    let mut matrix = RowMajorMatrix::zeros(rows, rows);
    for (i, &col) in perm.iter().enumerate() {
        if col < rows {
            matrix.set(i, col, Complex64::new(1.0, 0.0));
        }
    }
    matrix
}

const EPS: f64 = 1.0e-12;

fn matrix_to_value(matrix: &RowMajorMatrix) -> BuiltinResult<Value> {
    let mut has_imag = false;
    for val in &matrix.data {
        if val.im.abs() > EPS {
            has_imag = true;
            break;
        }
    }
    if has_imag {
        let mut data = Vec::with_capacity(matrix.rows * matrix.cols);
        for col in 0..matrix.cols {
            for row in 0..matrix.rows {
                let idx = row * matrix.cols + col;
                let v = matrix.data[idx];
                data.push((v.re, v.im));
            }
        }
        let tensor = if matrix.single {
            ComplexTensor::from_f32(
                data.into_iter()
                    .map(|(re, im)| (re as f32, im as f32))
                    .collect(),
                vec![matrix.rows, matrix.cols],
            )
        } else {
            ComplexTensor::new(data, vec![matrix.rows, matrix.cols])
        }
        .map_err(|e| lu_internal_error(format!("lu: {e}")))?;
        Ok(Value::ComplexTensor(tensor))
    } else {
        let mut data = Vec::with_capacity(matrix.rows * matrix.cols);
        for col in 0..matrix.cols {
            for row in 0..matrix.rows {
                let idx = row * matrix.cols + col;
                data.push(matrix.data[idx].re);
            }
        }
        let tensor = if matrix.single {
            Tensor::from_f32(
                data.into_iter().map(|value| value as f32).collect(),
                vec![matrix.rows, matrix.cols],
            )
        } else {
            Tensor::new(data, vec![matrix.rows, matrix.cols])
        }
        .map_err(|e| lu_internal_error(format!("lu: {e}")))?;
        Ok(Value::Tensor(tensor))
    }
}

fn pivot_vector_to_value(pivot: &[f64]) -> BuiltinResult<Value> {
    let rows = pivot.len();
    let tensor = Tensor::new(pivot.to_vec(), vec![rows, 1])
        .map_err(|e| lu_internal_error(format!("lu: {e}")))?;
    Ok(Value::Tensor(tensor))
}

#[derive(Clone)]
struct RowMajorMatrix {
    rows: usize,
    cols: usize,
    data: Vec<Complex64>,
    single: bool,
}

impl RowMajorMatrix {
    fn zeros(rows: usize, cols: usize) -> Self {
        Self {
            rows,
            cols,
            data: vec![Complex64::new(0.0, 0.0); rows.saturating_mul(cols)],
            single: false,
        }
    }

    fn identity(size: usize) -> Self {
        let mut matrix = Self::zeros(size, size);
        for i in 0..size {
            matrix.set(i, i, Complex64::new(1.0, 0.0));
        }
        matrix
    }

    fn from_scalar(value: Complex64) -> Self {
        Self {
            rows: 1,
            cols: 1,
            data: vec![value],
            single: false,
        }
    }

    fn from_tensor(tensor: &Tensor) -> BuiltinResult<Self> {
        if tensor.shape.len() > 2 {
            return Err(lu_invalid_input("lu: input must be 2-D"));
        }
        let rows = tensor.rows();
        let cols = tensor.cols();
        let values = tensor::tensor_values_f64_cow(tensor);
        let mut data = vec![Complex64::new(0.0, 0.0); rows.saturating_mul(cols)];
        for col in 0..cols {
            for row in 0..rows {
                let idx_col_major = row + col * rows;
                let idx_row_major = row * cols + col;
                data[idx_row_major] = Complex64::new(values[idx_col_major], 0.0);
            }
        }
        Ok(Self {
            rows,
            cols,
            data,
            single: tensor.numeric_dtype() == runmat_value::NumericDType::F32,
        })
    }

    fn from_complex_tensor(tensor: &ComplexTensor) -> BuiltinResult<Self> {
        if tensor.shape.len() > 2 {
            return Err(lu_invalid_input("lu: input must be 2-D"));
        }
        let rows = tensor.rows;
        let cols = tensor.cols;
        let mut data = vec![Complex64::new(0.0, 0.0); rows.saturating_mul(cols)];
        for col in 0..cols {
            for row in 0..rows {
                let idx_col_major = row + col * rows;
                let idx_row_major = row * cols + col;
                let (re, im) = tensor.materialize_f64()[idx_col_major];
                data[idx_row_major] = Complex64::new(re, im);
            }
        }
        Ok(Self {
            rows,
            cols,
            data,
            single: tensor.numeric_dtype() == runmat_value::NumericDType::F32,
        })
    }

    fn get(&self, row: usize, col: usize) -> Complex64 {
        self.data[row * self.cols + col]
    }

    fn set(&mut self, row: usize, col: usize, value: Complex64) {
        self.data[row * self.cols + col] = value;
    }

    fn swap_rows(&mut self, r1: usize, r2: usize) {
        if r1 == r2 {
            return;
        }
        for col in 0..self.cols {
            self.data.swap(r1 * self.cols + col, r2 * self.cols + col);
        }
    }
}

#[cfg(test)]
pub(crate) mod tests {
    use super::*;
    use crate::builtins::common::test_support;
    use futures::executor::block_on;
    #[cfg(feature = "wgpu")]
    use runmat_accelerate_api::AccelProvider;
    use runmat_builtins::{ResolveContext, Type};
    use runmat_value::{ComplexTensor as CMatrix, IntegerStorage, Tensor as Matrix};

    fn error_message(err: RuntimeError) -> String {
        err.message().to_string()
    }

    fn tensor_from_value(value: Value) -> Matrix {
        match value {
            Value::Tensor(t) => t,
            other => panic!("expected dense tensor, got {other:?}"),
        }
    }

    fn row_major_from_value(value: Value) -> RowMajorMatrix {
        match value {
            Value::Tensor(t) => RowMajorMatrix::from_tensor(&t).expect("row-major tensor"),
            Value::ComplexTensor(ct) => {
                RowMajorMatrix::from_complex_tensor(&ct).expect("row-major complex tensor")
            }
            other => panic!("expected tensor value, got {other:?}"),
        }
    }

    #[test]
    fn lu_type_preserves_matrix_shape() {
        let out = matrix_unary_type(
            &[Type::Tensor {
                shape: Some(vec![Some(2), Some(3)]),
            }],
            &ResolveContext::new(Vec::new()),
        );
        assert_eq!(
            out,
            Type::Tensor {
                shape: Some(vec![Some(2), Some(3)])
            }
        );
    }

    #[test]
    fn lu_descriptor_signatures_cover_core_forms() {
        let labels: Vec<&str> = LU_DESCRIPTOR
            .signatures
            .iter()
            .map(|signature| signature.label)
            .collect();
        assert!(labels.contains(&"LU = lu(A)"));
        assert!(labels.contains(&"LU = lu(A, pivotMode)"));
        assert!(labels.contains(&"[L, U] = lu(A)"));
        assert!(labels.contains(&"[L, U] = lu(A, pivotMode)"));
        assert!(labels.contains(&"[L, U, P] = lu(A)"));
        assert!(labels.contains(&"[L, U, P] = lu(A, pivotMode)"));
    }

    #[test]
    fn lu_descriptor_errors_have_stable_codes() {
        let codes: Vec<&str> = LU_DESCRIPTOR.errors.iter().map(|err| err.code).collect();
        assert!(codes.contains(&"RM.LU.INVALID_ARGUMENT"));
        assert!(codes.contains(&"RM.LU.INVALID_INPUT"));
        assert!(codes.contains(&"RM.LU.INTERNAL"));
    }

    #[test]
    fn lu_matrix_conversion_reads_typed_integer_storage_exactly() {
        let tensor = Matrix::new_integer(IntegerStorage::I16(vec![4, 6, 3, 3]), vec![2, 2])
            .expect("typed integer tensor");

        let matrix = RowMajorMatrix::from_tensor(&tensor).expect("matrix");
        assert_eq!(matrix.rows, 2);
        assert_eq!(matrix.cols, 2);
        assert_eq!(
            matrix.data,
            vec![
                Complex64::new(4.0, 0.0),
                Complex64::new(3.0, 0.0),
                Complex64::new(6.0, 0.0),
                Complex64::new(3.0, 0.0),
            ]
        );
    }

    #[test]
    fn lu_compatibility_mode_gates_integer_and_logical_extensions() {
        let _compat = crate::compatibility::push_runmat_extensions_enabled(false);
        let integer = Matrix::new_integer(IntegerStorage::I16(vec![1]), vec![1, 1]).unwrap();
        let error = match evaluate(Value::Tensor(integer), &[]) {
            Err(error) => error,
            Ok(_) => panic!("integer extension must be gated"),
        };
        assert_eq!(
            error.identifier(),
            Some("RunMat:compatibility:LuIntegerInputExtension")
        );
        let error = match evaluate(Value::Bool(true), &[]) {
            Err(error) => error,
            Ok(_) => panic!("logical extension must be gated"),
        };
        assert_eq!(
            error.identifier(),
            Some("RunMat:compatibility:LuLogicalInputExtension")
        );
    }

    #[test]
    fn lu_rejects_integer_values_outside_exact_binary64_interval() {
        let _extensions = crate::compatibility::push_runmat_extensions_enabled(true);
        let integer =
            Matrix::new_integer(IntegerStorage::U64(vec![9_007_199_254_740_993]), vec![1, 1])
                .unwrap();
        let error = match evaluate(Value::Tensor(integer), &[]) {
            Err(error) => error,
            Ok(_) => panic!("wide integer must reject"),
        };
        assert_eq!(error.identifier(), LU_ERROR_INVALID_INPUT.identifier);
    }

    #[test]
    fn lu_preserves_single_precision_outputs() {
        let input = Matrix::from_f32(vec![2.0, 1.0, 1.0, 2.0], vec![2, 2]).unwrap();
        let eval = evaluate(Value::Tensor(input), &[]).expect("single LU");
        for output in [
            eval.combined(),
            eval.lower(),
            eval.upper(),
            eval.permutation_matrix(),
        ] {
            let Value::Tensor(tensor) = output else {
                panic!("expected real tensor");
            };
            assert_eq!(tensor.numeric_dtype(), runmat_value::NumericDType::F32);
        }
    }

    #[test]
    fn lu_does_not_treat_small_nonzero_pivot_as_zero() {
        let input = Matrix::new(vec![1.0e-14, 0.0, 0.0, 2.0e-14], vec![2, 2]).unwrap();
        let eval = evaluate(Value::Tensor(input.clone()), &[]).expect("small-scale LU");
        let l = tensor_from_value(eval.lower());
        let u = tensor_from_value(eval.upper());
        let p = tensor_from_value(eval.permutation_matrix());
        let pa = crate::builtins::common::matrix::matrix_mul(&p, &input).unwrap();
        let product = crate::builtins::common::matrix::matrix_mul(&l, &u).unwrap();
        assert_tensor_close(&pa, &product, 1e-28);
    }

    fn row_major_matmul(a: &RowMajorMatrix, b: &RowMajorMatrix) -> RowMajorMatrix {
        assert_eq!(a.cols, b.rows, "incompatible shapes for matmul");
        let mut out = RowMajorMatrix::zeros(a.rows, b.cols);
        for i in 0..a.rows {
            for k in 0..a.cols {
                let aik = a.get(i, k);
                for j in 0..b.cols {
                    let acc = out.get(i, j) + aik * b.get(k, j);
                    out.set(i, j, acc);
                }
            }
        }
        out
    }

    fn assert_tensor_close(a: &Matrix, b: &Matrix, tol: f64) {
        assert_eq!(a.shape, b.shape);
        for (lhs, rhs) in a.materialize_f64().iter().zip(&b.materialize_f64()) {
            assert!(
                (lhs - rhs).abs() <= tol,
                "mismatch: lhs={lhs}, rhs={rhs}, tol={tol}"
            );
        }
    }

    fn assert_row_major_close(a: &RowMajorMatrix, b: &RowMajorMatrix, tol: f64) {
        assert_eq!(a.rows, b.rows, "row mismatch");
        assert_eq!(a.cols, b.cols, "col mismatch");
        for row in 0..a.rows {
            for col in 0..a.cols {
                let lhs = a.get(row, col);
                let rhs = b.get(row, col);
                let diff = (lhs - rhs).norm();
                assert!(
                    diff <= tol,
                    "mismatch at ({row}, {col}): lhs={lhs:?}, rhs={rhs:?}, diff={diff}, tol={tol}"
                );
            }
        }
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_single_output_produces_combined_matrix() {
        let a = Matrix::new(
            vec![2.0, 4.0, -2.0, 1.0, -6.0, 7.0, 1.0, 0.0, 2.0],
            vec![3, 3],
        )
        .unwrap();
        let result = lu_builtin(Value::Tensor(a.clone()), Vec::new()).expect("lu");
        let lu = tensor_from_value(result);
        let eval = evaluate(Value::Tensor(a), &[]).expect("evaluate");
        let expected = tensor_from_value(eval.combined());
        assert_tensor_close(&lu, &expected, 1e-12);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_three_outputs_matches_factorization() {
        let data = vec![2.0, 4.0, -2.0, 1.0, -6.0, 7.0, 1.0, 0.0, 2.0];
        let a = Matrix::new(data.clone(), vec![3, 3]).unwrap();
        let eval = evaluate(Value::Tensor(a.clone()), &[]).expect("evaluate");
        let l = tensor_from_value(eval.lower());
        let u = tensor_from_value(eval.upper());
        let p = tensor_from_value(eval.permutation_matrix());

        let pa = crate::builtins::common::matrix::matrix_mul(&p, &a).expect("P*A");
        let lu_product = crate::builtins::common::matrix::matrix_mul(&l, &u).expect("L*U");
        assert_tensor_close(&pa, &lu_product, 1e-9);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_complex_matrix_factorization() {
        let data = vec![(1.0, 2.0), (3.0, -1.0), (2.0, -1.0), (4.0, 2.0)];
        let a = CMatrix::new(data.clone(), vec![2, 2]).expect("complex tensor");
        let eval = evaluate(Value::ComplexTensor(a.clone()), &[]).expect("evaluate complex");

        let l = row_major_from_value(eval.lower());
        let u = row_major_from_value(eval.upper());
        let p = row_major_from_value(eval.permutation_matrix());
        let input = RowMajorMatrix::from_complex_tensor(&a).expect("row-major input");

        let pa = row_major_matmul(&p, &input);
        let lu = row_major_matmul(&l, &u);
        assert_row_major_close(&pa, &lu, 1e-9);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_handles_singular_matrix() {
        let a = Matrix::new(vec![0.0, 0.0, 0.0, 0.0], vec![2, 2]).unwrap();
        let eval = evaluate(Value::Tensor(a.clone()), &[]).expect("evaluate singular");
        let l = tensor_from_value(eval.lower());
        let u = tensor_from_value(eval.upper());
        let p = tensor_from_value(eval.permutation_matrix());

        assert!(u.materialize_f64().iter().any(|&v| v.abs() <= 1e-12));

        let pa = crate::builtins::common::matrix::matrix_mul(&p, &a).expect("P*A");
        let lu_product = crate::builtins::common::matrix::matrix_mul(&l, &u).expect("L*U");
        assert_tensor_close(&pa, &lu_product, 1e-9);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_vector_option_returns_pivot_vector() {
        let a = Matrix::new(vec![4.0, 6.0, 3.0, 3.0], vec![2, 2]).unwrap();
        let eval =
            evaluate(Value::Tensor(a), &[Value::from("vector")]).expect("evaluate vector mode");
        assert_eq!(eval.pivot_mode(), PivotMode::Vector);
        let pivot = tensor_from_value(eval.pivot_vector());
        assert_eq!(pivot.shape, vec![2, 1]);
        assert_eq!(pivot.materialize_f64(), vec![2.0, 1.0]);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_vector_option_case_insensitive() {
        let a = Matrix::new(vec![4.0, 6.0, 3.0, 3.0], vec![2, 2]).unwrap();
        let eval =
            evaluate(Value::Tensor(a), &[Value::from("VECTOR")]).expect("evaluate vector option");
        assert_eq!(eval.pivot_mode(), PivotMode::Vector);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_matrix_option_returns_permutation_matrix() {
        let a = Matrix::new(vec![2.0, 1.0, 3.0, 4.0], vec![2, 2]).unwrap();
        let eval =
            evaluate(Value::Tensor(a), &[Value::from("matrix")]).expect("evaluate matrix option");
        assert_eq!(eval.pivot_mode(), PivotMode::Matrix);
        let perm_selected = tensor_from_value(eval.permutation());
        let perm_matrix = tensor_from_value(eval.permutation_matrix());
        assert_eq!(perm_selected.shape, perm_matrix.shape);
        assert_tensor_close(&perm_selected, &perm_matrix, 1e-12);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_handles_rectangular_matrices() {
        let a = Matrix::new(vec![3.0, 6.0, 1.0, 3.0, 2.0, 4.0], vec![2, 3]).unwrap();
        let eval = evaluate(Value::Tensor(a.clone()), &[]).expect("evaluate rectangular");
        let l = tensor_from_value(eval.lower());
        let u = tensor_from_value(eval.upper());
        let p = tensor_from_value(eval.permutation_matrix());
        assert_eq!(l.shape, vec![2, 2]);
        assert_eq!(u.shape, vec![2, 3]);
        assert_eq!(p.shape, vec![2, 2]);

        let pa = crate::builtins::common::matrix::matrix_mul(&p, &a).expect("P*A");
        let lu_product = crate::builtins::common::matrix::matrix_mul(&l, &u).expect("L*U");
        assert_tensor_close(&pa, &lu_product, 1e-9);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_rejects_unknown_option() {
        let a = Matrix::new(vec![1.0], vec![1, 1]).unwrap();
        let err = match evaluate(Value::Tensor(a), &[Value::from("invalid")]) {
            Ok(_) => panic!("expected option parse failure"),
            Err(err) => {
                assert_eq!(err.identifier(), LU_ERROR_INVALID_ARGUMENT.identifier);
                error_message(err)
            }
        };
        assert!(err.contains("unknown option"));
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_rejects_non_string_option() {
        let a = Matrix::new(vec![1.0], vec![1, 1]).unwrap();
        let err = match evaluate(Value::Tensor(a), &[Value::Num(2.0)]) {
            Ok(_) => panic!("expected option parse failure"),
            Err(err) => {
                assert_eq!(err.identifier(), LU_ERROR_INVALID_ARGUMENT.identifier);
                error_message(err)
            }
        };
        assert!(err.contains("unknown option"));
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_rejects_multiple_options() {
        let a = Matrix::new(vec![1.0], vec![1, 1]).unwrap();
        let err = match evaluate(
            Value::Tensor(a),
            &[Value::from("matrix"), Value::from("vector")],
        ) {
            Ok(_) => panic!("expected option arity failure"),
            Err(err) => {
                assert_eq!(err.identifier(), LU_ERROR_INVALID_ARGUMENT.identifier);
                error_message(err)
            }
        };
        assert!(err.contains("too many option arguments"));
    }

    #[test]
    fn lu_invalid_input_identifier_is_stable() {
        let tensor = Matrix::new(vec![1.0, 2.0, 3.0, 4.0], vec![1, 2, 2]).expect("tensor");
        let err = match evaluate(Value::Tensor(tensor), &[]) {
            Ok(_) => panic!("expected 2-D input failure"),
            Err(err) => err,
        };
        assert_eq!(err.identifier(), LU_ERROR_INVALID_INPUT.identifier);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_gpu_provider_roundtrip() {
        test_support::with_test_provider(|provider| {
            let host = Matrix::new(vec![10.0, 3.0, 7.0, 2.0], vec![2, 2]).unwrap();
            let view = runmat_accelerate_api::HostTensorView {
                data: &host.materialize_f64(),
                shape: &host.shape,
            };
            let handle = provider.upload(&view).expect("upload");
            let eval = evaluate(Value::GpuTensor(handle.clone()), &[]).expect("evaluate gpu input");
            let lower_val = eval.lower();
            let upper_val = eval.upper();
            let perm_val = eval.permutation_matrix();
            assert!(matches!(lower_val, Value::GpuTensor(_)));
            assert!(matches!(upper_val, Value::GpuTensor(_)));
            assert!(matches!(perm_val, Value::GpuTensor(_)));
            let l = test_support::gather(lower_val).expect("gather lower");
            let u = test_support::gather(upper_val).expect("gather upper");
            let p = test_support::gather(perm_val).expect("gather permutation");
            let pa = crate::builtins::common::matrix::matrix_mul(&p, &host).expect("P*A");
            let lu_product = crate::builtins::common::matrix::matrix_mul(&l, &u).expect("L*U");
            assert_tensor_close(&pa, &lu_product, 1e-9);
        });
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_gpu_vector_option_roundtrip() {
        test_support::with_test_provider(|provider| {
            let host = Matrix::new(vec![4.0, 6.0, 3.0, 3.0], vec![2, 2]).unwrap();
            let view = runmat_accelerate_api::HostTensorView {
                data: &host.materialize_f64(),
                shape: &host.shape,
            };
            let handle = provider.upload(&view).expect("upload");
            let eval =
                evaluate(Value::GpuTensor(handle), &[Value::from("vector")]).expect("gpu vector");
            let pivot_val = eval.permutation();
            assert!(matches!(pivot_val, Value::GpuTensor(_)));
            let pivot = test_support::gather(pivot_val).expect("gather pivot");
            assert_eq!(pivot.shape, vec![2, 1]);
            let expected = Matrix::new(vec![2.0, 1.0], vec![2, 1]).unwrap();
            assert_tensor_close(&pivot, &expected, 1e-12);
        });
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    fn lu_accepts_scalar_inputs() {
        let eval = evaluate(Value::Num(5.0), &[]).expect("evaluate scalar");
        let l = tensor_from_value(eval.lower());
        let u = tensor_from_value(eval.upper());
        let p = tensor_from_value(eval.permutation_matrix());
        assert_eq!(l.materialize_f64(), vec![1.0]);
        assert_eq!(u.materialize_f64(), vec![5.0]);
        assert_eq!(p.materialize_f64(), vec![1.0]);
    }

    #[cfg_attr(target_arch = "wasm32", wasm_bindgen_test::wasm_bindgen_test)]
    #[test]
    #[cfg(feature = "wgpu")]
    fn lu_wgpu_matches_cpu() {
        let Ok(provider) = runmat_accelerate::backend::wgpu::provider::register_wgpu_provider(
            runmat_accelerate::backend::wgpu::provider::WgpuProviderOptions::default(),
        ) else {
            return;
        };
        let host = Matrix::new(
            vec![2.0, 4.0, -2.0, 1.0, -6.0, 7.0, 1.0, 0.0, 2.0],
            vec![3, 3],
        )
        .unwrap();
        let cpu_eval = evaluate(Value::Tensor(host.clone()), &[]).expect("cpu evaluate");
        let view = runmat_accelerate_api::HostTensorView {
            data: &host.materialize_f64(),
            shape: &host.shape,
        };
        let handle = provider.upload(&view).expect("upload");
        let gpu_eval = evaluate(Value::GpuTensor(handle), &[]).expect("gpu evaluate");

        let l_cpu = tensor_from_value(cpu_eval.lower());
        let u_cpu = tensor_from_value(cpu_eval.upper());
        let p_cpu = tensor_from_value(cpu_eval.permutation_matrix());
        let lu_cpu = tensor_from_value(cpu_eval.combined());

        let l_gpu = test_support::gather(gpu_eval.lower()).expect("gather L");
        let u_gpu = test_support::gather(gpu_eval.upper()).expect("gather U");
        let p_gpu = test_support::gather(gpu_eval.permutation_matrix()).expect("gather P");
        let lu_gpu = test_support::gather(gpu_eval.combined()).expect("gather LU");

        assert_tensor_close(&l_cpu, &l_gpu, 1e-12);
        assert_tensor_close(&u_cpu, &u_gpu, 1e-12);
        assert_tensor_close(&p_cpu, &p_gpu, 1e-12);
        assert_tensor_close(&lu_cpu, &lu_gpu, 1e-12);

        let pivot_cpu = tensor_from_value(cpu_eval.pivot_vector());
        let pivot_gpu = test_support::gather(gpu_eval.pivot_vector()).expect("gather pivot vector");
        assert_tensor_close(&pivot_cpu, &pivot_gpu, 1e-12);

        let handle_vector = provider.upload(&view).expect("upload vector option");
        let gpu_vector_eval = evaluate(Value::GpuTensor(handle_vector), &[Value::from("vector")])
            .expect("gpu vector evaluate");
        let pivot_vector =
            test_support::gather(gpu_vector_eval.permutation()).expect("gather vector pivot");
        assert_tensor_close(&pivot_cpu, &pivot_vector, 1e-12);
    }

    fn lu_builtin(value: Value, rest: Vec<Value>) -> BuiltinResult<Value> {
        block_on(super::lu_builtin(value, rest))
    }

    fn evaluate(value: Value, args: &[Value]) -> BuiltinResult<LuEval> {
        block_on(super::evaluate(value, args))
    }
}