Help with this code (deallocation/pointers on sparse matrices)

General Usage Performance
, , ,
1

This is a follow up of a recent thread I made Sparse matrices copy arrays, do not view them
I sum up all the shared ideas, with 2 code examples. This is intended as a much clearer version and prosecution of that thread, using what we came up with.

I’m currently stuck with a simple iterative algorithm by the lack of manual deallocation of stuff.
This could be solved by using pointers, but I can’t find anything concrete regarding them.
I was lucky to find a guy (@foobar_lv2) which lead me to some pointer syntax (unsafe_wrap), but I can’t find any documentation on them and they still remain unclear to me.

The issue is that, written a sparse matrix A and used it in a function I will call H_super(), I need to create a new_A using A itself.
new_A is a matrix of different dimension than A and will be used in H_super() the next iteration.

I need to use A and new_A at the same time just when creating new_A.
Thus I split the loop into 2 sections which alternate at every loop iteration:
Odd iteration and even iteration.
This way I can call H_super(A) (iter 1), H_super(new_A) (iter 2), H_super(A) (iter 3)…

A and new_A are preallocated.
I wanted to view slices of them, but sparseArrays doesn’t have a built-in feature to do that. Thus the use of unsafe_wraps.
In additon I incur in an issue when trying to write new_A expression in its buffer.
A sparse matrix buffer consists in 3 arrays and I need to put the expression in these 3 arrays. To use the sparsematrix structure I need to call SparseMatrixCSC, but this only works for arrays that are not undef arrays, waiting to be filled.

In the following program i will call A as Sz and new_A as Sz_temp. It’s actually the Sz pauli matrix

using SparseArrays
using LinearAlgebra
using Kronecker

#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{Float64}(undef, N)
nzval_Sz[1] = 0.5
nzval_Sz[2] = -0.5
rowval_Sz = Vector{UInt32}(undef, N)
rowval_Sz[1] = 1
rowval_Sz[2] = 2
colptr_Sz = Vector{UInt32}(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{Float64}(undef, N)
rowval_Sz_temp = Vector{UInt32}(undef, N)
colptr_Sz_temp = Vector{UInt32}(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,)))    
        #the line above gives error, as some values are still not allocated. what I want to do is to save teh sparse matrix below in the Sz_temp sparse matrix created by the temporary arrays preallocated.
        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,)))
        #This shares the same issue of the other condition, but with Sz and Sz_temp switched
        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

Another, much clearer version, is the following.
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)

using SparseArrays
using LinearAlgebra
using Kronecker

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

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

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

for i = 1:NIter
    m = 2^i
    Sz_temp = Sz #copy of the data
    #somehow delete Sz which is defined "globally"
    #H_super(Sz_temp)
    Sz = kronecker(sparse(I, m, m), Sz_temp) #definition of new Sz
    #somehow delete Sz_temp
end
2

It might help to explain the big picture. What is the purpose?

3

I think your question is quite ill-posed and based on fundamental misunderstanding how things work in Julia (and many other languages that are not Fortran or Matlab).
Variable names in Julia just point to some object, they have no associated, fixed memory or anything.
As such Sz_temp = Sz #copy of the data is wrong. There is no copy performed here. Sz and Sz_temp now name the exactly same thing. Two pointers to the same memory if you wish.

Please have a look at the sections relevant for you of
https://docs.julialang.org/en/v1/manual/noteworthy-differences/

3 Likes
4

It’s a simple NRG program, where I connect the Hamiltonian of one site (divide the space in single points, a site is a point, here we have just 1D) to the Hamiltonian of the other.

The Hamiltonian is H_super and it depends on the Pauli matrix Sz. This is just one of the sparse matrices I need to treat. But the steps are the same. Thus the exponential growth of size

1 Like
5

Yes, you are right. I didn’t understand well that part, so you helped me quite a lot.

However the second example throws an error due to the fact that Sz is defined before the for loop. If I understood well (I’m still new to Julia, so I finished 2 guides on it yesterday and started them 5 days ago) that’s something related to global values. Indeed the error I was seeing from the terminal was related to global variables.

I can write a version which throws this error (I still have one).

In this program I need to start with some values in the loop iteration and change them at the end of the loop iteration. I need a way to define them for the first iteration and the only way I know is defining them before the loop. The following iteration will use “new_A” which is defined at the end of the first iteration.

Maybe this could be solved if placed in a function, instead of the “main”?

Regarding

Julia values are not copied when passed to a function. If a function modifies an array, the changes will be visible in the caller.

That’s not always true, for example for SparseMatrixCSC, as discussed in the previous thread, more geared towards this aspect.

1 Like
6

Hm?

julia> using SparseArrays

julia> x = sparse([1, 3], [1, 3], [1, 1])
3×3 SparseMatrixCSC{Int64, Int64} with 2 stored entries:
 1  ⋅  ⋅
 ⋅  ⋅  ⋅
 ⋅  ⋅  1

julia> function change_sparse(s); s[1, 1] = 5; return s; end
change_sparse (generic function with 1 method)

julia> change_sparse(x);

julia> x
3×3 SparseMatrixCSC{Int64, Int64} with 2 stored entries:
 5  ⋅  ⋅
 ⋅  ⋅  ⋅
 ⋅  ⋅  1
7 
N::Int32 = 4
m::Int32 = 2

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

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

colptr_Sz = Vector{UInt32}(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 changed one of the arrays used in SparseMatrixCSC. If you display the sparse matrix you will see no change. It will be as if nzval_Sz[1]=0.5.

That’s cause I didn’t use @view. But if you try to use @view or @views it throws an error

P.S @nilshg: it seems that if I instead write the whole array inside it, instead of slicing, it does not copy.

However the issue of the the first example is that even if we just use 2 sparse matrix buffers, i need to change their dimensions too, as they increase each iteration.

8

As I told you in the other thread: slicing performs a copy in Julia. So there is an explicit copy done before the function call to SparseMatrixCSC.

And I honestly have no idea how this relates to code above. Do you somehow want to change values of the final Hamiltonian by changing something in the initial Sz? Well that just won’t at all.

10

You are right, but I thought that @view[slice] was treated just like a “pointer to array” as what happens when passing the whole array. Thus my misconception. Really sorry.

About the code above:
I have thought of those 2 versions.
The first is described in the text section. The second is quite simple.

In the second version I incur in a warning due to the fact I’m using the same name of Sz, defined outside the loop, for something created in the loop. Or at least it’s what I understood.

In the first version I’m not incurring in this issue, but there are other issues, like the fact that I’m trying to define an undefined sparse matrix of certain dimensions, through undefined arrays. I do this to fill it later with the right hand side part of the equation which is the kronecker expression.

11

Scoping in Julia is likely a bit different from what you expect. Scope of Variables · The Julia Language

12

To sum up, extending the second example, if I write this:

using SparseArrays
using LinearAlgebra
using Kronecker

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

#sparsization of Sz
nzval_Sz = Float64[0.5, -0.5]
rowval_Sz = UInt32[1,2]
colptr_Sz = UInt32[1,2,3]

Sz = SparseMatrixCSC(init_N, init_N, colptr_Sz, rowval_Sz, nzval_Sz)

for i = 1:NIter
    m = 2^i
    if i==1
        Sz_temp = Sz
    end
    #H_super(Sz_temp)
    Sz = kronecker(sparse(I, m, m), Sz_temp) #definition of new Sz
end

It gives:

┌ Warning: Assignment to `Sz` in soft scope is ambiguous because a global variable by the same name exists: `Sz` will be treated as a new local. Disambiguate by using `local Sz` to suppress this warning or `global Sz` to assign to the existing global variable.
└ @ ~/path/file.jl:167
ERROR: LoadError: UndefVarError: `Sz` not defined

If I instead write:

using SparseArrays
using LinearAlgebra
using Kronecker

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

#sparsization of Sz
nzval_Sz = Float64[0.5, -0.5]
rowval_Sz = UInt32[1,2]
colptr_Sz = UInt32[1,2,3]

Sz = SparseMatrixCSC(init_N, init_N, colptr_Sz, rowval_Sz, nzval_Sz)

for i = 1:NIter
    m = 2^i
    if i==1
        Sz_temp1 = Sz
    end
    if i%2 == 1
        println("odd")
        #H_super(Sz_temp1)
        Sz_temp2 = kronecker(sparse(I, m, m), Sz_temp1) #definition of new Sz
    else
        println("even")
        #H_super(Sz_temp2)
        Sz_temp1 = kronecker(sparse(I, m, m), Sz_temp2) #definition of new Sz
    end
end

It gives:

ERROR: LoadError: UndefVarError: `Sz_temp2` not defined

I don’t know other ways of writing this.

13

If you need Sz_temp2 to survive across iterations, then it needs to be defined in the scope outside the loop, something like this:

function foo()
     local tmp2
     for i=1:2
       if i==1
         tmp2="a"
       else
         print(tmp2)
       end
     end
end
1 Like
14

So, I think this is the solution. this is what i wrote:

function foo()
    #definition of constants
    NIter::Int32 = 2
    init_N::Int32 = 2

    #sparsization of Sz
    nzval_Sz = Float64[0.5, -0.5]
    rowval_Sz = UInt32[1,2]
    colptr_Sz = UInt32[1,2,3]

    Sz = SparseMatrixCSC(init_N, init_N, colptr_Sz, rowval_Sz, nzval_Sz)
    local Sz_temp
    local Sz_temp1
    display(Sz)

    for i = 1:NIter
        m = 2^i

        if i==1
            Sz_temp = Sz
        end

        if i%2 == 1
            #H_super(Sz_temp)
            Sz_temp1 = kronecker(sparse(I, m, m), Sz_temp) #definition of new Sz
            display(Sz_temp1)
        else
            #H_super(Sz_temp1)
            Sz_temp = kronecker(sparse(I, m, m), Sz_temp1) #definition of new Sz
            display(Sz_temp)
        end
    end
end

#main
foo()

I merged the idea of the first example, with your info about “local” and the idea of using a simple if in the first iteration.

Still, there is a lot of garbage all over the floor. Don’t know if it can be optimized with some preallocation

15

You can just initialize Sz_temp = Sz outside the loop instead of local Sz_temp, and omit the if i==1 block.

I don’t understand why your are forcing your NIter and init_N to be Int32. It’s not like it would save any significant space or make a performance difference here.

I also don’t understand why you need two, or even any, Sz_temp variables.

16

[quote=“Jeff_Emanuel, post:15, topic:109574”]
I also don’t understand why you need two, or even any, Sz_temp variables.
[/quote]

I need to write the kronecker product on something that’s not one of the matrices that are used by the product itself

I like to give types to everything and here I have a lot of shiny types

I need to try it. It will most probably work. I have no computer now

17

There is a big issue and i don’t know why.

If you write

display(issparse(Sz_temp1))

It will say: false.

For some reason it’s converted into a dense matrix

18

That’s because kronecker does not return a sparse matrix. It returns some lazy container for the product. I think that’s been explained to you earlier in other threads.

Compute a lazy Kronecker products between two matrices A and B
Basic use · Kronecker.jl

Using kronecker is beyond my expertise, (I just stepped in to provide the pointer to variable scoping.), but it looks like you want collect! to realize the product into a matrix of your choice. Presumably it works with your preallocated sparse matrices.

19

I don’t think this code does what you think it does. You just kronecker identity matrices onto your Sz. Also the matrix sizes do not really work out. Let’s forget Kronecker.jl for a moment and try to get a simple code that does what it is supposed to do:

function foo()
    #definition of constants
    NIter::Int32 = 2

    #sparsization of Sz
    Sz = sparse([0.5 0; 0 -0.5])

    for i = 1:NIter
        m = 2^i
        Sz = kron(I(m), Sz)
    end
    return Sz
end

#main
foo()

That is the condensed (and working) version of what you wrote. Note that you could do away with the loop and just kron a single identity matrix of the correct size kron(I(2^sum(1:Niter)), Sz) which is why I believe you probably wanted to do something else.

20

Yes, I know this much. However I would expected that the lazy evaluation retains the sparsity of the matrix. I don’t know why it doesn’t do this.

kron retains it instead. What I understand of lazy evaluation is that it allocates stuff when needed, but if the result is sparse, then why should it allocate the 0 i didn’t ask him to?

If I use it in the program, i expect to occupy the space of a sparse matrx. But if that’s not true, well, it could be an issue. I don’t think the issue is in the library, but rather in my understanding of it, but still, I’m confused.

Yes, I need to do stuff in each loop iteration.

21

Can you build on the simple code I wrote an add more stuff you need to do?

The result probably will be sparse (I actually don’t know - I never used Kronecker.jl). I think it’s issparse cannot detect this. See issparse is from SparseArrays.jl and probably does not know about Kronecker.jl so it just sees that it is not a SparseMatrixCSC from SparseArrays.jl so it gives false. Still let’s just omit Kronecker.jl for now let’s try to get a working version of your code and worry about optimization later. SparseMatrixCSC will likely be good enough for your needs.

2 Likes