# How to implement hybrid parallelized programs
**URL:** https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604
**Category:** Performance
**Tags:** parallel, multithreading, distributed
**Created:** [February 2, 2024, 4:34am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604 "2024-02-02T04:34:46Z")
**Posts on this page:** 10
**Page:** 2
### Author:  [@jsjie](https://discourse.julialang.org/u/jsjie)
#### Post date: [February 8, 2024, 3:11am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/21 "2024-02-08T03:11:15Z")
If you are simply diagonalizing matrices, then iterative techniques should normally be of not much help. It is a must in planewave-DFT because the matrix can be as large as 1m \* 1m, while you only need 1k eigenvectors.
And, in most iterative techniques, you will make a reduced space, which is still larger than the number you need (most of the time it will be twice as large), and then perform direct diagonalization in the reduced space. Because of this, it may be a waste of time to calculate all eigenvectors with `zhegv`. Instead, with `zhegvx` you can set the number of eigenvectors you need. `zhegvd` will also calculate all eigenvectors, but it uses a divide-and-conquer algorithm and is generally faster than `zhegv`.
I’m using the `gv` variants because we need that in DFT, but it should be the same for the `ev` variants. That being said, you may want to have a try to `zheevd`, and even `zheevx` if you need only some of the eigenvectors. However, if your matrix is small, there should be little difference.
---
### Author:  [@jsjie](https://discourse.julialang.org/u/jsjie)
#### Post date: [February 8, 2024, 3:19am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/22 "2024-02-08T03:19:38Z")
If the assignment of off-diagonal parts are the most time-consuming, it may be helpful to rearrange arr\_1 so that it can be achieved as arr\_i[1, j, i]. But since it’s only approximately 10 \* 10 \* 3, it should have little impact.
And I’ve tried the following:
```julia
M = 10; N = 1000; y = zeros(ComplexF64, M*N, M*N); a = zeros(Int, M, M); x = zeros(ComplexF64, N, N); rand!(x);
for i = 1 : M * M
a[i] = Int(round(rand()*5))
end
function set_1!(y, a, x, M, N)
for i = 1 : M, j = 1 : i
y[1+N*(i-1):N*i, 1+N*(j-1):N*j] = a[i, j] .* x
end
end
function set_2!(y, a, x, M, N)
for i = 1 : M, j = 1 : i
y[1+N*(i-1):N*i, 1+N*(j-1):N*j] .= a[i, j] .* x
end
end
```
and get the following results:
```julia
julia> @btime set_1!(y, a, x, M, N)
338.428 ms (110 allocations: 839.24 MiB)
julia> @btime set_2!(y, a, x, M, N)
90.999 ms (0 allocations: 0 bytes)
```
I then tried directly calling BLAS:
```julia
function set_4!(y, a, x, M, N)
for p = 1 : N
ptr_x = pointer(x, 1+N*(p-1))
for i = 1 : M, j = 1 : i
coeff = a[i, j]
ptr_y = pointer(y, 1+N*M*(N*(j-1)+p-1)+N*(i-1))
ccall((BLAS.@blasfunc("zaxpby_"), BLAS.libblas), Cvoid, (Ref{BlasInt}, Ref{ComplexF64}, Ptr{ComplexF64}, Ref{BlasInt}, Ref{ComplexF64}, Ptr{ComplexF64}, Ref{BlasInt}), N, coeff, ptr_x, 1, zero(ComplexF64), ptr_y, 1)
end
end
end
```
and the result is
```julia
julia> @btime set_4!(y, a, x, M, N)
69.349 ms (0 allocations: 0 bytes)
```
why 4? because `set_3!` is just like it, but `p = 1 : N` is the innermost loop. It is slightly slower then `set_2!`, so I don’t put it here.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [February 8, 2024, 3:49am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/23 "2024-02-08T03:49:16Z")
One minor note that `Int(round(rand()*5))` could just be `rand(0:5)`
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [February 8, 2024, 11:22am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/24 "2024-02-08T11:22:04Z")
It appears to me that `mat_h_plus!` and `mat_h_minus!` construct spin raising and lowering operators for spins of size `M` (the factor of `im` there seems a bit non-standard to me and `sqrt(ind)` is not quite right, but I think you simplified for the purpose of posting the code). Anyways the point is `mat_plus` and `mat_minus` are rather sparse with (M-1)/M^2 \propto \mathcal{O}(1/M) entries. So `mat_x` and `mat_y` are similarly sparse. The following code
```julia
for i = 1:M, j=1:i
res[1+N*(i-1):N*i, 1+N*(j-1):N*j] = mat_x .* arr_1[i,j,1] .+ mat_y .* arr_1[i,j,2] .+ coeff_b .* arr_1[i,j,3]
end
```
essentially performs a few [Kronecker products](https://en.wikipedia.org/wiki/Kronecker_product) and should be equivalent to
```julia
res = @views kron(arr_1[:,:,1], mat_x) + kron(arr_1[:,:,2], mat_y) + kron(arr_1[:,:,3], coeff_b * I(M))
```
The second loop can be simplified similarly and should be equivalent to:
```julia
res += kron(Diagonal(mat_1), I(M)) + I(M*N)*coeff_a
```
So for better performance, I would use sparse matrices. Here is an attempt:
```julia
using SparseArrays
# HELPER FUNCTIONS
function mat_h_plus(size, coeff::Real)
I = 1:size-1
J = 2:size
V = im*coeff*sqrt.(I)
return sparse(I,J,V, size,size) # constructs a sparse size x size matrix
# size_o = size(original, 1)
# for ind in 1:size_o-1
# original[ind,ind+1] = im * coeff * sqrt(ind) # procedure like this
# end
end
function mat_h_minus() # similar to mat_h_plus
end
# MAIN FUNCTION
function make_matrix(
N::Integer, # large integer (~1000)
M::Integer, # small integer (~10)
mat_1::Vector{ComplexF64}, # size M
arr_1::Array{ComplexF64, 3}, # MxMx3 array
coeff::Real # coefficients for mat_hs
)
res = zeros(ComplexF64, M*N, M*N)
mat_plus = mat_h_plus(M, coeff)
mat_minus = mat_h_minus(M, coeff) # actually shouldn't that just be the Hermitian conjugate of mat_plus?
# mat_minus = mat_plus' # hermitian conjugate
mat_x = (mat_plus - mat_minus)/sqrt(2) # something like this
mat_y = # similar to calculation above
# diag = coeff_a * I(M) # diagonal matrix
# assign values for off-diagonal parts
res .+= @views kron(arr_1[i,j,1], mat_x)
res .+= @views kron(arr_1[i,j,2], mat_y)
res .+= @views kron(arr_1[i,j,3], coeff_b*I(M))
# for i = 1:M, j=1:i
# (
# res[1+N*(i-1):N*i, 1+N*(j-1):N*j]
# = mat_x .* arr_1[i,j,1] .+ mat_y .* arr_1[i,j,2] .+ coeff_b .* arr_1[i,j,3]
# )
# end
# diagonal parts
res .+= kron(Diagonal(mat_1), I(M))
res .+= coeff_a*I(M*N)
# for i=1:M
# (
# res[1+N*(i-1):N*i, 1+N*(i-1):N*i]
# .+= diag .+ mat_1[i] * I(M)
# )
# end
return Hermitian(res)
end
```
The next low-hanging fruit here is to move the allocation of `res` outside of `make_matrix` to only allocate it once per thread as I wrote [earlier](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/13).
Another improvement would be use a better in-place addition by replacing the `.+=` with a custom addition function like:
```julia
function myadd!(dense::Matrix, sparse::SparseMatrixCSC)
for col in 1:size(sparse, 2)
for r in nzrange(sparse, col)
dense[rowvals(sparse)[r], col] += nonzeros(sparse)[r]
end
end
return dense
end
```
---
### Author:  [@mdyuki](https://discourse.julialang.org/u/mdyuki)
#### Post date: [February 10, 2024, 6:12am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/25 "2024-02-10T06:12:25Z")
Thank you.
> [@jsjie](#):
>
> If you are simply diagonalizing matrices, then iterative techniques should normally be of not much help. It is a must in planewave-DFT because the matrix can be as large as 1m \* 1m, while you only need 1k eigenvectors.
> [@jsjie](#):
>
> you may want to have a try to `zheevd`, and even `zheevx` if you need only some of the eigenvectors.
For our final result, we should obtain all the eigenvalues of given matrices. So it might not be useful to use iterative techinques or subroutines like `zhegvx`.
> [@jsjie](#):
>
> it may be helpful to rearrange arr\_1 so that it can be achieved as arr\_i[1, j, i]
Actually, the `arr_1[:,:,i]` is passed by a function in a form of matrices. Rearranging them might increase numbers of assignment.
I will try your last idea to call blas subroutines directry. Thank you.
---
### Author:  [@mdyuki](https://discourse.julialang.org/u/mdyuki)
#### Post date: [February 10, 2024, 6:27am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/26 "2024-02-10T06:27:57Z")
Thank you for the reply.
> [@abraemer](#):
>
> The following code
>
> ```julia
> for i = 1:M, j=1:i
> res[1+N*(i-1):N*i, 1+N*(j-1):N*j] = mat_x .* arr_1[i,j,1] .+ mat_y .* arr_1[i,j,2] .+ coeff_b .* arr_1[i,j,3]
> end
>
> ```
>
> essentially performs a few [Kronecker products](https://en.wikipedia.org/wiki/Kronecker_product) and should be equivalent to
>
> ```julia
> res = @views kron(arr_1[:,:,1], mat_x) + kron(arr_1[:,:,2], mat_y) + kron(arr_1[:,:,3], coeff_b * I(M))
>
> ```
> [@abraemer](#):
>
> So for better performance, I would use sparse matrices.
I tried three versions from your advices.
1. original code. it takes 13 minutes, and make 38 GB allocations.
2. use kronecker product. 10 minutes, 69 GB
3. kronecker product + sparse arrays. 9 minutes, 69 GB
4. only sparse arrays. 9 minutes, 22 GB
Considering those experiments, I decided to use sparse arrrays, and not to use kronecker product methods. It reduces roughly 30 % of time and 40 % of memory allocations.
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [February 10, 2024, 6:34am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/27 "2024-02-10T06:34:47Z")
When you say “without kronecker products” you mean you build the sparse array directly from a loop? That would of course be more efficient but also a bit more difficult, so I didn’t propose it originally. If you post the current code again, I’ll have another look and see whether we can make it even faster 🙂
Edit: Does M stay constant the whole time? If so it would make sense to construct the matrices `mat_plus`,… only once and not on every call to `make_matrix`.
---
### Author:  [@mdyuki](https://discourse.julialang.org/u/mdyuki)
#### Post date: [February 10, 2024, 7:01am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/28 "2024-02-10T07:01:32Z")
I wanted to say that I don’t compute elements of `res_mat` with kronecker product.
> [@abraemer](#):
>
> Does M stay constant the whole time? If so it would make sense to construct the matrices `mat_plus`,… only once and not on every call to `make_matrix`.
Certainly, M is fixed for whole the iterations, but make\_matrix has other argument called `x`, which influences the mat\_plus through `coeff` argument. So I should make mat\_plus for every iterations.
Below is my current code for make\_matrix:
```julia
using SparseArrays
# HELPER FUNCTIONS
function mat_h_plus(size, coeff::Real)
I = 1:size-1
J = 2:size
V = coeff*sqrt.(I)
return sparse(I,J,V, size,size) # constructs a sparse size x size matrix
end
# MAIN FUNCTION
function make_matrix(
N::Integer, # large integer (~1000)
M::Integer, # small integer (~10)
mat_1::Vector{ComplexF64}, # size M
arr_1::Array{ComplexF64, 3}, # MxMx3 array
coeff::Real # coefficients for mat_hs
)
res = zeros(ComplexF64, M*N, M*N)
mat_plus = mat_h_plus(M, coeff)
mat_minus = adjoint(mat_plus) # transpose conjugate of mat_plus
mat_x = (mat_plus - mat_minus)/sqrt(2) # something like this
mat_y = # similar to calculation above
# assign values for off-diagonal parts
for i = 1:M, j=1:i
(
res[1+N*(i-1):N*i, 1+N*(j-1):N*j]
= mat_x .* arr_1[i,j,1] .+ mat_y .* arr_1[i,j,2] .+ coeff_b .* arr_1[i,j,3]
)
end
# diagonal parts
for i=1:M
(
res[1+N*(i-1):N*i, 1+N*(i-1):N*i]
.+= coef.*I(M) .+ mat_1[i] * I(M)
)
end
return Hermitian(res)
end
```
And this is main part:
```julia
# calculate needed values for x
function calc_function(M, N, x)
kvalue_range = range(0, 0.5, 201)
result_vector = zeros(201)
mat_1, arr_1 = # diagonalization of size M(~10)
coeff = coeff_func(x) # some arithmetics
Threads.@threads for k_ind in 1:201
kvalue = kvalue_range[k_ind]
energy = eigvals!( make_matrix(N, M, mat_1, arr_1, mat_1, coeff)
# we need different matrices for each kvalues
result_vector[k_ind] = sub_calculation(energy)
end
result = integrate(result_vector)
# integrate uses the trapezoidal rule: just adding values
retun result
end
const dx = 0.01 # for differentiation
container = zeros(201)
for ind in eachindex(container)
x = 0.1 * (ind-1)
Res_minus = calc_func(x-dx)
Res_plus = calc_func(x+dx)
container[ind] = (Res_plus - Res_minus) / (2*dx)
end
# some IOs
```
---
### Author:  [@abraemer](https://discourse.julialang.org/u/abraemer)
#### Post date: [February 12, 2024, 6:38am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/29 "2024-02-12T06:38:59Z")
> [@mdyuki](#):
>
> Certainly, M is fixed for whole the iterations, but make\_matrix has other argument called `x`, which influences the mat\_plus through `coeff` argument. So I should make mat\_plus for every iterations.
Note though that in this code. `mat_plus`, `mat_minus` and thus also `mat_x` and `mat_y` are directly proportional to `coeff`. So you could still build them once and then reuse. Also you really should reuse the space for `res` in `make_matrix` [like I suggested earlier](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/13).
I have the strong feeling that the code you posted is not your true code (spelling mistakes, syntax errors,…), so I won’t make the effort of restructuring it, because I don’t know whether these transformations would apply to your true code. If you post what you real code does here, then I’ll gladly show you how to restructure it to gain more performance.
---
### Author:  [@mdyuki](https://discourse.julialang.org/u/mdyuki)
#### Post date: [February 12, 2024, 8:31am UTC](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604/30 "2024-02-12T08:31:30Z")
@abraemer
I’m sorry, but it is hard to provide a more detailed code example since the real code involves information we have not published yet.
Thank you so much for all your advice. It was helpful for more efficient calculation by methods other than parallelization, which was my first purpose.
[Previous page](https://discourse.julialang.org/t/how-to-implement-hybrid-parallelized-programs/109604.md?page=1)