diff --git a/docs/src/basics/common_solver_opts.md b/docs/src/basics/common_solver_opts.md index b86139935..c87fec58a 100644 --- a/docs/src/basics/common_solver_opts.md +++ b/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) \ No newline at end of file +- 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 \ No newline at end of file diff --git a/ext/LinearSolveBLISExt.jl b/ext/LinearSolveBLISExt.jl index fc4ac6bc5..1ab6d96f5 100644 --- a/ext/LinearSolveBLISExt.jl +++ b/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 diff --git a/src/LinearSolve.jl b/src/LinearSolve.jl index e8bcc5b3e..37c5e2d43 100644 --- a/src/LinearSolve.jl +++ b/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") diff --git a/src/blas_logging.jl b/src/blas_logging.jl new file mode 100644 index 000000000..a9e4f911f --- /dev/null +++ b/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 + diff --git a/src/verbosity.jl b/src/verbosity.jl index e4b73988f..3285e8d14 100644 --- a/src/verbosity.jl +++ b/src/verbosity.jl @@ -9,14 +9,21 @@ const linear_defaults = Dict{Symbol, Verbosity.Type}( :KrylovKit_verbosity => Verbosity.Warn(), :KrylovJL_verbosity => Verbosity.None(), :HYPRE_verbosity => Verbosity.Level(1), - :pardiso_verbosity => Verbosity.None() + :pardiso_verbosity => Verbosity.None(), + :blas_errors => Verbosity.Warn(), + :blas_info => Verbosity.None(), + :blas_success => Verbosity.None(), + :blas_invalid_args => Verbosity.Error(), + :condition_number => Verbosity.None() ) mutable struct LinearErrorControlVerbosity default_lu_fallback::Verbosity.Type + blas_invalid_args::Verbosity.Type function LinearErrorControlVerbosity(; - default_lu_fallback = linear_defaults[:default_lu_fallback]) - new(default_lu_fallback) + default_lu_fallback = linear_defaults[:default_lu_fallback], + blas_invalid_args = linear_defaults[:blas_invalid_args]) + new(default_lu_fallback, blas_invalid_args) end function LinearErrorControlVerbosity(verbose::Verbosity.Type) @@ -80,6 +87,10 @@ mutable struct LinearNumericalVerbosity KrylovJL_verbosity::Verbosity.Type HYPRE_verbosity::Verbosity.Type pardiso_verbosity::Verbosity.Type + blas_errors::Verbosity.Type + blas_info::Verbosity.Type + blas_success::Verbosity.Type + condition_number::Verbosity.Type function LinearNumericalVerbosity(; using_IterativeSolvers = linear_defaults[:using_IterativeSolvers], @@ -87,9 +98,14 @@ mutable struct LinearNumericalVerbosity KrylovKit_verbosity = linear_defaults[:KrylovKit_verbosity], KrylovJL_verbosity = linear_defaults[:KrylovJL_verbosity], HYPRE_verbosity = linear_defaults[:HYPRE_verbosity], - pardiso_verbosity = linear_defaults[:pardiso_verbosity]) + pardiso_verbosity = linear_defaults[:pardiso_verbosity], + blas_errors = linear_defaults[:blas_errors], + blas_info = linear_defaults[:blas_info], + blas_success = linear_defaults[:blas_success], + condition_number = linear_defaults[:condition_number]) new(using_IterativeSolvers, IterativeSolvers_iterations, - KrylovKit_verbosity, KrylovJL_verbosity, HYPRE_verbosity, pardiso_verbosity) + KrylovKit_verbosity, KrylovJL_verbosity, HYPRE_verbosity, pardiso_verbosity, + blas_errors, blas_info, blas_success, condition_number) end function LinearNumericalVerbosity(verbose::Verbosity.Type) diff --git a/test/test_blas_logging.jl b/test/test_blas_logging.jl new file mode 100644 index 000000000..117fedca0 --- /dev/null +++ b/test/test_blas_logging.jl @@ -0,0 +1,107 @@ +using LinearSolve +using LinearAlgebra +using Test +using SciMLLogging: Verbosity + +@testset "BLAS Return Code Interpretation" begin + # Test interpretation of various BLAS return codes + @testset "Return Code Interpretation" begin + # Test successful operation + category, message, details = LinearSolve.interpret_blas_code(:dgetrf, 0) + @test category == :success + @test message == "Operation completed successfully" + + # Test invalid argument + category, message, details = LinearSolve.interpret_blas_code(:dgetrf, -3) + @test category == :invalid_argument + @test occursin("Argument 3", details) + + # Test singular matrix in LU + category, message, details = LinearSolve.interpret_blas_code(:dgetrf, 2) + @test category == :singular_matrix + @test occursin("U(2,2)", details) + + # Test not positive definite in Cholesky + category, message, details = LinearSolve.interpret_blas_code(:dpotrf, 3) + @test category == :not_positive_definite + @test occursin("minor of order 3", details) + + # Test SVD convergence failure + category, message, details = LinearSolve.interpret_blas_code(:dgesvd, 5) + @test category == :convergence_failure + @test occursin("5 off-diagonal", details) + end + + @testset "BLAS Operation Info" begin + # Test getting operation info without condition number + A = rand(10, 10) + b = rand(10) + + # Test with condition_number disabled (default) + verbose_no_cond = LinearVerbosity(condition_number = Verbosity.None()) + info = LinearSolve.get_blas_operation_info(:dgetrf, A, b, verbose_no_cond) + + @test info[:matrix_size] == (10, 10) + @test info[:element_type] == Float64 + @test !haskey(info, :condition_number) # Should not compute by default + @test info[:memory_usage_MB] >= 0 # Memory can be 0 for very small matrices + + # Test with condition number computation enabled via verbosity + verbose_with_cond = LinearVerbosity(condition_number = Verbosity.Info()) + info_with_cond = LinearSolve.get_blas_operation_info(:dgetrf, A, b, verbose_with_cond) + @test haskey(info_with_cond, :condition_number) + end + + @testset "Verbosity Integration" begin + # Test with singular matrix + A = [1.0 2.0; 2.0 4.0] # Singular matrix + b = [1.0, 2.0] + + # Test with warnings enabled + verbose = LinearVerbosity( + blas_errors = Verbosity.Warn(), + blas_info = Verbosity.None() + ) + + prob = LinearProblem(A, b) + + # This should fail due to singularity but not throw + sol = solve(prob, LUFactorization(); verbose=verbose) + @test sol.retcode == ReturnCode.Failure + + # Test with all logging enabled + verbose_all = LinearVerbosity( + blas_errors = Verbosity.Info(), + blas_info = Verbosity.Info() + ) + + # Non-singular matrix for successful operation + A_good = [1.0 0.0; 0.0 1.0] + b_good = [1.0, 1.0] + prob_good = LinearProblem(A_good, b_good) + + sol_good = solve(prob_good, LUFactorization(); verbose=verbose_all) + @test sol_good.retcode == ReturnCode.Success + @test sol_good.u ≈ b_good + end + + @testset "Error Categories" begin + # Test different error categories are properly identified + test_cases = [ + (:dgetrf, 1, :singular_matrix), + (:dpotrf, 2, :not_positive_definite), + (:dgeqrf, 3, :numerical_issue), + (:dgesdd, 4, :convergence_failure), + (:dsyev, 5, :convergence_failure), + (:dsytrf, 6, :singular_matrix), + (:dgetrs, 1, :unexpected_error), + (:unknown_func, 1, :unknown_error) + ] + + for (func, code, expected_category) in test_cases + category, _, _ = LinearSolve.interpret_blas_code(func, code) + @test category == expected_category + end + end +end +