# 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.

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)
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
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]
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]
_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 ./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?

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)
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)