Enhanced BLAS/LAPACK return code interpretation and logging

docs/src/basics/common_solver_opts.md

@@ -108,10 +108,62 @@ sol = solve(prob; verbose=verbose)
 #### Verbosity Levels
 ##### Error Control Settings
 - default_lu_fallback: Controls messages when falling back to LU factorization (default: Warn)
+- blas_invalid_args: Controls messages when BLAS/LAPACK receives invalid arguments (default: Error)
 ##### Performance Settings
 - no_right_preconditioning: Controls messages when right preconditioning is not used (default: Warn)
 ##### Numerical Settings
 - using_IterativeSolvers: Controls messages when using the IterativeSolvers.jl package (default: Warn)
 - IterativeSolvers_iterations: Controls messages about iteration counts from IterativeSolvers.jl (default: Warn)
 - KrylovKit_verbosity: Controls messages from the KrylovKit.jl package (default: Warn)
-- KrylovJL_verbosity: Controls verbosity of the KrylovJL.jl package (default: None)
+- KrylovJL_verbosity: Controls verbosity of the KrylovJL.jl package (default: None)
+- blas_errors: Controls messages for BLAS/LAPACK errors like singular matrices (default: Warn)
+- blas_info: Controls informational messages from BLAS/LAPACK operations (default: None)
+- blas_success: Controls success messages from BLAS/LAPACK operations (default: None)
+- condition_number: Controls computation and logging of matrix condition numbers (default: None)
+
+### BLAS/LAPACK Return Code Interpretation
+
+LinearSolve.jl now provides detailed interpretation of BLAS/LAPACK return codes (info parameter) to help diagnose numerical issues. When BLAS/LAPACK operations encounter problems, the logging system will provide human-readable explanations based on the specific return code and operation.
+
+#### Common BLAS/LAPACK Return Codes
+
+- **info = 0**: Operation completed successfully
+- **info < 0**: Argument -info had an illegal value (e.g., wrong dimensions, invalid parameters)
+- **info > 0**: Operation-specific issues:
+  - **LU factorization (getrf)**: U(info,info) is exactly zero, matrix is singular
+  - **Cholesky factorization (potrf)**: Leading minor of order info is not positive definite
+  - **QR factorization (geqrf)**: Numerical issues with Householder reflector info
+  - **SVD (gesdd/gesvd)**: Algorithm failed to converge, info off-diagonal elements did not converge
+  - **Eigenvalue computation (syev/heev)**: info off-diagonal elements did not converge to zero
+  - **Bunch-Kaufman (sytrf/hetrf)**: D(info,info) is exactly zero, matrix is singular
+
+#### Example Usage with Enhanced BLAS Logging
+
+```julia
+using LinearSolve
+
+# Enable detailed BLAS error logging with condition numbers
+verbose = LinearVerbosity(
+    blas_errors = Verbosity.Info(),     # Show detailed error interpretations
+    blas_info = Verbosity.Info(),       # Show all BLAS operation details
+    condition_number = Verbosity.Info() # Compute and log condition numbers
+)
+
+# Solve a potentially singular system
+A = [1.0 2.0; 2.0 4.0]  # Singular matrix
+b = [1.0, 2.0]
+prob = LinearProblem(A, b)
+sol = solve(prob, LUFactorization(); verbose=verbose)
+
+# The logging will show:
+# - BLAS/LAPACK dgetrf: Matrix is singular
+# - Details: U(2,2) is exactly zero. The factorization has been completed, but U is singular
+# - Return code (info): 2
+# - Additional context: matrix size, condition number, memory usage
+```
+
+The enhanced logging also provides:
+- Matrix properties (size, type, element type)
+- Memory usage estimates
+- Detailed context for debugging numerical issues
+- Condition number computation controlled by the condition_number verbosity setting

ext/LinearSolveBLISExt.jl

@@ -9,8 +9,11 @@ using LinearSolve
 using LinearAlgebra: BlasInt, LU
 using LinearAlgebra.LAPACK: require_one_based_indexing, chkfinite, chkstride1, 
                             @blasfunc, chkargsok
-using LinearSolve: ArrayInterface, BLISLUFactorization, @get_cacheval, LinearCache, SciMLBase
+using LinearSolve: ArrayInterface, BLISLUFactorization, @get_cacheval, LinearCache, SciMLBase,
+                   interpret_blas_code, log_blas_info, get_blas_operation_info, 
+                   check_and_log_lapack_result, LinearVerbosity
 using SciMLBase: ReturnCode
+using SciMLLogging: Verbosity, @SciMLMessage
 
 const global libblis = blis_jll.blis
 const global liblapack = LAPACK_jll.liblapack
@@ -220,11 +223,31 @@ function SciMLBase.solve!(cache::LinearCache, alg::BLISLUFactorization;
     kwargs...)
     A = cache.A
     A = convert(AbstractMatrix, A)
+    verbose = cache.verbose
+
     if cache.isfresh
         cacheval = @get_cacheval(cache, :BLISLUFactorization)
+
+        # Perform the factorization
         res = getrf!(A; ipiv = cacheval[1].ipiv, info = cacheval[2])
+
         fact = LU(res[1:3]...), res[4]
         cache.cacheval = fact
+
+        # Log BLAS return code with detailed interpretation if logging is enabled
+        info_value = res[3]
+        if info_value != 0
+            # Only get operation info if we need to log
+            if verbose.numerical.blas_errors !== Verbosity.None()
+                op_info = get_blas_operation_info(:dgetrf, A, cache.b, verbose)
+                log_blas_info(:dgetrf, info_value, verbose; extra_context=op_info)
+            end
+        elseif verbose.numerical.blas_success !== Verbosity.None()
+            # Only get operation info if we need to log success
+            op_info = get_blas_operation_info(:dgetrf, A, cache.b, verbose)
+            success_msg = "BLAS LU factorization (dgetrf) completed successfully for $(op_info[:matrix_size]) matrix"
+            @SciMLMessage(success_msg, verbose.numerical.blas_success, :blas_success, :numerical)
+        end
 
         if !LinearAlgebra.issuccess(fact[1])
             return SciMLBase.build_linear_solution(
@@ -236,6 +259,8 @@ function SciMLBase.solve!(cache::LinearCache, alg::BLISLUFactorization;
     A, info = @get_cacheval(cache, :BLISLUFactorization)
     require_one_based_indexing(cache.u, cache.b)
     m, n = size(A, 1), size(A, 2)
+
+    # Perform the solve
     if m > n
         Bc = copy(cache.b)
         getrs!('N', A.factors, A.ipiv, Bc; info)
@@ -244,6 +269,11 @@ function SciMLBase.solve!(cache::LinearCache, alg::BLISLUFactorization;
         copyto!(cache.u, cache.b)
         getrs!('N', A.factors, A.ipiv, cache.u; info)
     end
+
+    # Log solve operation result if there was an error
+    if info[] != 0 && verbose.numerical.blas_errors !== Verbosity.None()
+        log_blas_info(:dgetrs, info[], verbose)
+    end
 
     SciMLBase.build_linear_solution(alg, cache.u, nothing, cache; retcode = ReturnCode.Success)
 end

src/LinearSolve.jl

@@ -339,6 +339,7 @@ const BLASELTYPES = Union{Float32, Float64, ComplexF32, ComplexF64}
 function defaultalg_symbol end
 
 include("verbosity.jl")
+include("blas_logging.jl")
 include("generic_lufact.jl")
 include("common.jl")
 include("extension_algs.jl")

src/blas_logging.jl

@@ -0,0 +1,213 @@
+# BLAS and LAPACK Return Code Interpretation
+
+using SciMLLogging: Verbosity, @match, @SciMLMessage
+using LinearAlgebra: cond
+
+"""
+    interpret_blas_code(func::Symbol, info::Integer)
+
+Interpret BLAS/LAPACK return codes (info parameter) to provide human-readable error messages.
+Returns a tuple of (category::Symbol, message::String, details::String)
+"""
+function interpret_blas_code(func::Symbol, info::Integer)
+    if info == 0
+        return (:success, "Operation completed successfully", "")
+    elseif info < 0
+        return (:invalid_argument, 
+                "Invalid argument error", 
+                "Argument $(-info) had an illegal value")
+    else
+        # info > 0 means different things for different functions
+        return interpret_positive_info(func, info)
+    end
+end
+
+function interpret_positive_info(func::Symbol, info::Integer)
+    func_str = string(func)
+
+    # LU factorization routines
+    if occursin("getrf", func_str)
+        return (:singular_matrix,
+                "Matrix is singular",
+                "U($info,$info) is exactly zero. The factorization has been completed, but U is singular and division by U will produce infinity.")
+
+    # Cholesky factorization routines
+    elseif occursin("potrf", func_str)
+        return (:not_positive_definite,
+                "Matrix is not positive definite",
+                "The leading minor of order $info is not positive definite, and the factorization could not be completed.")
+
+    # QR factorization routines
+    elseif occursin("geqrf", func_str) || occursin("geqrt", func_str)
+        return (:numerical_issue,
+                "Numerical issue in QR factorization",
+                "Householder reflector $info could not be formed properly.")
+
+    # SVD routines
+    elseif occursin("gesdd", func_str) || occursin("gesvd", func_str)
+        return (:convergence_failure,
+                "SVD did not converge",
+                "The algorithm failed to compute singular values. $info off-diagonal elements of an intermediate bidiagonal form did not converge to zero.")
+
+    # Symmetric/Hermitian eigenvalue routines
+    elseif occursin("syev", func_str) || occursin("heev", func_str)
+        return (:convergence_failure,
+                "Eigenvalue computation did not converge",
+                "$info off-diagonal elements of an intermediate tridiagonal form did not converge to zero.")
+
+    # Bunch-Kaufman factorization
+    elseif occursin("sytrf", func_str) || occursin("hetrf", func_str)
+        return (:singular_matrix,
+                "Matrix is singular", 
+                "D($info,$info) is exactly zero. The factorization has been completed, but the block diagonal matrix D is singular.")
+
+    # Solve routines (should not have positive info)
+    elseif occursin("getrs", func_str) || occursin("potrs", func_str) || 
+           occursin("sytrs", func_str) || occursin("hetrs", func_str)
+        return (:unexpected_error,
+                "Unexpected positive return code from solve routine",
+                "Solve routine $func returned info=$info which should not happen.")
+
+    # General eigenvalue problem
+    elseif occursin("ggev", func_str) || occursin("gges", func_str)
+        if info <= size
+            return (:convergence_failure,
+                    "QZ iteration failed",
+                    "The QZ iteration failed to compute all eigenvalues. Elements 1:$(info-1) converged.")
+        else
+            return (:unexpected_error,
+                    "Unexpected error in generalized eigenvalue problem",
+                    "Info value $info is unexpected for $func.")
+        end
+
+    # LDLT factorization
+    elseif occursin("ldlt", func_str)
+        return (:singular_matrix,
+                "Matrix is singular",
+                "The $(info)-th pivot is zero. The factorization has been completed but division will produce infinity.")
+
+    # Default case
+    else
+        return (:unknown_error,
+                "Unknown positive return code",
+                "Function $func returned info=$info. Consult LAPACK documentation for details.")
+    end
+end
+
+"""
+    log_blas_info(func::Symbol, info::Integer, verbose::LinearVerbosity; 
+                  extra_context::Dict{Symbol,Any} = Dict())
+
+Log BLAS/LAPACK return code information with appropriate verbosity level.
+"""
+function log_blas_info(func::Symbol, info::Integer, verbose::LinearVerbosity;
+                       extra_context::Dict{Symbol,Any} = Dict())
+    category, message, details = interpret_blas_code(func, info)
+
+    # Determine appropriate verbosity field based on category
+    verbosity_field = if category in [:singular_matrix, :not_positive_definite, :convergence_failure]
+        verbose.numerical.blas_errors
+    elseif category == :invalid_argument
+        verbose.error_control.blas_invalid_args
+    else
+        verbose.numerical.blas_info
+    end
+
+    # Build structured message components
+    msg_main = "BLAS/LAPACK $func: $message"
+    msg_details = !isempty(details) ? details : nothing
+    msg_info = info
+
+    # Build complete message with all details
+    full_msg = if !isempty(extra_context) || msg_details !== nothing
+        parts = String[msg_main]
+        if msg_details !== nothing
+            push!(parts, "Details: $msg_details")
+        end
+        push!(parts, "Return code (info): $msg_info")
+        if !isempty(extra_context)
+            for (key, value) in extra_context
+                push!(parts, "$key: $value")
+            end
+        end
+        join(parts, "\n  ")
+    else
+        "$msg_main (info=$msg_info)"
+    end
+
+    # Use proper @SciMLMessage syntax
+    @SciMLMessage(full_msg, verbosity_field, :blas_return_code, :numerical)
+end
+
+"""
+    check_and_log_lapack_result(func::Symbol, result, verbose::LinearVerbosity;
+                                extra_context::Dict{Symbol,Any} = Dict())
+
+Check LAPACK operation result and log appropriately based on verbosity settings.
+Returns true if operation was successful, false otherwise.
+"""
+function check_and_log_lapack_result(func::Symbol, result, verbose::LinearVerbosity;
+                                     extra_context::Dict{Symbol,Any} = Dict())
+    # Extract info code from result
+    info = if isa(result, Tuple) && length(result) >= 3
+        # Standard LAPACK return format: (A, ipiv, info, ...)
+        result[3]
+    elseif isa(result, LinearAlgebra.Factorization) && hasfield(typeof(result), :info)
+        result.info
+    else
+        0  # Assume success if we can't find info
+    end
+
+    if info != 0
+        log_blas_info(func, info, verbose; extra_context=extra_context)
+    elseif verbose.numerical.blas_success !== Verbosity.None()
+        success_msg = "BLAS/LAPACK $func completed successfully"
+        @SciMLMessage(success_msg, verbose.numerical.blas_success, :blas_success, :numerical)
+    end
+
+    return info == 0
+end
+
+# Extended information for specific BLAS operations
+"""
+    get_blas_operation_info(func::Symbol, A, b, verbose::LinearVerbosity)
+
+Get additional information about a BLAS operation for enhanced logging.
+Condition number is computed based on the condition_number verbosity setting.
+"""
+function get_blas_operation_info(func::Symbol, A, b, verbose::LinearVerbosity)
+    info = Dict{Symbol,Any}()
+
+    # Matrix properties
+    info[:matrix_size] = size(A)
+    info[:matrix_type] = typeof(A)
+    info[:element_type] = eltype(A)
+
+    # Condition number (based on verbosity setting)
+    should_compute_cond = verbose.numerical.condition_number !== Verbosity.None()
+    if should_compute_cond && size(A, 1) == size(A, 2)
+        try
+            cond_num = cond(A)
+            info[:condition_number] = cond_num
+
+            # Log the condition number if enabled  
+            cond_msg = "Matrix condition number: $(round(cond_num, sigdigits=4)) for $(size(A, 1))×$(size(A, 2)) matrix in $func"
+            @SciMLMessage(cond_msg, verbose.numerical.condition_number, :condition_number, :numerical)
+        catch
+            # Skip if condition number computation fails
+        end
+    end
+
+    # RHS properties if provided
+    if b !== nothing
+        info[:rhs_size] = size(b)
+        info[:rhs_type] = typeof(b)
+    end
+
+    # Memory usage estimate
+    mem_bytes = prod(size(A)) * sizeof(eltype(A))
+    info[:memory_usage_MB] = round(mem_bytes / 1024^2, digits=2)
+
+    return info
+end
+