# Normalize a large matrix by row

**URL:** <https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997>\
**Category:** GPU\
**Tags:** cuda\
**Created:** [August 20, 2023, 7:05am UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997 "2023-08-20T07:05:49Z")\
**Posts on this page:** 11\
**Page:** 1

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 20, 2023, 7:05am UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/1 "2023-08-20T07:05:49Z")

</div>

How can I use a CUDA kernel to normalize a big matrix A by row and use the same scaling factors to a vector b?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [August 21, 2023, 3:40pm UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/2 "2023-08-21T15:40:41Z")

</div>

I don’t have a lot of practical experience but I guess you could use a `Diagonal` matrix with the normalization constants in both cases. CUDA has fast linear algebra kernels (cuBLAS)

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 22, 2023, 12:44am UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/3 "2023-08-22T00:44:43Z")

</div>

Thanks, @josuagrw  
First I tried

```julia
function row_scale_1_kernel(
    scales::CuDeviceVector{Tv},
    rowPtr::CuDeviceVector{Ti},
    colVal::CuDeviceVector{Ti},
    nzVal::CuDeviceVector{Tv}
) where {Tv, Ti}
    row = threadIdx().x + (blockIdx().x - 1) * blockDim().x
    if row <= length(rowPtr) - 1 && row <= size(scales, 1)
        row_start = rowPtr[row]
        row_end = rowPtr[row + 1] - 1
        sum = zero(Tv)
        for j = row_start:row_end
            sum += nzVal[j] * nzVal[j]
        end
        scales[row] = 1.0 / sqrt(sum)
    end
    return
end

function A_matrix(A_cpu::SparseArrays.SparseMatrixCSC{Float64,Int})
    A = CuSparseMatrixCSR(A_cpu)
    rows, cols = size(A) # A is in CSR format in GPU

    scale_constants = CUDA.zeros(Float64, size(A, 1))
    rowPtr = A.rowPtr
    colVal = A.colVal
    nzVal = A.nzVal

    num_rows = size(A, 1) #Number of rows in the matrix
    threads_per_block = 256 # Number of threads per block (can be 256 or any other suitable value)
    blocks = cld(num_rows, threads_per_block) # Number of blocks needed to cover all rows
    
    @cuda threads=threads_per_block blocks=blocks row_scale_1_kernel(
        scale_constants,
        rowPtr,
        colVal,
        nzVal,
    )
    
    D = Diagonal(scale_constants)
    A = D * A

    a_i = Array(rowPtr)
    a_j = Array(colVal)
    a_v = Array(nzVal)
    scale = Array(scale_constants)

    return a_i, a_j, a_v, scale, rows, cols
end

```

It returns the right values, but before Julia stops, I get tones of

```julia
WARNING: Error while freeing DeviceBuffer(37.820 KiB at 0x0000000302000000):
CUDA.CuError(code=CUDA.cudaError_enum(0x000002c5), meta=nothing)

```

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [August 22, 2023, 9:05am UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/4 "2023-08-22T09:05:30Z")

</div>

I’m not qualified to comment on this pointer arithmetic.

But are you certain you cannot do this using the library functions? It looks like all you need is

1. a dense element-wise multiply to square the values
2. a sparse matrix vector multiply to sum each row

Both of these definitely exist in cuSparse and cuBLAS. I would not reinvent the wheel here…

---

<div class="post-metadata">

**Author:** ![Per](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/per/32/10387_2.png) [@Per](https://discourse.julialang.org/u/Per)\
**Post date:** [August 22, 2023, 9:45am UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/5 "2023-08-22T09:45:24Z")

</div>

Does it have to be a kernel? It could easily be done via broadcast. Something like:

```julia
n = sqrt.(sum(abs2, A, dims=2)))
A ./= n
b ./= n

```

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 23, 2023, 11:30pm UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/6 "2023-08-23T23:30:43Z")

</div>

@Per

`n = sqrt.(sum(abs2, A, dims=2))` causes scalar indexing of a GPU array, so I used Kernel.

Try this code

```julia
using CUDA
using SparseArrays
using CUDA.CUSPARSE
using LinearAlgebra

n = 1000
A_sparse = sprand(n, n, 0.5)
A = CuSparseMatrixCSR(A_sparse)
n = sqrt.(sum(abs2, A, dims=2))

```

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [August 24, 2023, 6:05am UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/7 "2023-08-24T06:05:50Z")

</div>

Looks like support for reductions like `sum(abs2, A; dims=2)` for CUSPARSE arrays was recently implemented but not yet released: [CUSPARSE: Implement out-of-place reductions (#1987) · JuliaGPU/CUDA.jl@7660edc · GitHub](https://github.com/JuliaGPU/CUDA.jl/commit/7660edc616d841b76c74e9489b71d37907d950f2)

So you can either wait for the next release of CUDA.jl (ping @maleadt), or repurpose some of the code removed in that commit, specifically the functions highlighted here: [CUSPARSE: Implement out-of-place reductions (#1987) · JuliaGPU/CUDA.jl@7660edc · GitHub](https://github.com/JuliaGPU/CUDA.jl/commit/7660edc616d841b76c74e9489b71d37907d950f2#diff-e1fa24914853993ba813e2b2ee908a06f18ebcf7282c9dc9fc3e47842e91aac5L4-L49)

---

<div class="post-metadata">

**Author:** ![maleadt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maleadt/32/10097_2.png) [@maleadt](https://discourse.julialang.org/u/maleadt)\
**Post date:** [August 24, 2023, 6:56am UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/8 "2023-08-24T06:56:47Z")

</div>

> [@danielwe](#):
>
> repurpose some of the code removed in that commit, specifically the functions highlighted here: [CUSPARSE: Implement out-of-place reductions (#1987) · JuliaGPU/CUDA.jl@7660edc · GitHub](https://github.com/JuliaGPU/CUDA.jl/commit/7660edc616d841b76c74e9489b71d37907d950f2#diff-e1fa24914853993ba813e2b2ee908a06f18ebcf7282c9dc9fc3e47842e91aac5L4-L49)

That’s the removed code; the implementation of `mapreduce` below should be more generic.

> [@danielwe](#):
>
> wait for the next release of CUDA.jl

Just FYI, there’s a couple of minor but breaking changes incoming (such as the deprecation of Julia 1.6-1.7 and CUDA 11.0-11.3), so the next release will take a while. If you desperately need a feature in a release, we can consider backporting it to the previous stable one.

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 24, 2023, 8:13pm UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/9 "2023-08-24T20:13:32Z")

</div>

Thanks @maleadt

I really need that functionality because CPU calculation is very lengthy.

---

<div class="post-metadata">

**Author:** ![A\_M](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/a_m/32/46038_2.png) [@A\_M](https://discourse.julialang.org/u/A_M)\
**Post date:** [August 24, 2023, 9:11pm UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/10 "2023-08-24T21:11:02Z")

</div>

I just added the `lp/sum-dim-fix` branch, I could use it like

```julia
using CUDA
using CUDA.CUSPARSE
using SparseArrays
using LinearAlgebra

function scale!(coefficients)
    A = CuSparseMatrixCSR(convert(SparseArrays.SparseMatrixCSC{Float64,Int}, coefficients))
    scaling_factors = 1.0 ./ sqrt.(sum(abs2, A, dims=2))
    D = Diagonal(scaling_factors)
    A = D * A
    return A
end

m,n = 5,6
p = 0.5
x = sprand(Float64, m, n, p)
scale!(x)

```

But I am not sure; when I use the same function inside my code, I get the same error

```julia
WARNING: Error while freeing DeviceBuffer(223.113 KiB at 0x0000000302000000):
CUDA.CuError(code=CUDA.cudaError_enum(0x000002c5), meta=nothing)

```

---

<div class="post-metadata">

**Author:** ![maleadt](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/maleadt/32/10097_2.png) [@maleadt](https://discourse.julialang.org/u/maleadt)\
**Post date:** [August 25, 2023, 8:21am UTC](https://discourse.julialang.org/t/normalize-a-large-matrix-by-row/102997/11 "2023-08-25T08:21:51Z")

</div>

That error is unrelated. The exception, `CUDA_ERROR_CONTEXT_IS_DESTROYED`, indicates that some memory is incorrectly being freed during process teardown. If you have an MWE, please file an issue.