What range of LinearAlgebra does Enzyme support?

Hi all,

I would like to use the LinearAlgebra package in Enzyme.
Since Julia’s sparse matrix package, SparseArrays, did not seem to be available, I created my own sparse-matrix library.

When I looked at the Enzyme testing program, I found a BLAS test.
Therefore, using that library, I created a program for the ICCG method to solve the simultaneous equations.
I used BLAS.dot(), but in some cases the program crashed when used.
In the program, resl0 = sqrt(BLAS.dot(res, res)) is not a problem. But sk = BLAS.dot(res, zk2) and pkapk = BLAS.dot(pk, zk) crash the program.

I would like to know if Enzyme currently supports up to the LinearAlgebra part.
I am not even able to use basic functions such as det(), although I have not used them in this ICCG program.
If you know more about this, please let me know.

Some of the programs are as follows:
#------------------------------------------------------------------------------------

#------------------------------------------------------------------------------------
function ICCG(A::MySparseMatrixCRS, b::Vector{Float64}, x0::Vector{Float64}=zeros(Float64, length(b)), itr_max::Int64=100000, tol_abs::Float64=1.0e-13, tol_rel::Float64=1.0e-13)
#
small = 1.0e-15

# Initalize working value and arrays
n = A.nrows
pk = zeros(Float64, n)
x = zeros(Float64, n)

# Calculate initial residual vector and the norm
res = copy(b) - A * x0

# L^1-norm of residual
#resl0 = sqrt(my_dot_product(res, res))
resl0 = sqrt(BLAS.dot(res, res))
tol = tol_abs + tol_rel * resl0

# Initially converged solution
if resl0 < tol
    println("Iteration: 0, Final residual = $(resl0)")
    return x
end

# create diag
diag = Vector{Int64}(undef, n)
Threads.@threads for i in 1 : n
     for k in A.row_ptr[i] : A.row_ptr[i+1] - 1
         if A.col_ind[k] == i
             diag[i] = k
        end
    end
end

# Calculate elements of diagonal preconditioning matrix
d = Vector{Float64}(undef, n)
for i in 1 : n
     d[i] = A.values[diag[i]]
     for k in A.row_ptr[i] : diag[i] - 1
         d[i] -= A.values[k]^2.0 * d[A.col_ind[k]]
    end
     d[i] = 1.0 / (d[i] + small)
end

# Start iterations
s0 = 1.0e+20
itr_used = 0
 for l in 1 : itr_max
    # forward substitution
    zk1 = forward_substitution(A, d, diag, res)

    # Backward substitution
    zk2 = backward_substitution(A, d, diag, zk1)

    # Inner product
    sk = my_dot_product(res, zk2)
    #sk = BLAS.dot(res, zk2)

    # Calculate beta
    bet = sk / s0

    # Calculate new search vector pk
    pk .= zk2 .+ bet .* pk

    # Calculate scalar product (pk.a pk) and alpha (overwrite zk)
    zk = A * pk

    # Inner product
    pkapk = my_dot_product(pk, zk)
    #pkapk = BLAS.dot(pk, zk)

    # Calculate alpha
    alf = sk / pkapk

    # Update solution vector
    x .+= alf .* pk

    # Update residual vector
    res .-= alf .* zk

    # L1-norm of residual
    #resl = my_dot_product(res, res)
    resl = BLAS.dot(res, res)

    # update
    s0 = sk
    itr_used = itr_used + 1

    # Check convergence
    rsm = sqrt(resl)
    if rsm < tol
        #println("Convergence achieved after $(itr_used) iterations, residual = $(rsm).")
        return x
    end
end

#
println("Max iteration. Not converged")
return x

end

gemm (matmul), gemv(matvec), dot (dot product), axpy (vector combine), and others exist with fast blas libraries for float32/64. Fallback derivatives for other BLAS libraries is also included. Lapack support (with the exception of some additions and current work underway, such as cholesky, etc) is presently limited. I believe det is a lapack call in julia. It shouldn’t be hard to add if you’re interested.

As for whatever else you saw, it would be hard to help without seeing the error message / backtrace.

1 Like

For example, the following changes were made in the ICCC program:
sk = my_dot_product(res, zk2) → sk = BLAS.dot(res, zk2)

my_dot_product() is a function of the linear algebra library I created in Julia.
The types of “res” and “zk2” are both Vector{Float64}.

When run with the normal REPL, Julia itself crashes and stops before checking for errors.
Therefore, we ran it in debug mode.
Sorry it is very long. Here are the errors.


[51880] signal (11.2): Segmentation fault: 11
in expression starting at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/debugger/run_debugger.jl:13
_ZN13GradientUtils15cacheForReverseERN4llvm9IRBuilderINS0_14ConstantFolderENS0_24IRBuilderDefaultInserterEEEPNS0_5ValueEib at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator10handle_dotE8BlasInfoRN4llvm8CallInstEPNS1_8FunctionERKNSt3__16vectorIbNS6_9allocatorIbEEEEPNS1_4TypeE at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator10handleBLASERN4llvm8CallInstEPNS0_8FunctionE8BlasInfoRKNSt3__16vectorIbNS6_9allocatorIbEEEE at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
ZN16AdjointGenerator26handleKnownCallDerivativesERN4llvm8CallInstEPNS0_8FunctionENS0_9StringRefERKNSt3__16vectorIbNS6_9allocatorIbEEEEPS1 at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator13visitCallInstERN4llvm8CallInstE at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN11EnzymeLogic23CreatePrimalAndGradientE14RequestContextOK15ReverseCacheKeyR12TypeAnalysisPK15AugmentedReturnb at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator28recursivelyHandleSubfunctionERN4llvm8CallInstEPNS0_8FunctionERKNSt3__16vectorIbNS5_9allocatorIbEEEEb10DIFFE_TYPEb at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator13visitCallInstERN4llvm8CallInstE at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN11EnzymeLogic23CreatePrimalAndGradientE14RequestContextOK15ReverseCacheKeyR12TypeAnalysisPK15AugmentedReturnb at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator28recursivelyHandleSubfunctionERN4llvm8CallInstEPNS0_8FunctionERKNSt3__16vectorIbNS5_9allocatorIbEEEEb10DIFFE_TYPEb at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator13visitCallInstERN4llvm8CallInstE at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN11EnzymeLogic23CreatePrimalAndGradientE14RequestContextOK15ReverseCacheKeyR12TypeAnalysisPK15AugmentedReturnb at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator28recursivelyHandleSubfunctionERN4llvm8CallInstEPNS0_8FunctionERKNSt3__16vectorIbNS5_9allocatorIbEEEEb10DIFFE_TYPEb at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN16AdjointGenerator13visitCallInstERN4llvm8CallInstE at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
_ZN11EnzymeLogic23CreatePrimalAndGradientE14RequestContextOK15ReverseCacheKeyR12TypeAnalysisPK15AugmentedReturnb at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
EnzymeCreatePrimalAndGradient at /Users/iberico_mac/.julia/artifacts/f9ff3c847f566c574f14defcc79f62b90d0c2cc8/lib/libEnzyme-15.dylib (unknown line)
EnzymeCreatePrimalAndGradient at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/api.jl:154
unknown function (ip: 0x28c54c25b)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
enzyme! at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/compiler.jl:3113
unknown function (ip: 0x28c4fb153)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
#codegen#488 at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/compiler.jl:4825
codegen at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/compiler.jl:4340 [inlined]
_thunk at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/compiler.jl:5430
_thunk at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/compiler.jl:5430 [inlined]
cached_compilation at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/compiler.jl:5464 [inlined]
#532 at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/compiler.jl:5530
JuliaContext at /Users/iberico_mac/.julia/packages/GPUCompiler/U36Ed/src/driver.jl:47
unknown function (ip: 0x281030263)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
#s1883#531 at /Users/iberico_mac/.julia/packages/Enzyme/wR2t7/src/compiler.jl:5482 [inlined]
#s1883#531 at ./none:0
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
GeneratedFunctionStub at ./boot.jl:602
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
jl_call_staged at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/method.c:540
ijl_code_for_staged at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/method.c:593
get_staged at ./compiler/utilities.jl:123
get_staged at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/utils.jl:216 [inlined]
#prepare_framecode#40 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/construct.jl:157
prepare_framecode at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/construct.jl:138 [inlined]
#prepare_call#43 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/construct.jl:249
prepare_call at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/construct.jl:235 [inlined]
#get_call_framecode#58 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/localmethtable.jl:63
get_call_framecode at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/localmethtable.jl:10 [inlined]
#evaluate_call_recurse!#65 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:274
evaluate_call_recurse! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:249 [inlined]
eval_rhs at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:422
unknown function (ip: 0x12549e8bb)
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:495
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:633
finish! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:14
finish_and_return! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:30
#evaluate_call_recurse!#65 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:288
evaluate_call_recurse! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:249 [inlined]
eval_rhs at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:422
unknown function (ip: 0x12549e8bb)
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:578
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:633
finish! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:14
finish_and_return! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:30
#evaluate_call_recurse!#65 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:288
evaluate_call_recurse! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:249 [inlined]
eval_rhs at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:422
unknown function (ip: 0x12549e8bb)
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:578
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:633
finish! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:14
finish_and_return! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:30
#evaluate_call_recurse!#65 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:288
evaluate_call_recurse! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:249 [inlined]
eval_rhs at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:422
unknown function (ip: 0x12549e8bb)
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:578
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:633
finish! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:14
finish_and_return! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:30
#evaluate_call_recurse!#65 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:288
evaluate_call_recurse! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:249 [inlined]
eval_rhs at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:422
unknown function (ip: 0x12549e8bb)
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:573
step_expr! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/interpret.jl:633
finish! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:14
finish_and_return! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:30
unknown function (ip: 0x1254c409b)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
finish_stack! at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:60
unknown function (ip: 0x1255d4a5b)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
#debug_command#85 at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:516
debug_command at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/JuliaInterpreter/src/commands.jl:448
unknown function (ip: 0x1255cc06f)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
jl_apply at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/./julia.h:1982 [inlined]
jl_f__call_latest at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/builtins.c:812
#invokelatest#2 at ./essentials.jl:892 [inlined]
invokelatest at ./essentials.jl:889 [inlined]
our_debug_command at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/DebugAdapter/src/debugger_core.jl:67
startdebug at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/DebugAdapter/src/packagedef.jl:106
startdebugger at /Users/iberico_mac/.vscode/extensions/julialang.language-julia-1.75.2/scripts/packages/VSCodeDebugger/src/VSCodeDebugger.jl:43
unknown function (ip: 0x10c400563)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
jl_apply at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/./julia.h:1982 [inlined]
do_call at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/interpreter.c:126
eval_body at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/interpreter.c:0
jl_interpret_toplevel_thunk at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/interpreter.c:775
jl_toplevel_eval_flex at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/toplevel.c:934
jl_toplevel_eval_flex at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/toplevel.c:877
ijl_toplevel_eval at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/toplevel.c:943 [inlined]
ijl_toplevel_eval_in at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/toplevel.c:985
eval at ./boot.jl:385 [inlined]
include_string at ./loading.jl:2076
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
_include at ./loading.jl:2136
include at ./Base.jl:495
jfptr_include_46545 at /Users/iberico_mac/.julia/juliaup/julia-1.10.2+0.aarch64.apple.darwin14/lib/julia/sys.dylib (unknown line)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
exec_options at ./client.jl:318
_start at ./client.jl:552
jfptr__start_82857 at /Users/iberico_mac/.julia/juliaup/julia-1.10.2+0.aarch64.apple.darwin14/lib/julia/sys.dylib (unknown line)
_jl_invoke at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:0 [inlined]
ijl_apply_generic at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/gf.c:3076
jl_apply at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/./julia.h:1982 [inlined]
true_main at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/jlapi.c:582
jl_repl_entrypoint at /Users/julia/.julia/scratchspaces/a66863c6-20e8-4ff4-8a62-49f30b1f605e/agent-cache/default-honeycrisp-HL2F7YQ3XH.0/build/default-honeycrisp-HL2F7YQ3XH-0/julialang/julia-release-1-dot-10/src/jlapi.c:731
Allocations: 282966118 (Pool: 282828385; Big: 137733); GC: 205

I have incorporated the entire finite element method program into Enzyme’s automatic differentiation.
As a linear algebra library I need the following:
dot(), inv(), det(), transpose() (or I want to calculate B’ * C * B )

The sparse arrays sparse matrix creation and direct method linear solver are perfect.
However, it still takes time?

Please format your messages (see Please read: make it easier to help you)

Also it helps tremendously if you post minimal working examples so that we can more easier reproduce what you are seeing

Sorry about that. I will be careful.
Also, thanks for the information on how to write.

Once again, I have created a program to test independently.

Using dot() in the calculation of variable sk causes the program to crash.
On the other hand, the calculation can be done well using my_dot_product(), a function of my own creation.

resl=dot(res, res) can be computed successfully.
Why is it sometimes calculated well and sometimes not?

using LinearAlgebra
using Enzyme
import Base.*
#----------------------------------------------------------------------------
# sparse matrix (CRS format)
#----------------------------------------------------------------------------
mutable struct MySparseMatrixCRS
    values::Vector{Float64}        # non zero value
    col_ind::Vector{Int64}         # index
    row_ptr::Vector{Int64}         # pointer
    nrows::Int64
    ncols::Int64
end
#----------------------------------------------------------------------------
# A * b : for my sparse matrix type
#----------------------------------------------------------------------------
function *(A::MySparseMatrixCRS, x::Vector{Float64})
    z = Vector{Float64}(undef, A.nrows)
    for i in 1 : A.nrows
        z[i] = 0.0
        for j in A.row_ptr[i] : A.row_ptr[i+1] - 1
            z[i] += A.values[j] * x[A.col_ind[j]]
        end
    end
    return z
end
#----------------------------------------------------------------------------
# dot product by myself
#----------------------------------------------------------------------------
function my_dot_product(a::Vector{Float64}, b::Vector{Float64})
    value = 0.0
    for i in 1 : length(a)
        value += a[i] * b[i]
    end
    return value
end
#----------------------------------------------------------------------------
# dense matrix -> sparse matrix
#----------------------------------------------------------------------------
function convert_to_sparse_matrix_crs(A::Matrix{Float64})
    nrows, ncols = size(A)
    values = Vector{Float64}()
    col_ind = Vector{Int64}()
    row_ptr = Vector{Int64}()
    push!(row_ptr, 1)
     for i in 1:nrows
        for j in 1:ncols
            if A[i, j] != 0.0
                push!(values, A[i, j])
                push!(col_ind, j)
            end
        end
        push!(row_ptr, length(values) + 1)
    end
    return MySparseMatrixCRS(values, col_ind, row_ptr, nrows, ncols)
end
#------------------------------------------------------------------------------------
# linear solver
#------------------------------------------------------------------------------------
function ICCG(A::MySparseMatrixCRS, b::Vector{Float64}, itr_max::Int64=100000, tol_abs::Float64=1.0e-13, tol_rel::Float64=1.0e-13)    
    #
    small = 1.0e-15
    
    # Initalize working value and arrays
    n = A.nrows
    pk = zeros(Float64, n)
    x = zeros(Float64, n)

    # Calculate initial residual vector and the norm
    res = copy(b) - A * x0

    # L^1-norm of residual
    resl0 = sqrt(dot(res, res))
    tol = tol_abs + tol_rel * resl0

    # Initially converged solution
    if resl0 < tol
        println("Iteration: 0, Final residual = $(resl0)")
        return x
    end

    # create diag
    diag = Vector{Int64}(undef, n)
    Threads.@threads for i in 1 : n
         for k in A.row_ptr[i] : A.row_ptr[i+1] - 1
             if A.col_ind[k] == i
                 diag[i] = k
            end
        end
    end

    # Calculate elements of diagonal preconditioning matrix
    d = Vector{Float64}(undef, n)
    for i in 1 : n
         d[i] = A.values[diag[i]]
         for k in A.row_ptr[i] : diag[i] - 1
             d[i] -= A.values[k]^2.0 * d[A.col_ind[k]]
        end
         d[i] = 1.0 / (d[i] + small)
    end

    # Start iterations
    s0 = 1.0e+10
    itr_used = 0
     for l in 1 : itr_max
        # forward substitution
        zk1 = forward_substitution(A, d, diag, res)

        # Backward substitution
        zk2 = backward_substitution(A, d, diag, zk1)

        # Inner product
        #sk = my_dot_product(res, zk2)
        sk = dot(res, zk2)

        # Calculate beta
        bet = sk / s0

        # Calculate new search vector pk
        pk .= zk2 .+ bet .* pk

        # Calculate scalar product (pk.a pk) and alpha (overwrite zk)
        zk = A * pk

        # Inner product
        pkapk = my_dot_product(pk, zk)
        #pkapk = dot(pk, zk)

        # Calculate alpha
        alf = sk / pkapk

        # Update solution vector
        x .+= alf .* pk

        # Update residual vector
        res .-= alf .* zk

        # L1-norm of residual
        resl = dot(res, res)

        # update
        s0 = sk
        itr_used = itr_used + 1

        # Check convergence
        rsm = sqrt(resl)
        if rsm < tol
            #println("Convergence achieved after $(itr_used) iterations, residual = $(rsm).")
            return x
        end
    end
    #
    println("Max iteration. Not converged")
    return x
end
#------------------------------------------------------------------------------------
#
#------------------------------------------------------------------------------------
function forward_substitution(A::MySparseMatrixCRS, d::Vector{Float64}, diag::Vector{Int64}, res::Vector{Float64})
    small = 1.0e-15
    zk = similar(res)
    for i in 1 : A.nrows
        zk[i] = res[i]
        for k in A.row_ptr[i] : diag[i] - 1
            zk[i] = zk[i] - A.values[k] * zk[A.col_ind[k]]
        end
        zk[i] = zk[i] * d[i]
    end
    for i in 1 : A.nrows
         zk[i] /= (d[i] + small)
    end
    return zk
end
#------------------------------------------------------------------------------------
#
#------------------------------------------------------------------------------------
function backward_substitution(A::MySparseMatrixCRS, d::Vector{Float64}, diag::Vector{Int64}, zk1::Vector{Float64})
    zk = copy(zk1)
    for i in A.nrows : -1 : 1
        for k in diag[i] + 1 : A.row_ptr[i+1] - 1
            zk[i] -= A.values[k] * zk[A.col_ind[k]]
        end
        zk[i] *= d[i]
    end
    return zk
end
#----------------------------------------------------------------------------
# test function
#----------------------------------------------------------------------------
function iccg_test(x)

    # make sparse matrix
    A_dense = Matrix{Float64}(undef, 3, 3)
    for i in 1 : 3
        for j in 1 : 3
            A_dense[i, j] = 2.0 * i + 3.0 * j^2.0 + x[i]^2.0 + x[j]^3.0
        end
    end
    # for symmetric matrix
    for i in 1 : 3
        for j in i : 3
            A_dense[i, j] = A_dense[j, i]
        end
    end
    A_sparse = convert_to_sparse_matrix_crs(A_dense)

    # right hand side
    b = Vector{Float64}(undef, 3)
    b[1] = 1.0
    b[2] = 0.0 
    b[3] = x[2]^3.0

    # solve
    u = ICCG(A_sparse, b)

    # for scalar
    return dot(u, u)
end
#----------------------------------------------------------------------------
# main
#----------------------------------------------------------------------------
x0 = rand(Float64, 3)

val = iccg_test(x0)
display(val)

df = gradient(Reverse, iccg_test, x0)
display(df)