# Using PreallocationTools with LinearSolve

**URL:** <https://discourse.julialang.org/t/using-preallocationtools-with-linearsolve/98508>\
**Category:** Modelling & Simulations\
**Created:** [May 8, 2023, 9:36pm UTC](https://discourse.julialang.org/t/using-preallocationtools-with-linearsolve/98508 "2023-05-08T21:36:16Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![clm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/clm/32/45923_2.png) [@clm](https://discourse.julialang.org/u/clm)\
**Post date:** [May 8, 2023, 9:36pm UTC](https://discourse.julialang.org/t/using-preallocationtools-with-linearsolve/98508/1 "2023-05-08T21:36:16Z")

</div>

I’m trying to solve a large 2D PDE using [this example](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/) as a template. I’m effectively trying to solve some simplified (though still nonlinear) pseudo-Boussinesq Navier Stokes equations where I evolve u and temperature and solve for w using \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0. This becomes a system of equations that I have to solve for w at each time step (with my problem geometry I can assume variations in w are small, so divergence in u should be the only thing driving w)

I went through many iterations before realizing that the minimal non-working example is actually very straightforward. This problem is trivialized for simplicity, but it highlights the line that is giving me an error.

```julia
using DifferentialEquations, SparseArrays, LinearSolve, PreallocationTools

function ddt(dy, y, p, t)
    Nz,Nx,A_w,linsolve,w,divu = p

    ### I'm not certain I understand how get_tmp works, but I know it doesn't work without it.
    divu = get_tmp(divu,y)
    w = get_tmp(w,y)
    

    ### Ultimately, I would be computing "divu" before and then doing something with "w"
    
    ### This approach works when I comment out the line that makes "A_w_" sparse.
    # w[2:end-1,2:end-1] = reshape(solve(LinearProblem(A_w,reshape(divu,Nx*Nz))).u,(Nz,Nx))

    ### I can't make this approach work regardless of sparisty
    w[2:end-1,2:end-1] = reshape(solve(LinearSolve.set_b(linsolve,reshape(divu,Nx*Nz))).u,(Nz,Nx))

    ### Set all to zero for simplicity
    dy[:,:,:] .= 0.
end

Nx_= 10
Nz_= 10

A_w_ = Matrix(1.0I, Nx_*Nz_, Nx_*Nz_)
A_w_ = sparse(A_w_)

### Always going to use the same A matrix, so LU factorize and cache.
bw = ones(Nz_*Nx_)
linsolve_ = init(LinearProblem(A_w_, bw))
sol1 = solve(linsolve_)

w_ = DiffCache(zeros(Nz_+2, Nx_+2))
divu_ = DiffCache(zeros(Nz_, Nx_))

p = Nz_,Nx_,A_w_,linsolve_,w_,divu_

y0 = zeros(Nz_+2, Nx_+2, 2)

prob = ODEProblem(ddt, y0, (0.0, 1.0), p);
sol = solve(prob, FBDF(),save_everystep=false);"

```

I get an error that resembles:

```julia
First call to automatic differentiation for the Jacobian
failed. This means that the user `f` function is not compatible
with automatic differentiation. Methods to fix this include:...

```

As I’ve indicated in my comments, I can make the first method work if my A matrix is not sparse, but I cannot get the second approach (caching the LU decomposition of A) to work at all. If I allocate divu and w (and don’t pass them in through p), the code also works fine, but is far slower.

I think my issue has something to do with a misuse of PreallocationTools.jl (likely that I’m not applying it to linsolve\_?) but I’m having a hard time understanding how to do it. It looks like LazyBufferCache might be what I need somehow, but can’t wrap my head around what it does. Any help would be great.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [May 11, 2023, 8:30am UTC](https://discourse.julialang.org/t/using-preallocationtools-with-linearsolve/98508/2 "2023-05-11T08:30:40Z")

</div>

> [@clm](#):
>
> `This approach works when I comment out the line that makes "A_w_" sparse.`

Try forcing it to `SparspakFactorization()`, i.e. `w[2:end-1,2:end-1] = reshape(solve(LinearProblem(A_w,reshape(divu,Nx*Nz)),SparspakFactorization()).u,(Nz,Nx))`. If that works I can update the automated algorithm to handle this case.

> [@clm](#):
>
> ```julia
> ### I can't make this approach work regardless of sparisty
> w[2:end-1,2:end-1] = reshape(solve(LinearSolve.set_b(linsolve,reshape(divu,Nx*Nz))).u,(Nz,Nx))
> 
> ```

Yeah for that you should probably use the `GeneralLazyBufferCache`.

```julia
lbc = GeneralLazyBufferCache(function (y)
    init(LinearProblem(A_w,similar(y,Nx*Nz), SparspakFactorization())
end)

```

and then `prob = lbc[y]` gives you a dual/non-dual prob, and then

```julia
w[2:end-1,2:end-1] = reshape(solve(LinearSolve.set_b(linsolve,reshape(divu,Nx*Nz))).u,(Nz,Nx))

```

should work. Note I didn’t try all of this code and just writing this on the fly since I’m short on time, but thought you might want the untested pointers instead of nothing.

---

<div class="post-metadata">

**Author:** ![clm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/clm/32/45923_2.png) [@clm](https://discourse.julialang.org/u/clm)\
**Post date:** [May 11, 2023, 5:46pm UTC](https://discourse.julialang.org/t/using-preallocationtools-with-linearsolve/98508/3 "2023-05-11T17:46:55Z")

</div>

I appreciate you taking the time. It’s giving me more avenues to try. Somewhat blindly permuting different parameters seems to have worked!

The first approach gave me the same error as before.

I wasn’t able to get the second approach to work, but “SparspakFactorization()” seems like it was the cause. Once I removed it and just defined lbc as follows, it worked.

```julia
lbc = GeneralLazyBufferCache(function (y)
    init(LinearProblem(A_w_,y))
end)
linsolve_ = lbc[bw]
sol1 = solve(linsolve_)

```

Thanks, Chris. Very much appreciated.

---

<div class="post-metadata">

**Author:** ![clm](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/clm/32/45923_2.png) [@clm](https://discourse.julialang.org/u/clm)\
**Post date:** [August 12, 2024, 8:54pm UTC](https://discourse.julialang.org/t/using-preallocationtools-with-linearsolve/98508/4 "2024-08-12T20:54:48Z")

</div>

While I’ve had this code working for a while, I’ve encountered an argument error. I don’t know if it’s been doing this all along and I haven’t noticed, or if this started when I updated a couple of packages recently. Admittedly, it doesn’t seem to be causing any trouble, but I wanted to see if it was something I was doing wrong (in case it starts causing trouble later). Here’s a very simple MWE:

```julia
def_mat(M) = [diagm(1 => ones(M-1))-I [zeros(M-1,1); 1.0]]

n = 4
A = sparse(def_mat(n))

x_vec = ones(n)

lbc_ = GeneralLazyBufferCache(function (y)
    init(LinearProblem(A,y))
end)

linsolve_w_ = lbc_[x_vec]
sol1 = solve(linsolve_w_)
linsolve_w_.b = [4.0; 3.0; 2.0; 1.0]
sol2 = solve(linsolve_w_)

```

Which works as I would hope. However, if I comment out the last 3 lines (or otherwise end a line with “linsolve\_w\_” in any way), I get the following:

```julia
ArgumentError: pointer to the SparseArrays.LibSuiteSparse.cholmod_factor_struct object is null. This can happen if the object has been serialized.

```

Should I be structuring my GeneralLazyBufferCache function differently, or is this something within SparseArrays?