# Sparse matrices copy arrays, do not view them

**URL:** <https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505>\
**Category:** Performance\
**Tags:** question, sparse, matrices\
**Created:** [January 31, 2024, 9:19am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505 "2024-01-31T09:19:49Z")\
**Posts on this page:** 19\
**Page:** 1

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 9:19am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/1 "2024-01-31T09:19:49Z")

</div>

I wrote this program, which preallocates the csc arrays.  
I have to do a lot of kronecker products and i wanted to re-use this arrays and to write in them.  
Each iteration would have a different @view inside them, but if I use @view, it gives me an error.

```julia
N::Int32 = 4
m::Int32 = 2

nzval_Sz = Vector{Int8}(undef, N)
nzval_Sz[1] = 1
nzval_Sz[2] = -1

rowval_Sz = Vector{Int32}(undef, N)
rowval_Sz[1] = 1
rowval_Sz[2] = 2

colptr_Sz = Vector{Int32}(undef, N+1)
colptr_Sz[1] = 1
colptr_Sz[2] = 2
colptr_Sz[3] = 3 #last one is tot_nnz+1

sparsematrix = SparseMatrixCSC(m, m, colptr_Sz[1:m+1], rowval_Sz[1:m], nzval_Sz[1:m])

nzval_Sz[1] = 22

display(sparsematrix)

```

I tried to see if they were making an automatical view of these vectors, but after changing a value of nzval, the matrix didn’t change.

How can i make views of these vectors?  
Can I use array instead of vector, or does it lead to internal conversions?  
Why did they choose vectors?

related unanswered questions:

> [@Extend SparseMatrixCSC interface (nzrange - rowvals - nonzeros) to views of sparse matrices](https://discourse.julialang.org/t/extend-sparsematrixcsc-interface-nzrange-rowvals-nonzeros-to-views-of-sparse-matrices/19837):
>
> I want to extend the nzrange - rowvals - nonzeros interface from SparseMatrixCSC to views of sparse matrices. That would allow to re-formulate many of the standard algorithms for SparseMatrixCSC and extend them to the views without further effort. The existing exported functions SparseArrays.nzrange … are currently only defined for SparseMatrixCSC. Most existing algorithms for SparseMatrixCSC access the data fields colptr, rowval and nzval of SparseMatrixCSC directly. It would be easy to repla…

P.S: If i write my own sparsematrixcsc subtype, just like i did fo some matrices before, do I need to rewrite most of the operations, or do they already work using views etc?

If I want to save the kronecker product:

```julia
result = kronecker(sparse_Identity, sparsematrix)

```

not as a standalone sparse matrix, but instead re-using sparsematrix itself, is it possible, or is it going to dodge what i write, to instead do what he wants?

related question:

> [@Kronecker doesn't work, kron yes, in this case](https://discourse.julialang.org/t/kronecker-doesnt-work-kron-yes-in-this-case/109474/6):
>
> Also, your pre-allocation of matr\_C isn’t effective or needed, since that allocation is not used. Unlike in Fortran or Matlab, when you write matr\_C = ..., this does not copy the right-hand side into matr\_C. Instead, the label matr\_C in effect is “re-pointed” to whatever is on the right-hand side, in this case a KroneckerProduct type as defined in the Kronecker package. So the memory you allocated for matr\_C is discarded and eventually cleaned up by Julia’s garbage collector. This is an impor…

---

<div class="post-metadata">

**Author:** ![nilshg](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nilshg/32/2283_2.png) [@nilshg](https://discourse.julialang.org/u/nilshg)\
**Post date:** [January 31, 2024, 9:41am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/2 "2024-01-31T09:41:08Z")

</div>

I don’t actually understand your question, but one thing that might help clarify, you say:

> [@Enlil50](#):
>
> Can I use array instead of vector

but:

```julia
julia> Vector{Int32}
Vector{Int32} (alias for Array{Int32, 1})

```

`Vector{T}` is just `Array{T, 1}`, i.e. an `Array` with the dimension parameter set to 1. So saying “can I use an array instead of a vector” doesn’t really make sense, unless you are asking for an array with specific dimensions (e.g. “can I use an `Array{Int32, 2}`” - which is a `Matrix` - or `Array{Int32, 7}`)

---

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 9:42am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/3 "2024-01-31T09:42:10Z")

</div>

I heard somewhere that it would lead to a conversion, thus my confusion regarding the topic (I asked to stack overflow before asking here, cause i didn’t want to flod this place with questions, but it seems it’s more accurate to just ask here to not incur in false informations)

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [January 31, 2024, 9:51am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/4 "2024-01-31T09:51:50Z")

</div>

> [@Enlil50](#):
>
> If I want to save the kronecker product:
> 
> ```julia
> result = kronecker(sparse_Identity, sparsematrix)
> 
> ```
> 
> not as a standalone sparse matrix, but instead re-using sparsematrix itself, is it possible, or is it going to dodge what i write, to instead do what he wants?

You cannot really reuse the memory of a sparse matrix - especially not when the result has much more entries.

> [@Enlil50](#):
>
> I tried to see if they were making an automatical view of these vectors, but after changing a value of nzval, the matrix didn’t change.

You can change the internals of a `SparseMatrixCSC` directly if you wish - however that might be dangerous if you don’t truly know what you are doing. So I’d recommend not doing that. It looks like:

```julia
sparsematrix.nzval[1] = 22

```

If you know the coordinates of the value you want to change, just do it directly:

```julia
sparsematrix[1,1] = 22

```

This does not allocate a new matrix (or anything at all) if there was a value at that place before hand.

By the way:

```julia
sparsematrix = SparseMatrixCSC(m, m, colptr_Sz[1:m+1], rowval_Sz[1:m], nzval_Sz[1:m])

```

Here you explicitely **copy** the matrices `colptr_Sz`, `rowval_Sz` and `nzval_Sz`. Slicing in Julia copies! You could try again with `@views` in front. This macro converts every slicing operation to a view instead:

```julia
@views sparsematrix = SparseMatrixCSC(m, m, colptr_Sz[1:m+1], rowval_Sz[1:m], nzval_Sz[1:m])

```

---

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 9:55am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/5 "2024-01-31T09:55:06Z")

</div>

I esplicitely changed an element of those arrays, to see if sparsematrixcsc was copying them or not.

I tried with @views, but it still syas:

```julia
ERROR: LoadError: MethodError: no method matching SparseMatrixCSC(::Int32, ::Int32, ::SubArray{Int32, 1, Vector{Int32}, Tuple{UnitRange{Int64}}, true}, ::SubArray{Int32, 1, Vector{Int32}, Tuple{UnitRange{Int64}}, true}, ::SubArray{Int8, 1, Vector{Int8}, Tuple{UnitRange{Int64}}, true})

Closest candidates are:
  SparseMatrixCSC(::Integer, ::Integer, ::Vector, ::Vector, ::Vector)
   @ SparseArrays ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/SparseArrays/src/sparsematrix.jl:34
  SparseMatrixCSC(::Any, ::Any, ::SparseArrays.ReadOnly, ::SparseArrays.ReadOnly, ::Vector)
   @ SparseArrays ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/SparseArrays/src/sparsematrix.jl:47
  SparseMatrixCSC(::UniformScaling, ::Integer, ::Integer)
   @ SparseArrays ~/.julia/juliaup/julia-1.10.0+0.x64.linux.gnu/share/julia/stdlib/v1.10/SparseArrays/src/sparsematrix.jl:2190

```

About

> You cannot really reuse the memory of a sparse matrix - especially not when the result has much more entries.

I thought of preallocating the array, so that it has still space to accomodate the result. This would lead in less allocations.

Note that I’m not going to append values at the end of the slicing, still inside the preallocated array. Instead I need to rewrite all the values of the array (if we have just 1 kronecker product and a preallocated array of 4 elements. If we have 1 kronceker product and the preallocated array was of 8 elements, then i don’t need to rewrite all its elements, but just half of them).

I just wanted some stuff to do it in place, but if it needs a temporary array, what I’m doing (avoiding excessive allocations) could be meaningless

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [January 31, 2024, 10:42am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/6 "2024-01-31T10:42:36Z")

</div>

> [@Enlil50](#):
>
> I tried with @views, but it still syas:

Ok that is an unfortunate error because the constructor expects proper arrays. You could use `resize!` to shrink the arrays instead. This does not change the underlying memory, so when the code uses `resize!` again to make room for new elements, no actual memory needs to be allocated. More Julia would be to use `sizehint!` ([manual entry](https://docs.julialang.org/en/v1/base/collections/#Base.sizehint!)) instead.

```julia
N::Int32 = 4
m::Int32 = 2
nzval_Sz = [1,-1]
rowval_Sz = [1,2]
colptr_Sz = [1,2,3]#last one is tot_nnz+1

# this could cause allocations because space is preallocated
Base.sizehint!(nzval_Sz, N)
Base.sizehint!(rowval_Sz, N)
Base.sizehint!(colptr_Sz, N+1)

sparsematrix = SparseMatrixCSC(m, m, colptr_Sz, rowval_Sz, nzval_Sz)

nzval_Sz[1] = 22 # now changes the contents of sparsematrix

```

> [@Enlil50](#):
>
> I just wanted some stuff to do it in place, but if it needs a temporary array, what I’m doing (avoiding excessive allocations) could be meaningless

Most operations with `SparseArrayCSC` are not in-place. This whole preallocation buisness done above only works for operations where `SparseArrayCSC` uses `resize!` (or `push!` or similar) to increase the size of the vectors, which is the case for setting the value at an index but not much else I think.

I still think you are worrying too much about performance details before you have an actual functioning code. I think a good approach to (numerical) computing is

1. Make it work (i.e. have something that produces some output)
2. Make it right (i.e. the output is correct)
3. Make it fast

in that order.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [January 31, 2024, 10:52am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/7 "2024-01-31T10:52:38Z")

</div>

> [@Enlil50](#):
>
> ```julia
> Closest candidates are:
> SparseMatrixCSC(::Integer, ::Integer, ::Vector, ::Vector, ::Vector)
> 
> ```

This tells you all of the problem: Sparse matrices expect `Vector`. This is not a case of overly restrictive dispatch – this is really expected in the sparse matrix type.

Now julia `Vector` is in some sense a glorified pointer. So you can do the following:

```julia
julia> sparsematrix = SparseMatrixCSC(m, m, unsafe_wrap(Array, pointer(colptr_Sz), (m+1,)) , unsafe_wrap(Array, pointer(rowval_Sz), (m,)) , unsafe_wrap(Array, pointer(nzval_Sz), (m,)))
2×2 SparseMatrixCSC{Int8, Int32} with 2 stored entries:
 1 ⋅
 ⋅ -1

julia> nzval_Sz[1] = 22
22

julia> display(sparsematrix)
2×2 SparseMatrixCSC{Int8, Int32} with 2 stored entries:
 22 ⋅
  ⋅ -1

```

Use of `unsafe_wrap` implies that you are now responsible for memory management. This especially means that you must make sure to keep `colptr_Sz`, `rowval_Sz` and `nzval_Sz` alive until the last access to the wrapped array (on pain of use-after-free bugs). Furthermore, you cannot safely resize these arrays anymore until the last access to the wrapped array (that would `realloc` their memory, leading to the same kind of use-after-free bug).

Otherwise, this should behave like a view for most purposes (including the small allocation for the wrapper object).

PS. It might be possible to punt on the memory management by installing a finalizer on the wrapped array that keeps the parent alive. Does anyone here know whether that (1) works in practice and (2) works in spec (i.e. is officially supported)? The docs for `finalizer` already suggest that finalizers on arrays must throw, but that is probably a bug in the docs and not the implementation?

---

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 11:18am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/8 "2024-01-31T11:18:40Z")

</div>

Is it then possible to overwrite the contents of these arrays, which are the storage of SparseMatrix (and in our case only m elements of these N-element-arrays are used) with the result of kronecker (which creates an N-element-array) without temporary arrays?

If it’s not possible, then, even if this is kind of useful because we are not duplicating the inputs, I don’t think it gives me a real advantage as I need to reach really big matrices after multiple kronecker products (and consequent multiple internal allocations)

---

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 11:33am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/9 "2024-01-31T11:33:07Z")

</div>

I already have an algorithm which works, so I’m just testing the different options, before creating something that is hard to rewrite from the beginning.

However I agree with those points. It’s just an easy NRG algorithm the one I’m implementing. I would then pass to tensor networks, so before doing it I would like to consolidate my Julia basis

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [January 31, 2024, 11:39am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/10 "2024-01-31T11:39:05Z")

</div>

You can look at what `kron` does when called with `SparseMatrixCSC`:

```julia-repl
julia> using LinearAlgebra, SparseArrays;
julia> sparsematrix = sprand(5,5,0.1); 
julia> @less kron(sparsematrix, sparsematrix) # opens text viewer at code that would run

```

There we can see this code:

```julia
function kron(A::AbstractSparseMatrixCSC, B::AbstractSparseMatrixCSC)
    mA, nA = size(A)
    mB, nB = size(B)
    mC, nC = mA*mB, nA*nB
    Tv = typeof(oneunit(eltype(A))*oneunit(eltype(B)))
    Ti = promote_type(indtype(A), indtype(B))
    C = spzeros(Tv, Ti, mC, nC)
    sizehint!(C, nnz(A)*nnz(B))
    return @inbounds kron!(C, A, B)
end

```

This:

1. computes the size of the new matrix,
2. creates a new `SparseMatrixCSC` and uses `sizehint!` to allocate the appropriate size and then
3. uses `kron!` to perform the kronecker product into the newly allocated array.

This is interesting because you can of course allocate the output array yourself and call `kron!` to overwrite it. Note that this exclamation mark `!` at the of function names indicates that this functions mutates (some of) its inputs and you need to take special care when using them! Most importantly, in this case you cannot reuse one of the input arrays as output array like `kron!(A,B,A)`. So you need to allocate a 3rd array.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [January 31, 2024, 11:43am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/11 "2024-01-31T11:43:05Z")

</div>

If you already have a working code, you should [profile](https://docs.julialang.org/en/v1/manual/profile/) it to see which part is slowest and optimize that. There is no point in optimizing a part of the program that makes up only 1% of the runtime (because you can at most gain 1% of runtime).

If you need help with profiling, or interpreting the results or optimizing the relevant parts of the code, you are welcome to post those things here and usually a couple of people will give you helpful comments on how to optimize the hell out of your code snippet 😃

---

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 11:48am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/12 "2024-01-31T11:48:02Z")

</div>

If I know I have to do x kronecker products, I could preallocate a single array with the sole purpose of storing all the temporary arrays. Thus allocating 1 time the biggest one (and slicing/viewing it for the various iterations), instead of allocating and garbage collecting x temporary arrays.

Is it faster? How can I say to him to use the preallocated one (viewing a slice [1:m] of it)?

( I need to see how kronecker of Kronecker.jl works and if I could do it the exact same way I would do it for kron)

Example:

I need to define:

```
Sparse_matrix

```

And:

```
A = kronecker(Identity, sparse_matrix) 

```

In the actual algorithm there are even other operations involved

Now A has double the dimensions of sparse matrix. In the next iteration I need to do:

```
B = kronecker(Identity, A) 

```

That means that I’m no more using sparsematrix.  
That also means that if I allocate 2 big arrays and view into them, I could not spend time in new allocation and garbage collections

> [@abraemer](#):
>
> If you need help with profiling, or interpreting the results or optimizing the relevant parts of the code, you are welcome

Thanks, I will do it later in this month. I still need to understand how to profile well. Up until now I’m just judging the complexity of the various algorithms inside the program and doing it manually, but I know I need to switch to a proper method. I just don’t have time right now.

---

<div class="post-metadata">

**Author:** ![abraemer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/abraemer/32/51403_2.png) [@abraemer](https://discourse.julialang.org/u/abraemer)\
**Post date:** [January 31, 2024, 11:56am UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/13 "2024-01-31T11:56:01Z")

</div>

> [@Enlil50](#):
>
> Is it faster?

This only a measurement can tell you. Thus I keep recommending to actually implement something before wasting hours on improving something that might not need improvement (and the risk your “improvement” actually turns out worse!).

> [@Enlil50](#):
>
> How can I say to him to use the preallocated one (viewing a slice [1:m] of it)?

You don’t do that directly. You just create an empty `SparseMatrixCSC` and use `sizehint!` to tell it the maximal size it will need. Then you keep using that and it won’t allocate until you exceed that size.

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [January 31, 2024, 12:06pm UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/14 "2024-01-31T12:06:48Z")

</div>

If your matrices are non-tiny, I think you should think really hard about whether materializing kronecker products is a productive way forward.

The amount of data stored in `C = kron(A,B)` is the product of the amounts of data stored in A and B; but in reality you only need the sum. So you have just pissed away an entire polynomial order of complexity. It would be simple to write a wrapper matrix type that lazily evaluates the kronecker – and that would probably be faster for read-access (memory transfers more important than compute).

So this depends very much on what you plan to do with the kron products!

So whatever you want to do with `C`, think hard whether that thing can work directly on `A` and `B`. Maybe ask around here. Maybe take another look at why you thought kronecker product was the right way forward – people who prove theorems often write something convenient, with the implicit understanding that nobody would ever implement that in this stupid way (e.g. Leibniz formula for determinants with a sum over all permutations, lol)…

Apologies if you already conclusively know that you must materialize the `kron`.

---

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 1:20pm UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/15 "2024-01-31T13:20:20Z")

</div>

There are 2 ways:

Kronecker for matrices  
Outer product for tensors.

I’m now using the first way.

I will provide the naïve program in few hours

---

<div class="post-metadata">

**Author:** ![foobar\_lv2](https://avatars.discourse-cdn.com/v4/letter/f/ee59a6/32.png) [@foobar\_lv2](https://discourse.julialang.org/u/foobar_lv2)\
**Post date:** [January 31, 2024, 1:34pm UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/16 "2024-01-31T13:34:30Z")

</div>

> [@Enlil50](#):
>
> Outer product for tensors.

Outer products have the same issue.

Take the example where you need to work with projections onto a subspace. That is, the linear map `x -> x - v * (v' * x)` for some unit vector `v`.

This is a fine linear operator, and in most cases you would be stupid to represent that as a matrix `LinearAlgebra.Diagonal(fill(1.0, length(v))) - v * v'`. Not just because of construction speeds, but for everything.

(it makes sense for small known dimensions using StaticArrays, though – and there surely are other settings where you want to represent that as a matrix. Just think very hard whether that is good!)

---

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 4:32pm UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/17 "2024-01-31T16:32:50Z")

</div>

I’ve wrote this simplified piece of code, but I’m stuck, as julia has no pointers like C,C++

```julia
using SparseArrays
using LinearAlgebra
using Kronecker

#Definition of starting parameters
Delta::Float64 = 0
max_states::Int16 = 10 
NIter::Int32 = 2 #Final lattice size: 2*NIter + 2
init_rowsize::Int32 = 2

# Sz = Array{Float64}(undef, 2, 2)
# Sp = Array{Float64}(undef, 2, 2)
# Sm = Array{Float64}(undef, 2, 2)

Sz = [[0.5, 0] [0, -0.5]]
Sp = [[0, 0] [1.0, 0]]
Sm = [[0, 1.0] [0, 0]]
Ham = zeros(2,2)
# If I first preallocate them, and then check their type, the last statement overwrites the preallocation. 
#I would like to use this syntax and to control their type. Didn't found anything. Same for Ham.

#In C I would have wrote these Block_stuff as pointers ,see in the loop
Block_Sp = Sp
Block_Sm = Sm
Block_Sz = Sz
Block_Ham = Ham

for i = 1:max_states
    actual_rowsize = 2^i

    #Hamiltonian constructor
    H_super = sparse( kroneckersum(Block_Ham, Block_Ham) .- Delta*kronecker(Block_Sz, Block_Sz) .+ 0.5*(kronecker(SBlock_Spp, Block_Sm) .+ kronecker(Block_Sm, Block_Sp)) )
    H_dense_super = collect(H_super) #I make a dense matrix as this is a dummy example with full diagonalization.

    #Higher dimension kron
    Block_Ham = H_super
    Block_Sz = sparse(kronecker(sparse(I, actual_rowsize, actual_rowsize), Block_Sz))
    Block_Sp = sparse(kronecker(sparse(I, actual_rowsize, actual_rowsize), Block_Sp))
    Block_Sm = sparse(kronecker(sparse(I, actual_rowsize, actual_rowsize), Block_Sm))

end

```

Now I’ve made Block\_“stuff” ambigous, cause I globally associated it with “stuff”  
What I wanted to do is to create a pointer to this stuff, create a higher dimensional “stuff” with kron, and change the address of the pointers to the new “stuff”.  
This way I’m able to use this pointer inside the construction of H\_super iteratively

I don’t know if there is a good way or if the clearest way is to preallocate 2 big sparse arrays for each “stuff” and go with them, without removing them from memory

P.S: the initial Ham is a 0 matrix, but doing the kronecker product out of 0 matrix leds to an error.

---

<div class="post-metadata">

**Author:** ![pitsianis](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/pitsianis/32/26588_2.png) [@pitsianis](https://discourse.julialang.org/u/pitsianis)\
**Post date:** [January 31, 2024, 6:20pm UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/18 "2024-01-31T18:20:21Z")

</div>

You should think of the mathematical notion of the Kronecker product as a means to compose a linear operator without forming it.

Remember, `kronecker(A,B) * x` does not do `kron(A,B)*x` but `reshape(B * reshape(x,n,n) * A', n*n,1)` which is much cheaper in time and space. Make sure you understand why.

I see two distinct approaches here:

1. Build a sparse matrix. The `kron` of sparse matrices is sparse, so if you make the initial arrays `Sz`, etc, sparse, then the result will be sparse. But sparse forces the formation of the Kronecker product, so do not do `sparse(kronecker())` use simply `kron`

2. Build a hierarchy of `kronecker` operators. Because `kronecker` does not form the Kronecker product, you do get what you want with the pointers to matrices that cannot be done with approach 1. See

```julia
julia> A = ones(2,2)
2×2 Matrix{Float64}:
 1.0 1.0
 1.0 1.0

julia> M = kronecker(A,A)
4×4 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}:
 1.0 1.0 1.0 1.0
 1.0 1.0 1.0 1.0
 1.0 1.0 1.0 1.0
 1.0 1.0 1.0 1.0

julia> A .= 2  
2×2 Matrix{Float64}:
 2.0 2.0
 2.0 2.0

julia> M
4×4 Kronecker.KroneckerProduct{Float64, Matrix{Float64}, Matrix{Float64}}:
 4.0 4.0 4.0 4.0
 4.0 4.0 4.0 4.0
 4.0 4.0 4.0 4.0
 4.0 4.0 4.0 4.0

```

Finally, initialize your matrices as sparse, so the result will be sparse. Do not use `sparse` within the for loop.

---

<div class="post-metadata">

**Author:** ![Enlil50](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/enlil50/32/206539_2.png) [@Enlil50](https://discourse.julialang.org/u/Enlil50)\
**Post date:** [January 31, 2024, 9:06pm UTC](https://discourse.julialang.org/t/sparse-matrices-copy-arrays-do-not-view-them/109505/19 "2024-01-31T21:06:01Z")

</div>

The main issue is that, written A, I need to create a new\_A using A.  
In addition new\_A is of different dimension than A.

I didn’t really understand your point 2, even if I understand that kronecker products produce an expression: not allocated stuff. However I’m writing a simple example. In general I would like to diagonalize it, so I need it fully allocated.

I would have created something like this, but i didn’t run it. just to sketch the idea. Don’t know if it’s functional.

```julia
#Definition of constant
NIter::Int32 = 3
init_N::Int32 = 2 #inizial size of matrices
N::Int32 = init_N ^ NIter #final size of matrices

#sparsization of Sz
nzval_Sz = Vector{Int8}(undef, N)
nzval_Sz[1] = 0.5
nzval_Sz[2] = -0.5
rowval_Sz = Vector{Int32}(undef, N)
rowval_Sz[1] = 1
rowval_Sz[2] = 2
colptr_Sz = Vector{Int32}(undef, N+1)
colptr_Sz[1] = 1
colptr_Sz[2] = 2
colptr_Sz[3] = 3 #last one is tot_nnz+1

#definition of temporary arrays
nzval_Sz_temp = Vector{Int8}(undef, N)
rowval_Sz_temp = Vector{Int32}(undef, N)
colptr_Sz_temp = Vector{Int32}(undef, N+1)

for i = 1:NIter
    m = 2^i
    m_next = 2^(i+1)

    if i%2 != 0
        Sz = SparseMatrixCSC(m, m, unsafe_wrap(Array, pointer(colptr_Sz), (m+1,)) , unsafe_wrap(Array, pointer(rowval_Sz), (m,)) , unsafe_wrap(Array, pointer(nzval_Sz), (m,)))
        #H_super(Sz)  
        Sz_temp = SparseMatrixCSC(m_next, m_next, unsafe_wrap(Array, pointer(colptr_Sz), (m_next+1,)) , unsafe_wrap(Array, pointer(rowval_Sz), (m_next,)) , unsafe_wrap(Array, pointer(nzval_Sz), (m_next,)))    
        Sz_temp .= kronecker(sparse(I, m, m), Sz) #I'm writing on Sz_temp arrays #that means that next Sz I need to use is Sz_temp
    else
        Sz_temp = SparseMatrixCSC(m, m, unsafe_wrap(Array, pointer(colptr_Sz), (m+1,)) , unsafe_wrap(Array, pointer(rowval_Sz), (m,)) , unsafe_wrap(Array, pointer(nzval_Sz), (m,)))    
        #H_super(Sz_temp)
        Sz = SparseMatrixCSC(m_next, m_next, unsafe_wrap(Array, pointer(colptr_Sz), (m_next+1,)) , unsafe_wrap(Array, pointer(rowval_Sz), (m_next,)) , unsafe_wrap(Array, pointer(nzval_Sz), (m_next,)))
        Sz .= kronecker(sparse(I, m, m), Sz_temp) #I'm writing on Sz arrays #that means that next Sz I need to use is Sz
    end

end

```

This way I alternated the use of the temporary arrays, given the iteration index being odd or even.

A clearer version is the following, but it requires to somehow delete some data (I still don’t know how it’s done on julia and can’t find anything regarding manual deallocation of stuff)

```julia
#definition of constants
NIter::Int32 = 3
init_N::Int32 = 2

#sparsization of Sz
nzval_Sz = Vector{Int8}(undef, init_N)
nzval_Sz[1] = 0.5
nzval_Sz[2] = -0.5
rowval_Sz = Vector{Int32}(undef, init_N)
rowval_Sz[1] = 1
rowval_Sz[2] = 2
colptr_Sz = Vector{Int32}(undef, init_N+1)
colptr_Sz[1] = 1
colptr_Sz[2] = 2
colptr_Sz[3] = 3

Sz = SparseMatrixCSC(init_N, init_N, colptr_Sz, rowval_Sz, nzval_Sz) #copy of Sz

for i = 1:NIter
    m = 2^i
    Block_Sz = Sz #copy of the data
    #somehow delete Sz
    #H_super(Block_Sz)
    Sz = kronecker(sparse(I, m, m), Block_Sz) #definition of new Sz
    #somehow delete Block_Sz
end

```

however it involves a lot of allocations and deallocation if NIter is big.  
The data involved is of the same order, as there are always 2 sparse matrices at the same time.