rocSPARSE examples#
rocSPARSE is the AMD
implementation of sparse linear algebra for AMD GPUs. hipFORT exposes it through
the hipfort_rocsparse module, which mirrors the rocSPARSE C API one to one.
rocSPARSE reuses the rocBLAS handle type, so the examples that create a handle
call rocsparse_create_handle from the same module.
Every program on this page is a complete, self-contained example that is built
and run as part of the hipFORT test suite. The Fortran 2008 sources live in
test/f2008/rocsparse and the equivalent Fortran 2003 sources, which use
type(c_ptr) device pointers and explicit byte counts instead of Fortran
array pointers, live in test/f2003/rocsparse.
Where a routine has the four precisions, the example is provided for each: s
(real single), d (real double), c (complex single), and z (complex
double). This page shows the double-precision program of each group; the other
precisions differ only in the host data type and the rocsparse_ prefix
letter.
Conventions#
rocSPARSE follows a small number of conventions that recur in every example:
Sparse matrix formats. Most examples store the sparse matrix in CSR (compressed sparse row): a row-pointer array, a column-index array, and a values array. Block variants use BSR, and a few routines take COO (coordinate) row/column arrays.
Zero-based indexing. The examples use
rocsparse_index_base_zero, so CSR row pointers and column indices start at 0, matching the C samples. The Fortran host arrays that hold them are ordinary 1-based arrays whose values are 0-based.The generic API is descriptor- and stage-based. The newer routines (SpMV, SpMM, SDDMM, SpSV, SpSM) wrap the operands in matrix/vector descriptors (
rocsparse_create_csr_descr,rocsparse_create_dnmat_descr,rocsparse_create_dnvec_descr) and run in stages: query a workspace size, optionally preprocess/analyze, then compute. The descriptor constructors arec_ptr-only (no array overloads), so device buffers are passed viac_loc(...)even in the Fortran 2008 examples.Scalars.
alphaandbetaare passed by address (c_loc(alpha)) in the generic API, and as host scalars by reference in the older level-2/level-3 routines such asbsrmvandgemvi.Every call returns a status code. The examples wrap rocSPARSE calls in
rocsparseCheckand HIP calls inhipCheckfrom thehipfort_checkmodule, both of which abort on failure.
Building an example#
The examples need the rocsparse, rocblas, and hip hipFORT
components:
find_package(hipfort REQUIRED COMPONENTS hip rocblas rocsparse)
add_executable(my_sparse rocsparse_dspmv.f08)
target_link_libraries(my_sparse PRIVATE hipfort::rocsparse hipfort::rocblas hipfort::hip)
See Using hipFORT in your application for the full set of build options.
Sparse matrix-vector and matrix-matrix products#
spmv multiplies a sparse matrix by a dense vector, y = alpha*A*x + beta*y,
using the generic API: a CSR descriptor for A and dense-vector descriptors
for x and y, run through the buffer_size and compute stages.
!!!!!!!!!!!!!/
! dspmv example (double-precision sparse-matrix dense-vector multiply, y = alpha*A*x + beta*y)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/generic.html
!
! Uses the generic API: build a CSR descriptor for A and dense-vector
! descriptors for x and y, then run the three spmv stages
! (buffer_size -> preprocess -> compute). Result is checked against A*x.
!
! NOTE: the descriptor constructors are c_ptr-only (no array overloads), so
! device buffers are passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dspmv
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
! Sparse A (3x3) in CSR (0-based): A = [[1,0,2],[0,3,0],[4,0,5]]
integer(c_int), parameter :: M = 3, N = 3, nnz = 5
integer(c_int) :: h_csr_row_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_csr_col_ind(5) = (/0, 2, 1, 0, 2/)
real(c_double) :: h_csr_val(5) = (/1, 2, 3, 4, 5/)
real(c_double) :: h_x(3) = (/1, 2, 3/)
real(c_double) :: h_y(3)
real(c_double) :: h_expected(3) = (/7, 6, 19/) ! A*x
real(c_double), target :: alpha = 1.0_c_double, beta = 0.0_c_double
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:)
real(c_double), pointer :: d_x(:), d_y(:)
type(c_ptr) :: handle, matA, vecX, vecY, d_buffer
integer(c_size_t) :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dspmv' (Fortran 2008 interfaces) - "
! Allocate device memory and copy inputs to device
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_csr_col_ind, source=h_csr_col_ind))
call hipCheck(hipMalloc(d_csr_val, source=h_csr_val))
call hipCheck(hipMalloc(d_x, source=h_x))
call hipCheck(hipMalloc(d_y, mold=h_y))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
! Descriptors: CSR for A, dense vectors for x and y
call rocsparseCheck(rocsparse_create_csr_descr(matA, int(M,c_int64_t), int(N,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), c_loc(d_csr_val), &
rocsparse_indextype_i32, rocsparse_indextype_i32, &
rocsparse_index_base_zero, rocsparse_datatype_f64_r))
call rocsparseCheck(rocsparse_create_dnvec_descr(vecX, int(N,c_int64_t), c_loc(d_x), rocsparse_datatype_f64_r))
call rocsparseCheck(rocsparse_create_dnvec_descr(vecY, int(M,c_int64_t), c_loc(d_y), rocsparse_datatype_f64_r))
! Stage 1: query workspace size
call rocsparseCheck(rocsparse_spmv(handle, rocsparse_operation_none, c_loc(alpha), matA, vecX, c_loc(beta), vecY, &
rocsparse_datatype_f64_r, rocsparse_spmv_alg_default, &
rocsparse_spmv_stage_buffer_size, buffer_size, c_null_ptr))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Stage 2: preprocess
call rocsparseCheck(rocsparse_spmv(handle, rocsparse_operation_none, c_loc(alpha), matA, vecX, c_loc(beta), vecY, &
rocsparse_datatype_f64_r, rocsparse_spmv_alg_default, &
rocsparse_spmv_stage_preprocess, buffer_size, d_buffer))
! Stage 3: compute
call rocsparseCheck(rocsparse_spmv(handle, rocsparse_operation_none, c_loc(alpha), matA, vecX, c_loc(beta), vecY, &
rocsparse_datatype_f64_r, rocsparse_spmv_alg_default, &
rocsparse_spmv_stage_compute, buffer_size, d_buffer))
! Copy result back to host
call hipCheck(hipMemcpy(h_y, d_y, hipMemcpyDeviceToHost))
! Verify y == A*x
do i = 1,M
error = abs(h_y(i) - h_expected(i)) / max(abs(h_expected(i)), 1.0_c_double)
if(error .gt. error_max) then
write(*,*) "FAILED! Error bigger than max! Error = ", error, " y(", i, ") = ", h_y(i)
call exit
end if
end do
! Clean up
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_csr_col_ind))
call hipCheck(hipFree(d_csr_val))
call hipCheck(hipFree(d_x))
call hipCheck(hipFree(d_y))
call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program dspmv
spmm multiplies a sparse matrix by a dense matrix, C = alpha*A*B + beta*C,
with dense-matrix descriptors for B and C and the three
buffer_size / preprocess / compute stages.
!!!!!!!!!!!!!/
! dspmm example (double-precision sparse-matrix times dense-matrix, C = alpha*A*B + beta*C)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/generic.html
!
! Uses the generic API: build a CSR descriptor for A and dense-matrix
! descriptors for B and C, then run the three spmm stages
! (buffer_size -> preprocess -> compute). Result is checked against A*B.
!
! NOTE: the descriptor constructors are c_ptr-only (no array overloads), so
! device buffers are passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dspmm
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i, j
! Sparse A (3x3) in CSR (0-based):
! A = [[1,0,2],[0,3,0],[4,0,5]]
integer(c_int), parameter :: M = 3, K = 3, Ncol = 2, nnz = 5
integer(c_int) :: h_csr_row_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_csr_col_ind(5) = (/0, 2, 1, 0, 2/)
real(c_double) :: h_csr_val(5) = (/1, 2, 3, 4, 5/)
! Dense B (3x2), column-major
real(c_double) :: h_B(3,2) = reshape((/1, 2, 3, 4, 5, 6/), (/3, 2/))
real(c_double) :: h_C(3,2)
! Expected C = A*B
real(c_double) :: h_expected(3,2) = reshape((/7, 6, 19, 16, 15, 46/), (/3, 2/))
real(c_double), target :: alpha = 1.0_c_double, beta = 0.0_c_double
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:)
real(c_double), pointer :: d_B(:,:), d_C(:,:)
type(c_ptr) :: handle, matA, matB, matC, d_buffer
integer(c_size_t) :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dspmm' (Fortran 2008 interfaces) - "
! Allocate device memory and copy inputs to device
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_csr_col_ind, source=h_csr_col_ind))
call hipCheck(hipMalloc(d_csr_val, source=h_csr_val))
call hipCheck(hipMalloc(d_B, source=h_B))
call hipCheck(hipMalloc(d_C, mold=h_C))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
! Descriptors: CSR for A, dense (column-major) for B and C
call rocsparseCheck(rocsparse_create_csr_descr(matA, int(M,c_int64_t), int(K,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), c_loc(d_csr_val), &
rocsparse_indextype_i32, rocsparse_indextype_i32, &
rocsparse_index_base_zero, rocsparse_datatype_f64_r))
call rocsparseCheck(rocsparse_create_dnmat_descr(matB, int(K,c_int64_t), int(Ncol,c_int64_t), int(K,c_int64_t), &
c_loc(d_B), rocsparse_datatype_f64_r, rocsparse_order_column))
call rocsparseCheck(rocsparse_create_dnmat_descr(matC, int(M,c_int64_t), int(Ncol,c_int64_t), int(M,c_int64_t), &
c_loc(d_C), rocsparse_datatype_f64_r, rocsparse_order_column))
! Stage 1: query workspace size
call rocsparseCheck(rocsparse_spmm(handle, rocsparse_operation_none, rocsparse_operation_none, c_loc(alpha), &
matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, rocsparse_spmm_alg_default, &
rocsparse_spmm_stage_buffer_size, buffer_size, c_null_ptr))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Stage 2: preprocess
call rocsparseCheck(rocsparse_spmm(handle, rocsparse_operation_none, rocsparse_operation_none, c_loc(alpha), &
matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, rocsparse_spmm_alg_default, &
rocsparse_spmm_stage_preprocess, buffer_size, d_buffer))
! Stage 3: compute
call rocsparseCheck(rocsparse_spmm(handle, rocsparse_operation_none, rocsparse_operation_none, c_loc(alpha), &
matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, rocsparse_spmm_alg_default, &
rocsparse_spmm_stage_compute, buffer_size, d_buffer))
! Copy result back to host
call hipCheck(hipMemcpy(h_C, d_C, hipMemcpyDeviceToHost))
! Verify C == A*B
do j = 1,Ncol
do i = 1,M
error = abs(h_C(i,j) - h_expected(i,j)) / max(abs(h_expected(i,j)), 1.0_c_double)
if(error .gt. error_max) then
write(*,*) "FAILED! Error bigger than max! Error = ", error, " at (", i, ",", j, ")"
call exit
end if
end do
end do
! Clean up
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_csr_col_ind))
call hipCheck(hipFree(d_csr_val))
call hipCheck(hipFree(d_B))
call hipCheck(hipFree(d_C))
call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program dspmm
Sampled dense-dense matrix multiplication#
sddmm is the transpose of the SpMM data flow: the dense product A*B is
evaluated only at the nonzero positions of a sparse C, giving
C = alpha * (A*B) .* spy(C) + beta*C. It is the core primitive behind
attention and graph-neural-network kernels. The example uses dense descriptors
for A and B, a CSR descriptor for C, and the three sddmm stages.
!!!!!!!!!!!!!/
! dsddmm example (double-precision sampled dense-dense matrix multiplication)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/generic.html
!
! SDDMM computes C = alpha * (A * B) .* spy(C) + beta * C, where A and B are
! dense and C is sparse (CSR): the dense product A*B is only evaluated at the
! nonzero positions of C. Uses the generic API with the three sddmm stages
! (buffer_size -> preprocess -> compute). Result is checked against a host
! reference that samples matmul(A,B) at C's pattern.
!
! NOTE: the descriptor constructors are c_ptr-only (no array overloads), so
! device buffers are passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dsddmm
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
! Dense A (M x K) and B (K x N), column-major. Sparse C (M x N) in CSR.
integer(c_int), parameter :: M = 3, N = 2, K = 3, nnz = 4
! A = [[1,2,3],[4,5,6],[7,8,10]] (column-major storage)
real(c_double) :: h_A(3,3) = reshape((/1, 4, 7, 2, 5, 8, 3, 6, 10/), (/3,3/))
! B = [[1,2],[3,4],[5,6]] (column-major storage)
real(c_double) :: h_B(3,2) = reshape((/1, 3, 5, 2, 4, 6/), (/3,2/))
! Sparse C pattern (0-based CSR), C = [[*,0],[0,*],[*,*]]: nonzeros at
! (0,0),(1,1),(2,0),(2,1). Values start at 0 (beta = 0 so they are ignored).
integer(c_int) :: h_csr_row_ptr(4) = (/0, 1, 2, 4/)
integer(c_int) :: h_csr_col_ind(4) = (/0, 1, 0, 1/)
real(c_double) :: h_csr_val(4) = (/0, 0, 0, 0/)
real(c_double) :: h_AB(3,2), h_expected(4)
real(c_double), target :: alpha = 1.0_c_double, beta = 0.0_c_double
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:)
real(c_double), pointer :: d_A(:,:), d_B(:,:)
type(c_ptr) :: handle, matA, matB, matC, d_buffer
integer(c_size_t) :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dsddmm' (Fortran 2008 interfaces) - "
! Host reference: sample matmul(A,B) at C's nonzero pattern (row-major walk of CSR)
h_AB = matmul(h_A, h_B)
h_expected(1) = h_AB(1,1) ! (0,0)
h_expected(2) = h_AB(2,2) ! (1,1)
h_expected(3) = h_AB(3,1) ! (2,0)
h_expected(4) = h_AB(3,2) ! (2,1)
! Allocate device memory and copy inputs to device
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_csr_col_ind, source=h_csr_col_ind))
call hipCheck(hipMalloc(d_csr_val, source=h_csr_val))
call hipCheck(hipMalloc(d_A, source=h_A))
call hipCheck(hipMalloc(d_B, source=h_B))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
! Descriptors: dense (column-major) for A and B, CSR for C
call rocsparseCheck(rocsparse_create_dnmat_descr(matA, int(M,c_int64_t), int(K,c_int64_t), int(M,c_int64_t), &
c_loc(d_A), rocsparse_datatype_f64_r, rocsparse_order_column))
call rocsparseCheck(rocsparse_create_dnmat_descr(matB, int(K,c_int64_t), int(N,c_int64_t), int(K,c_int64_t), &
c_loc(d_B), rocsparse_datatype_f64_r, rocsparse_order_column))
call rocsparseCheck(rocsparse_create_csr_descr(matC, int(M,c_int64_t), int(N,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), c_loc(d_csr_val), &
rocsparse_indextype_i32, rocsparse_indextype_i32, &
rocsparse_index_base_zero, rocsparse_datatype_f64_r))
! Stage 1: query workspace size
call rocsparseCheck(rocsparse_sddmm_buffer_size(handle, rocsparse_operation_none, rocsparse_operation_none, &
c_loc(alpha), matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, &
rocsparse_sddmm_alg_default, buffer_size))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Stage 2: preprocess
call rocsparseCheck(rocsparse_sddmm_preprocess(handle, rocsparse_operation_none, rocsparse_operation_none, &
c_loc(alpha), matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, &
rocsparse_sddmm_alg_default, d_buffer))
! Stage 3: compute
call rocsparseCheck(rocsparse_sddmm(handle, rocsparse_operation_none, rocsparse_operation_none, &
c_loc(alpha), matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, &
rocsparse_sddmm_alg_default, d_buffer))
! Copy the sparse values back to host
call hipCheck(hipMemcpy(h_csr_val, d_csr_val, hipMemcpyDeviceToHost))
! Verify sampled values
do i = 1,nnz
error = abs(h_csr_val(i) - h_expected(i)) / max(abs(h_expected(i)), 1.0_c_double)
if(error .gt. error_max) then
write(*,*) "FAILED! Error bigger than max! Error = ", error, " at nnz ", i, &
" got ", h_csr_val(i), " expected ", h_expected(i)
call exit
end if
end do
! Clean up
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_csr_col_ind))
call hipCheck(hipFree(d_csr_val))
call hipCheck(hipFree(d_A))
call hipCheck(hipFree(d_B))
call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program dsddmm
rocSPARSE also supports a batched SDDMM: the dense operands are strided-batched
with rocsparse_dnmat_set_strided_batch and the sparse C shares one
sparsity pattern across the batch with rocsparse_csr_set_strided_batch.
!!!!!!!!!!!!!/
! dsddmm_batched example (double-precision batched sampled dense-dense matmul)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/generic.html
!
! Batched SDDMM: for each batch b, C_b = alpha * (A_b * B_b) .* spy(C) + beta*C_b.
! The dense matrices A and B are strided-batched (rocsparse_dnmat_set_strided_batch);
! the sparse C shares one CSR sparsity pattern across the batch
! (offsets_batch_stride = 0) with per-batch values (columns_values_batch_stride
! = nnz), set with rocsparse_csr_set_strided_batch. Each batch's sampled values
! are checked against matmul(A_b, B_b).
!
! NOTE: the descriptor constructors are c_ptr-only (no array overloads), so
! device buffers are passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dsddmm_batched
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: b, i
integer(c_int), parameter :: M = 3, N = 2, K = 3, nnz = 4, batch = 2
! Two dense A matrices (M x K each), column-major, stored back-to-back.
! A_1 = [[1,2,3],[4,5,6],[7,8,10]] A_2 = 2 * A_1
real(c_double) :: h_A(3,3,2) = reshape((/ &
1, 4, 7, 2, 5, 8, 3, 6, 10, &
2, 8, 14, 4, 10, 16, 6, 12, 20/), (/3,3,2/))
! Two dense B matrices (K x N each), column-major.
! B_1 = [[1,2],[3,4],[5,6]] B_2 = B_1 + 1
real(c_double) :: h_B(3,2,2) = reshape((/ &
1, 3, 5, 2, 4, 6, &
2, 4, 6, 3, 5, 7/), (/3,2,2/))
! Shared sparse C pattern (0-based CSR), nonzeros at (0,0),(1,1),(2,0),(2,1).
integer(c_int) :: h_csr_row_ptr(4) = (/0, 1, 2, 4/)
integer(c_int) :: h_csr_col_ind(4) = (/0, 1, 0, 1/)
real(c_double) :: h_csr_val(nnz*batch) = 0 ! batch-major: [b1 nnz | b2 nnz]
real(c_double) :: h_AB(3,2), h_expected(nnz,batch)
real(c_double), target :: alpha = 1.0_c_double, beta = 0.0_c_double
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:)
real(c_double), pointer :: d_A(:,:,:), d_B(:,:,:)
type(c_ptr) :: handle, matA, matB, matC, d_buffer
integer(c_size_t) :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dsddmm_batched' (Fortran 2008 interfaces) - "
! Host reference per batch: sample matmul(A_b, B_b) at C's pattern.
do b = 1, batch
h_AB = matmul(h_A(:,:,b), h_B(:,:,b))
h_expected(1,b) = h_AB(1,1) ! (0,0)
h_expected(2,b) = h_AB(2,2) ! (1,1)
h_expected(3,b) = h_AB(3,1) ! (2,0)
h_expected(4,b) = h_AB(3,2) ! (2,1)
end do
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_csr_col_ind, source=h_csr_col_ind))
call hipCheck(hipMalloc(d_csr_val, source=h_csr_val))
call hipCheck(hipMalloc(d_A, source=h_A))
call hipCheck(hipMalloc(d_B, source=h_B))
call rocsparseCheck(rocsparse_create_handle(handle))
! Dense descriptors on the first batch slice, then attach the batch stride.
call rocsparseCheck(rocsparse_create_dnmat_descr(matA, int(M,c_int64_t), int(K,c_int64_t), int(M,c_int64_t), &
c_loc(d_A), rocsparse_datatype_f64_r, rocsparse_order_column))
call rocsparseCheck(rocsparse_create_dnmat_descr(matB, int(K,c_int64_t), int(N,c_int64_t), int(K,c_int64_t), &
c_loc(d_B), rocsparse_datatype_f64_r, rocsparse_order_column))
call rocsparseCheck(rocsparse_create_csr_descr(matC, int(M,c_int64_t), int(N,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), c_loc(d_csr_val), &
rocsparse_indextype_i32, rocsparse_indextype_i32, &
rocsparse_index_base_zero, rocsparse_datatype_f64_r))
! Batch config: dense strides are the per-matrix element counts; the sparse C
! reuses one pattern (offsets stride 0) with per-batch values (stride nnz).
call rocsparseCheck(rocsparse_dnmat_set_strided_batch(matA, batch, int(M*K,c_int64_t)))
call rocsparseCheck(rocsparse_dnmat_set_strided_batch(matB, batch, int(K*N,c_int64_t)))
call rocsparseCheck(rocsparse_csr_set_strided_batch(matC, batch, 0_c_int64_t, int(nnz,c_int64_t)))
! Stage 1: query workspace size
call rocsparseCheck(rocsparse_sddmm_buffer_size(handle, rocsparse_operation_none, rocsparse_operation_none, &
c_loc(alpha), matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, &
rocsparse_sddmm_alg_default, buffer_size))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Stage 2: preprocess
call rocsparseCheck(rocsparse_sddmm_preprocess(handle, rocsparse_operation_none, rocsparse_operation_none, &
c_loc(alpha), matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, &
rocsparse_sddmm_alg_default, d_buffer))
! Stage 3: compute
call rocsparseCheck(rocsparse_sddmm(handle, rocsparse_operation_none, rocsparse_operation_none, &
c_loc(alpha), matA, matB, c_loc(beta), matC, rocsparse_datatype_f64_r, &
rocsparse_sddmm_alg_default, d_buffer))
call hipCheck(hipMemcpy(h_csr_val, d_csr_val, hipMemcpyDeviceToHost))
do b = 1, batch
do i = 1, nnz
error = abs(h_csr_val((b-1)*nnz+i) - h_expected(i,b)) / max(abs(h_expected(i,b)), 1.0_c_double)
if(error .gt. error_max) then
write(*,*) "FAILED! Error bigger than max! Error = ", error, " at batch ", b, " nnz ", i, &
" got ", h_csr_val((b-1)*nnz+i), " expected ", h_expected(i,b)
call exit
end if
end do
end do
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_csr_col_ind))
call hipCheck(hipFree(d_csr_val))
call hipCheck(hipFree(d_A))
call hipCheck(hipFree(d_B))
call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program dsddmm_batched
Sparse triangular solves#
sptrsv solves a sparse triangular system op(A)*y = alpha*x for a single
right-hand side. The generic API adds an analysis stage between the buffer-size
query and the solve, which inspects the sparsity pattern once and can be reused.
!!!!!!!!!!!!!/
! dsptrsv example (double-precision sparse triangular solve, op(A)*y = alpha*x)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/generic.html
!
! Uses the generic spsv API on a lower-triangular L. Self-verifying: pick a
! known y, form x = L*y, solve L*y' = x, and confirm y' recovers y.
! The triangular structure is set via spmat_set_attribute (fill_mode +
! diag_type). Three stages: buffer_size -> preprocess -> compute.
!
! NOTE: descriptor/array arguments are c_ptr-only, so device buffers and the
! attribute values are passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dsptrsv
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
! Lower-triangular L (3x3) in CSR (0-based):
! L = [[2,0,0],[1,3,0],[4,5,6]]
integer(c_int), parameter :: M = 3, N = 3, nnz = 6
integer(c_int) :: h_csr_row_ptr(4) = (/0, 1, 3, 6/)
integer(c_int) :: h_csr_col_ind(6) = (/0, 0, 1, 0, 1, 2/)
real(c_double) :: h_csr_val(6) = (/2, 1, 3, 4, 5, 6/)
real(c_double) :: h_y(3) = (/1, 2, 3/) ! known solution
real(c_double) :: h_x(3) ! rhs = L*y
real(c_double) :: h_yout(3) ! recovered solution
real(c_double), target :: alpha = 1.0_c_double
integer(kind(rocsparse_fill_mode_lower)), target :: fill = rocsparse_fill_mode_lower
integer(kind(rocsparse_diag_type_non_unit)), target :: diag = rocsparse_diag_type_non_unit
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:)
real(c_double), pointer :: d_x(:), d_y(:)
type(c_ptr) :: handle, matL, vecX, vecY, d_buffer
integer(c_size_t) :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 100 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dsptrsv' (Fortran 2008 interfaces) - "
! Build a consistent rhs so that L*y = x: x1=2, x2=1+6=7, x3=4+10+18=32
h_x(1) = 2.0_c_double
h_x(2) = 7.0_c_double
h_x(3) = 32.0_c_double
! Allocate device memory and copy inputs
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_csr_col_ind, source=h_csr_col_ind))
call hipCheck(hipMalloc(d_csr_val, source=h_csr_val))
call hipCheck(hipMalloc(d_x, source=h_x))
call hipCheck(hipMalloc(d_y, mold=h_yout))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
! CSR descriptor for L, marked lower-triangular / non-unit-diagonal
call rocsparseCheck(rocsparse_create_csr_descr(matL, int(M,c_int64_t), int(N,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), c_loc(d_csr_val), &
rocsparse_indextype_i32, rocsparse_indextype_i32, &
rocsparse_index_base_zero, rocsparse_datatype_f64_r))
call rocsparseCheck(rocsparse_spmat_set_attribute(matL, rocsparse_spmat_fill_mode, c_loc(fill), int(4,c_size_t)))
call rocsparseCheck(rocsparse_spmat_set_attribute(matL, rocsparse_spmat_diag_type, c_loc(diag), int(4,c_size_t)))
! Dense-vector descriptors for the rhs (x) and the solution (y)
call rocsparseCheck(rocsparse_create_dnvec_descr(vecX, int(M,c_int64_t), c_loc(d_x), rocsparse_datatype_f64_r))
call rocsparseCheck(rocsparse_create_dnvec_descr(vecY, int(M,c_int64_t), c_loc(d_y), rocsparse_datatype_f64_r))
! Stage 1: workspace size
call rocsparseCheck(rocsparse_spsv(handle, rocsparse_operation_none, c_loc(alpha), matL, vecX, vecY, &
rocsparse_datatype_f64_r, rocsparse_spsv_alg_default, &
rocsparse_spsv_stage_buffer_size, buffer_size, c_null_ptr))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Stage 2: preprocess (analysis)
call rocsparseCheck(rocsparse_spsv(handle, rocsparse_operation_none, c_loc(alpha), matL, vecX, vecY, &
rocsparse_datatype_f64_r, rocsparse_spsv_alg_default, &
rocsparse_spsv_stage_preprocess, buffer_size, d_buffer))
! Stage 3: solve
call rocsparseCheck(rocsparse_spsv(handle, rocsparse_operation_none, c_loc(alpha), matL, vecX, vecY, &
rocsparse_datatype_f64_r, rocsparse_spsv_alg_default, &
rocsparse_spsv_stage_compute, buffer_size, d_buffer))
! Copy the recovered solution back
call hipCheck(hipMemcpy(h_yout, d_y, hipMemcpyDeviceToHost))
! Verify y' == y
do i = 1,M
error = abs(h_yout(i) - h_y(i)) / max(abs(h_y(i)), 1.0_c_double)
if(error .gt. error_max) then
write(*,*) "FAILED! Error bigger than max! Error = ", error, " y(", i, ") = ", h_yout(i)
call exit
end if
end do
! Clean up
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_csr_col_ind))
call hipCheck(hipFree(d_csr_val))
call hipCheck(hipFree(d_x))
call hipCheck(hipFree(d_y))
call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program dsptrsv
sptrsm solves the same kind of system with several right-hand sides at once,
taking a dense-matrix descriptor for the right-hand sides.
!!!!!!!!!!!!!/
! dsptrsm example (double-precision sparse triangular solve with multiple rhs)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/generic.html
!
! Uses the generic spsm API on a lower-triangular L with a dense rhs matrix.
! Self-verifying: pick a known Y, form X = L*Y, solve L*C = X, and confirm C
! recovers Y. Triangular structure set via spmat_set_attribute. Three stages:
! buffer_size -> preprocess -> compute.
!
! NOTE: descriptor/array arguments are c_ptr-only, so device buffers and the
! attribute values are passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dsptrsm
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i, j
! Lower-triangular L (3x3) in CSR (0-based): L = [[2,0,0],[1,3,0],[4,5,6]]
integer(c_int), parameter :: M = 3, nrhs = 2, nnz = 6
integer(c_int) :: h_csr_row_ptr(4) = (/0, 1, 3, 6/)
integer(c_int) :: h_csr_col_ind(6) = (/0, 0, 1, 0, 1, 2/)
real(c_double) :: h_csr_val(6) = (/2, 1, 3, 4, 5, 6/)
real(c_double) :: h_Y(3,2) = reshape((/1, 2, 3, 4, 5, 6/), (/3, 2/)) ! known solution
real(c_double) :: h_X(3,2) ! rhs = L*Y
real(c_double) :: h_C(3,2) ! recovered solution
real(c_double) :: L_dense(3,3)
real(c_double), target :: alpha = 1.0_c_double
integer(kind(rocsparse_fill_mode_lower)), target :: fill = rocsparse_fill_mode_lower
integer(kind(rocsparse_diag_type_non_unit)), target :: diag = rocsparse_diag_type_non_unit
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:)
real(c_double), pointer :: d_X(:,:), d_C(:,:)
type(c_ptr) :: handle, matL, matB, matC, d_buffer
integer(c_size_t) :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 100 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dsptrsm' (Fortran 2008 interfaces) - "
! Build dense L and the consistent rhs X = L*Y on the host
L_dense = 0.0_c_double
L_dense(1,1) = 2.0_c_double
L_dense(2,1) = 1.0_c_double; L_dense(2,2) = 3.0_c_double
L_dense(3,1) = 4.0_c_double; L_dense(3,2) = 5.0_c_double; L_dense(3,3) = 6.0_c_double
h_X = matmul(L_dense, h_Y)
! Allocate device memory and copy inputs
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_csr_col_ind, source=h_csr_col_ind))
call hipCheck(hipMalloc(d_csr_val, source=h_csr_val))
call hipCheck(hipMalloc(d_X, source=h_X))
call hipCheck(hipMalloc(d_C, mold=h_C))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
! CSR descriptor for L, marked lower-triangular / non-unit-diagonal
call rocsparseCheck(rocsparse_create_csr_descr(matL, int(M,c_int64_t), int(M,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), c_loc(d_csr_val), &
rocsparse_indextype_i32, rocsparse_indextype_i32, &
rocsparse_index_base_zero, rocsparse_datatype_f64_r))
call rocsparseCheck(rocsparse_spmat_set_attribute(matL, rocsparse_spmat_fill_mode, c_loc(fill), int(4,c_size_t)))
call rocsparseCheck(rocsparse_spmat_set_attribute(matL, rocsparse_spmat_diag_type, c_loc(diag), int(4,c_size_t)))
! Dense-matrix descriptors: B is the rhs (X), C is the solution
call rocsparseCheck(rocsparse_create_dnmat_descr(matB, int(M,c_int64_t), int(nrhs,c_int64_t), int(M,c_int64_t), &
c_loc(d_X), rocsparse_datatype_f64_r, rocsparse_order_column))
call rocsparseCheck(rocsparse_create_dnmat_descr(matC, int(M,c_int64_t), int(nrhs,c_int64_t), int(M,c_int64_t), &
c_loc(d_C), rocsparse_datatype_f64_r, rocsparse_order_column))
! Stage 1: workspace size
call rocsparseCheck(rocsparse_spsm(handle, rocsparse_operation_none, rocsparse_operation_none, c_loc(alpha), &
matL, matB, matC, rocsparse_datatype_f64_r, rocsparse_spsm_alg_default, &
rocsparse_spsm_stage_buffer_size, buffer_size, c_null_ptr))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Stage 2: preprocess (analysis)
call rocsparseCheck(rocsparse_spsm(handle, rocsparse_operation_none, rocsparse_operation_none, c_loc(alpha), &
matL, matB, matC, rocsparse_datatype_f64_r, rocsparse_spsm_alg_default, &
rocsparse_spsm_stage_preprocess, buffer_size, d_buffer))
! Stage 3: solve
call rocsparseCheck(rocsparse_spsm(handle, rocsparse_operation_none, rocsparse_operation_none, c_loc(alpha), &
matL, matB, matC, rocsparse_datatype_f64_r, rocsparse_spsm_alg_default, &
rocsparse_spsm_stage_compute, buffer_size, d_buffer))
! Copy the recovered solution back
call hipCheck(hipMemcpy(h_C, d_C, hipMemcpyDeviceToHost))
! Verify C == Y
do j = 1,nrhs
do i = 1,M
error = abs(h_C(i,j) - h_Y(i,j)) / max(abs(h_Y(i,j)), 1.0_c_double)
if(error .gt. error_max) then
write(*,*) "FAILED! Error bigger than max! Error = ", error, " at (", i, ",", j, ")"
call exit
end if
end do
end do
! Clean up
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_csr_col_ind))
call hipCheck(hipFree(d_csr_val))
call hipCheck(hipFree(d_X))
call hipCheck(hipFree(d_C))
call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program dsptrsm
Sparse matrix arithmetic#
csrgemm multiplies two sparse matrices, C = alpha*A*B. Because the
sparsity pattern of C is not known in advance, the routine runs in two
passes: nnz first computes the number of nonzeros and the row pointers of
C, then the values pass fills the columns and values.
!!!!!!!!!!!!!/
! dcsrgemm example (double-precision sparse-matrix sparse-matrix multiply, C = alpha*A*B)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/extra.html
!
! Two-phase flow: csrgemm_buffer_size -> csrgemm_nnz (fills row_ptr_C and
! nnz_C) -> allocate col_ind_C/val_C -> csrgemm (computes C). Here B = A, so
! C = A*A, checked against the known product. The optional D term is unused
! (descr_D / arrays passed as null, beta = 0).
!
! NOTE: the mat-descr/array arguments are c_ptr-only, so device buffers are
! passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dcsrgemm
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
! Sparse A (3x3) in CSR (0-based): A = [[1,0,2],[0,3,0],[4,0,5]]
! C = A*A = [[9,0,12],[0,9,0],[24,0,33]]
integer(c_int), parameter :: M = 3, N = 3, K = 3, nnz_A = 5
integer(c_int) :: h_csr_row_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_csr_col_ind(5) = (/0, 2, 1, 0, 2/)
real(c_double) :: h_csr_val(5) = (/1, 2, 3, 4, 5/)
! Expected C
integer(c_int) :: h_exp_row_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_exp_col_ind(5) = (/0, 2, 1, 0, 2/)
real(c_double) :: h_exp_val(5) = (/9, 12, 9, 24, 33/)
real(c_double), target :: alpha = 1.0_c_double, beta = 0.0_c_double
integer(c_int) :: h_row_ptr_C(4)
integer(c_int) :: nnz_C
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:)
integer(c_int), pointer :: d_row_ptr_C(:), d_col_ind_C(:)
real(c_double), pointer :: d_val_C(:)
integer(c_int), pointer :: d_nnz_C
type(c_ptr) :: handle, descr_A, descr_B, descr_C, info_C, d_buffer
integer(c_size_t) :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dcsrgemm' (Fortran 2008 interfaces) - "
! Create handle, matrix descriptors and the csrgemm info object
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_create_mat_descr(descr_A))
call rocsparseCheck(rocsparse_create_mat_descr(descr_B))
call rocsparseCheck(rocsparse_create_mat_descr(descr_C))
call rocsparseCheck(rocsparse_create_mat_info(info_C))
! Copy A to device (B aliases A) and allocate row_ptr_C + nnz_C
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_csr_col_ind, source=h_csr_col_ind))
call hipCheck(hipMalloc(d_csr_val, source=h_csr_val))
call hipCheck(hipMalloc(d_row_ptr_C, mold=h_row_ptr_C))
call hipCheck(hipMalloc(d_nnz_C, source=0))
! Phase 0: workspace size
call rocsparseCheck(rocsparse_dcsrgemm_buffer_size(handle, rocsparse_operation_none, rocsparse_operation_none, &
M, N, K, alpha, descr_A, nnz_A, c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), &
descr_B, nnz_A, c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), &
beta, c_null_ptr, 0, c_null_ptr, c_null_ptr, info_C, buffer_size))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Phase 1: compute the sparsity of C (row_ptr_C + total nnz_C)
call rocsparseCheck(rocsparse_csrgemm_nnz(handle, rocsparse_operation_none, rocsparse_operation_none, M, N, K, &
descr_A, nnz_A, c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), &
descr_B, nnz_A, c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), &
c_null_ptr, 0, c_null_ptr, c_null_ptr, &
descr_C, c_loc(d_row_ptr_C), c_loc(d_nnz_C), info_C, d_buffer))
call hipCheck(hipMemcpy(nnz_C, d_nnz_C, hipMemcpyDeviceToHost))
! Allocate C column indices and values now that nnz_C is known
call hipCheck(hipMalloc(d_col_ind_C, dims=(/nnz_C/)))
call hipCheck(hipMalloc(d_val_C, dims=(/nnz_C/)))
! Phase 2: compute the values of C
call rocsparseCheck(rocsparse_dcsrgemm(handle, rocsparse_operation_none, rocsparse_operation_none, M, N, K, alpha, &
descr_A, nnz_A, c_loc(d_csr_val), c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), &
descr_B, nnz_A, c_loc(d_csr_val), c_loc(d_csr_row_ptr), c_loc(d_csr_col_ind), &
beta, c_null_ptr, 0, c_null_ptr, c_null_ptr, c_null_ptr, &
descr_C, c_loc(d_val_C), c_loc(d_row_ptr_C), c_loc(d_col_ind_C), info_C, d_buffer))
call hipCheck(hipDeviceSynchronize())
! Copy the C structure back to host
call hipCheck(hipMemcpy(h_row_ptr_C, d_row_ptr_C, hipMemcpyDeviceToHost))
! Verify nnz and row pointers
if(nnz_C /= 5) then
write(*,*) "FAILED! nnz_C = ", nnz_C, " expected 5"
call exit
end if
do i = 1,N+1
if(h_row_ptr_C(i) /= h_exp_row_ptr(i)) then
write(*,*) "FAILED! row_ptr_C(", i, ") = ", h_row_ptr_C(i), " expected ", h_exp_row_ptr(i)
call exit
end if
end do
! Verify column indices and values
block
integer(c_int) :: h_col_ind_C(nnz_C)
real(c_double) :: h_val_C(nnz_C)
call hipCheck(hipMemcpy(h_col_ind_C, d_col_ind_C, hipMemcpyDeviceToHost))
call hipCheck(hipMemcpy(h_val_C, d_val_C, hipMemcpyDeviceToHost))
do i = 1,nnz_C
if(h_col_ind_C(i) /= h_exp_col_ind(i)) then
write(*,*) "FAILED! col_ind_C(", i, ") = ", h_col_ind_C(i), " expected ", h_exp_col_ind(i)
call exit
end if
error = abs(h_val_C(i) - h_exp_val(i)) / max(abs(h_exp_val(i)), 1.0_c_double)
if(error .gt. error_max) then
write(*,*) "FAILED! val_C(", i, ") = ", h_val_C(i), " expected ", h_exp_val(i)
call exit
end if
end do
end block
! Clean up
call hipCheck(hipFree(d_col_ind_C))
call hipCheck(hipFree(d_val_C))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_csr_col_ind))
call hipCheck(hipFree(d_csr_val))
call hipCheck(hipFree(d_row_ptr_C))
call hipCheck(hipFree(d_nnz_C))
call hipCheck(hipFree(d_buffer))
call rocsparseCheck(rocsparse_destroy_mat_info(info_C))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr_A))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr_B))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr_C))
call rocsparseCheck(rocsparse_destroy_handle(handle))
write(*,*) "PASSED!"
end program dcsrgemm
csrgeam adds two sparse matrices, C = alpha*A + beta*B, with the same
two-pass structure.
!!!!!!!!!!!!!/
! dcsrgeam example (double-precision sparse matrix addition, C = alpha*A + beta*B)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/extra.html
!
! Two-phase flow: csrgeam_nnz (fills row_ptr_C and nnz_C) -> allocate
! col_ind_C/val_C -> csrgeam (computes C). csrgeam needs no workspace buffer.
! Here C = A + B (alpha = beta = 1), checked against the known sum.
!
! NOTE: the mat-descr/array arguments are c_ptr-only, so device buffers are
! passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dcsrgeam
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
! A (3x3) CSR (0-based): A = [[1,0,2],[0,3,0],[4,0,5]]
! B (3x3) CSR (0-based): B = diag(10,20,30)
! C = A + B = [[11,0,2],[0,23,0],[4,0,35]] (diagonals merge)
integer(c_int), parameter :: M = 3, N = 3, nnz_A = 5, nnz_B = 3
integer(c_int) :: h_row_ptr_A(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_col_ind_A(5) = (/0, 2, 1, 0, 2/)
real(c_double) :: h_val_A(5) = (/1, 2, 3, 4, 5/)
integer(c_int) :: h_row_ptr_B(4) = (/0, 1, 2, 3/)
integer(c_int) :: h_col_ind_B(3) = (/0, 1, 2/)
real(c_double) :: h_val_B(3) = (/10, 20, 30/)
! Expected C
integer(c_int) :: h_exp_row_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_exp_col_ind(5) = (/0, 2, 1, 0, 2/)
real(c_double) :: h_exp_val(5) = (/11, 2, 23, 4, 35/)
real(c_double), target :: alpha = 1.0_c_double, beta = 1.0_c_double
integer(c_int) :: h_row_ptr_C(4)
integer(c_int) :: nnz_C
integer(c_int), pointer :: d_row_ptr_A(:), d_col_ind_A(:)
real(c_double), pointer :: d_val_A(:)
integer(c_int), pointer :: d_row_ptr_B(:), d_col_ind_B(:)
real(c_double), pointer :: d_val_B(:)
integer(c_int), pointer :: d_row_ptr_C(:), d_col_ind_C(:)
real(c_double), pointer :: d_val_C(:)
integer(c_int), pointer :: d_nnz_C
type(c_ptr) :: handle, descr_A, descr_B, descr_C
real(c_double) :: error
real(c_double), parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dcsrgeam' (Fortran 2008 interfaces) - "
! Create handle and matrix descriptors
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_create_mat_descr(descr_A))
call rocsparseCheck(rocsparse_create_mat_descr(descr_B))
call rocsparseCheck(rocsparse_create_mat_descr(descr_C))
! Copy A and B to device; allocate row_ptr_C + nnz_C
call hipCheck(hipMalloc(d_row_ptr_A, source=h_row_ptr_A))
call hipCheck(hipMalloc(d_col_ind_A, source=h_col_ind_A))
call hipCheck(hipMalloc(d_val_A, source=h_val_A))
call hipCheck(hipMalloc(d_row_ptr_B, source=h_row_ptr_B))
call hipCheck(hipMalloc(d_col_ind_B, source=h_col_ind_B))
call hipCheck(hipMalloc(d_val_B, source=h_val_B))
call hipCheck(hipMalloc(d_row_ptr_C, mold=h_row_ptr_C))
call hipCheck(hipMalloc(d_nnz_C, source=0))
! Phase 1: compute the sparsity of C (row_ptr_C + total nnz_C)
call rocsparseCheck(rocsparse_csrgeam_nnz(handle, M, N, &
descr_A, nnz_A, c_loc(d_row_ptr_A), c_loc(d_col_ind_A), &
descr_B, nnz_B, c_loc(d_row_ptr_B), c_loc(d_col_ind_B), &
descr_C, c_loc(d_row_ptr_C), c_loc(d_nnz_C)))
call hipCheck(hipMemcpy(nnz_C, d_nnz_C, hipMemcpyDeviceToHost))
! Allocate C column indices and values now that nnz_C is known
call hipCheck(hipMalloc(d_col_ind_C, dims=(/nnz_C/)))
call hipCheck(hipMalloc(d_val_C, dims=(/nnz_C/)))
! Phase 2: compute the values of C
call rocsparseCheck(rocsparse_dcsrgeam(handle, M, N, alpha, &
descr_A, nnz_A, c_loc(d_val_A), c_loc(d_row_ptr_A), c_loc(d_col_ind_A), &
beta, descr_B, nnz_B, c_loc(d_val_B), c_loc(d_row_ptr_B), c_loc(d_col_ind_B), &
descr_C, c_loc(d_val_C), c_loc(d_row_ptr_C), c_loc(d_col_ind_C)))
call hipCheck(hipDeviceSynchronize())
! Copy the C structure back to host
call hipCheck(hipMemcpy(h_row_ptr_C, d_row_ptr_C, hipMemcpyDeviceToHost))
! Verify nnz and row pointers
if(nnz_C /= 5) then
write(*,*) "FAILED! nnz_C = ", nnz_C, " expected 5"
call exit
end if
do i = 1,N+1
if(h_row_ptr_C(i) /= h_exp_row_ptr(i)) then
write(*,*) "FAILED! row_ptr_C(", i, ") = ", h_row_ptr_C(i), " expected ", h_exp_row_ptr(i)
call exit
end if
end do
! Verify column indices and values
block
integer(c_int) :: h_col_ind_C(nnz_C)
real(c_double) :: h_val_C(nnz_C)
call hipCheck(hipMemcpy(h_col_ind_C, d_col_ind_C, hipMemcpyDeviceToHost))
call hipCheck(hipMemcpy(h_val_C, d_val_C, hipMemcpyDeviceToHost))
do i = 1,nnz_C
if(h_col_ind_C(i) /= h_exp_col_ind(i)) then
write(*,*) "FAILED! col_ind_C(", i, ") = ", h_col_ind_C(i), " expected ", h_exp_col_ind(i)
call exit
end if
error = abs(h_val_C(i) - h_exp_val(i)) / max(abs(h_exp_val(i)), 1.0_c_double)
if(error .gt. error_max) then
write(*,*) "FAILED! val_C(", i, ") = ", h_val_C(i), " expected ", h_exp_val(i)
call exit
end if
end do
end block
! Clean up
call hipCheck(hipFree(d_col_ind_C))
call hipCheck(hipFree(d_val_C))
call hipCheck(hipFree(d_row_ptr_A))
call hipCheck(hipFree(d_col_ind_A))
call hipCheck(hipFree(d_val_A))
call hipCheck(hipFree(d_row_ptr_B))
call hipCheck(hipFree(d_col_ind_B))
call hipCheck(hipFree(d_val_B))
call hipCheck(hipFree(d_row_ptr_C))
call hipCheck(hipFree(d_nnz_C))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr_A))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr_B))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr_C))
call rocsparseCheck(rocsparse_destroy_handle(handle))
write(*,*) "PASSED!"
end program dcsrgeam
Block-sparse matrix-vector products#
bsrmv multiplies a matrix stored in BSR (block sparse row) format by a dense
vector. BSR groups the nonzeros into fixed-size dense blocks, which suits
matrices with a natural block structure. The routine takes a matrix descriptor
and a matrix-info handle.
!!!!!!!!!!!!!!
! rocsparse sbsrmv example (block-sparse matrix-vector multiply, single)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Computes y = alpha * A * x + beta * y for a BSR matrix A. Here A is
! block-diagonal with mb=nb=2 blocks of block_dim=2 (a 4x4 dense equivalent):
! A = [ B1 0 ] B1 = [1 2] B2 = [5 6]
! [ 0 B2 ] [3 4] [7 8]
! Blocks are stored row-major (rocsparse_direction_row). The result is checked
! against the dense reference y = A_dense * x.
!!!!!!!!!!!!!!
!
program dbsrmv
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: mb = 2, nb = 2, nnzb = 2, block_dim = 2
integer(c_int), parameter :: mdim = mb * block_dim
real(c_double) :: hVal(8) = (/1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0, 7.0d0, 8.0d0/)
integer(c_int) :: hRowPtr(3) = (/0, 1, 2/)
integer(c_int) :: hColInd(2) = (/0, 1/)
real(c_double) :: hX(4) = (/1.0d0, 2.0d0, 3.0d0, 4.0d0/)
real(c_double) :: hY(4) = (/0.0d0, 0.0d0, 0.0d0, 0.0d0/)
real(c_double) :: hRef(4)
real(c_double) :: alpha = 1.0d0, beta = 0.0d0
type(c_ptr) :: handle = c_null_ptr
type(c_ptr) :: descr = c_null_ptr, info = c_null_ptr
real(c_double), pointer :: dVal(:), dX(:), dY(:)
integer(c_int), pointer :: dRowPtr(:), dColInd(:)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dbsrmv' (Fortran 2008 interfaces) - "
hRef(1) = 1.0d0*hX(1) + 2.0d0*hX(2)
hRef(2) = 3.0d0*hX(1) + 4.0d0*hX(2)
hRef(3) = 5.0d0*hX(3) + 6.0d0*hX(4)
hRef(4) = 7.0d0*hX(3) + 8.0d0*hX(4)
call hipCheck(hipMalloc(dVal, source=hVal))
call hipCheck(hipMalloc(dRowPtr, source=hRowPtr))
call hipCheck(hipMalloc(dColInd, source=hColInd))
call hipCheck(hipMalloc(dX, source=hX))
call hipCheck(hipMalloc(dY, source=hY))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_create_mat_descr(descr))
call rocsparseCheck(rocsparse_create_mat_info(info))
call rocsparseCheck(rocsparse_dbsrmv(handle, rocsparse_direction_row, rocsparse_operation_none, &
mb, nb, nnzb, alpha, descr, dVal, dRowPtr, dColInd, block_dim, info, dX, beta, dY))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hY, dY, hipMemcpyDeviceToHost))
do i = 1, mdim
if (abs(hY(i) - hRef(i)) > 1.0d-11) then
write(*,*) "FAILED! y(", i, ") = ", hY(i), " expected ", hRef(i); call exit(1)
end if
end do
call rocsparseCheck(rocsparse_destroy_mat_info(info))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr))
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(dVal)); call hipCheck(hipFree(dRowPtr)); call hipCheck(hipFree(dColInd))
call hipCheck(hipFree(dX)); call hipCheck(hipFree(dY))
write(*,*) "PASSED!"
end program dbsrmv
gebsrmv is the general variant, allowing rectangular blocks with separate
row and column block dimensions.
!!!!!!!!!!!!!!
! rocsparse sgebsrmv example (general block-sparse matrix-vector multiply, single)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Computes y = alpha * A * x + beta * y for a BSR matrix A. Here A is
! block-diagonal with mb=nb=2 blocks of block_dim=2 (a 4x4 dense equivalent):
! A = [ B1 0 ] B1 = [1 2] B2 = [5 6]
! [ 0 B2 ] [3 4] [7 8]
! Blocks are stored row-major (rocsparse_direction_row). The result is checked
! against the dense reference y = A_dense * x.
!!!!!!!!!!!!!!
!
program dgebsrmv
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: mb = 2, nb = 2, nnzb = 2, block_dim = 2
integer(c_int), parameter :: mdim = mb * block_dim
real(c_double) :: hVal(8) = (/1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0, 7.0d0, 8.0d0/)
integer(c_int) :: hRowPtr(3) = (/0, 1, 2/)
integer(c_int) :: hColInd(2) = (/0, 1/)
real(c_double) :: hX(4) = (/1.0d0, 2.0d0, 3.0d0, 4.0d0/)
real(c_double) :: hY(4) = (/0.0d0, 0.0d0, 0.0d0, 0.0d0/)
real(c_double) :: hRef(4)
real(c_double) :: alpha = 1.0d0, beta = 0.0d0
type(c_ptr) :: handle = c_null_ptr
type(c_ptr) :: descr = c_null_ptr
real(c_double), pointer :: dVal(:), dX(:), dY(:)
integer(c_int), pointer :: dRowPtr(:), dColInd(:)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dgebsrmv' (Fortran 2008 interfaces) - "
hRef(1) = 1.0d0*hX(1) + 2.0d0*hX(2)
hRef(2) = 3.0d0*hX(1) + 4.0d0*hX(2)
hRef(3) = 5.0d0*hX(3) + 6.0d0*hX(4)
hRef(4) = 7.0d0*hX(3) + 8.0d0*hX(4)
call hipCheck(hipMalloc(dVal, source=hVal))
call hipCheck(hipMalloc(dRowPtr, source=hRowPtr))
call hipCheck(hipMalloc(dColInd, source=hColInd))
call hipCheck(hipMalloc(dX, source=hX))
call hipCheck(hipMalloc(dY, source=hY))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_create_mat_descr(descr))
call rocsparseCheck(rocsparse_dgebsrmv(handle, rocsparse_direction_row, rocsparse_operation_none, &
mb, nb, nnzb, alpha, descr, dVal, dRowPtr, dColInd, block_dim, block_dim, dX, beta, dY))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hY, dY, hipMemcpyDeviceToHost))
do i = 1, mdim
if (abs(hY(i) - hRef(i)) > 1.0d-11) then
write(*,*) "FAILED! y(", i, ") = ", hY(i), " expected ", hRef(i); call exit(1)
end if
end do
call rocsparseCheck(rocsparse_destroy_mat_descr(descr))
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(dVal)); call hipCheck(hipFree(dRowPtr)); call hipCheck(hipFree(dColInd))
call hipCheck(hipFree(dX)); call hipCheck(hipFree(dY))
write(*,*) "PASSED!"
end program dgebsrmv
Incomplete factorization preconditioners#
Incomplete factorizations produce approximate factors that keep the sparsity of
the input and are used as preconditioners. Each runs an analysis stage before
the compute stage. csrilu0 computes an incomplete LU factorization with zero
fill-in; on a matrix with no fill-in (such as a tridiagonal one) it reproduces
the exact LU, which the example checks.
!!!!!!!!!!!!!!
! rocsparse scsrilu0 example (incomplete LU, single)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Computes the ILU(0) factorization of a sparse matrix in place
! (buffer_size -> analysis -> compute, with a mat descriptor and a mat info
! object). For a tridiagonal matrix there is no fill-in, so ILU(0) equals the
! exact LU factorization and the overwritten CSR values are checked against the
! hand-computed factors (unit lower L, upper U packed into one array).
!
! A = [ 4 1 0 ] LU (in place) = [ 4 1 0 ]
! [ 1 4 1 ] [ 1/4 15/4 1 ]
! [ 0 1 4 ] [ 0 4/15 56/15]
!!!!!!!!!!!!!!
!
program dcsrilu0
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: m = 3, nnz = 7
integer(c_int) :: hRowPtr(4) = (/0, 2, 5, 7/)
integer(c_int) :: hColInd(7) = (/0, 1, 0, 1, 2, 1, 2/)
real(c_double) :: hVal(7) = (/4.0d0, 1.0d0, 1.0d0, 4.0d0, 1.0d0, 1.0d0, 4.0d0/)
real(c_double) :: hOut(7)
real(c_double) :: hExp(7)
type(c_ptr) :: handle = c_null_ptr
type(c_ptr) :: descr = c_null_ptr, info = c_null_ptr
integer(c_int), pointer :: dRowPtr(:), dColInd(:)
real(c_double), pointer :: dVal(:)
type(c_ptr) :: dBuf
integer(c_size_t) :: bufSize
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dcsrilu0' (Fortran 2008 interfaces) - "
hExp(1) = 4.0d0
hExp(2) = 1.0d0
hExp(3) = 1.0d0/4.0d0
hExp(4) = 4.0d0 - (1.0d0/4.0d0)*1.0d0
hExp(5) = 1.0d0
hExp(6) = 1.0d0/hExp(4)
hExp(7) = 4.0d0 - hExp(6)*1.0d0
call hipCheck(hipMalloc(dRowPtr, source=hRowPtr))
call hipCheck(hipMalloc(dColInd, source=hColInd))
call hipCheck(hipMalloc(dVal, source=hVal))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_create_mat_descr(descr))
call rocsparseCheck(rocsparse_create_mat_info(info))
call rocsparseCheck(rocsparse_dcsrilu0_buffer_size(handle, m, nnz, descr, dVal, dRowPtr, &
dColInd, info, bufSize))
call hipCheck(hipMalloc(dBuf, max(bufSize, 1_c_size_t)))
call rocsparseCheck(rocsparse_dcsrilu0_analysis(handle, m, nnz, descr, dVal, dRowPtr, dColInd, &
info, rocsparse_analysis_policy_reuse, rocsparse_solve_policy_auto, dBuf))
call rocsparseCheck(rocsparse_dcsrilu0(handle, m, nnz, descr, dVal, dRowPtr, dColInd, &
info, rocsparse_solve_policy_auto, dBuf))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hOut, dVal, hipMemcpyDeviceToHost))
do i = 1, nnz
if (abs(hOut(i) - hExp(i)) > 1.0d-11) then
write(*,*) "FAILED! val(", i, ") = ", hOut(i), " expected ", hExp(i); call exit(1)
end if
end do
call rocsparseCheck(rocsparse_destroy_mat_info(info))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr))
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(dRowPtr)); call hipCheck(hipFree(dColInd)); call hipCheck(hipFree(dVal))
call hipCheck(hipFree(dBuf))
write(*,*) "PASSED!"
end program dcsrilu0
csric0 is the incomplete Cholesky counterpart for a symmetric positive
definite matrix.
!!!!!!!!!!!!!!
! rocsparse scsric0 example (incomplete Cholesky, single)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Computes the IC(0) factorization of an SPD sparse matrix in place
! (buffer_size -> analysis -> compute, with a mat descriptor and a mat info
! object). For an SPD tridiagonal matrix there is no fill-in, so IC(0) equals the
! exact Cholesky factor L (A = L*L^T). csric0 overwrites the lower-triangular
! part (including the diagonal) with L; the strict upper part is left unchanged.
! The lower/diagonal CSR entries are checked against the hand-computed L.
!
! A = [ 4 1 0 ] L = [ 2 0 0 ]
! [ 1 4 1 ] [ 1/2 sqrt(15)/2 0 ]
! [ 0 1 4 ] [ 0 2/sqrt(15) l33 ]
!!!!!!!!!!!!!!
!
program dcsric0
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: m = 3, nnz = 7
integer(c_int) :: hRowPtr(4) = (/0, 2, 5, 7/)
integer(c_int) :: hColInd(7) = (/0, 1, 0, 1, 2, 1, 2/)
real(c_double) :: hVal(7) = (/4.0d0, 1.0d0, 1.0d0, 4.0d0, 1.0d0, 1.0d0, 4.0d0/)
real(c_double) :: hOut(7)
real(c_double) :: L11, L21, L22, L32, L33
type(c_ptr) :: handle = c_null_ptr
type(c_ptr) :: descr = c_null_ptr, info = c_null_ptr
integer(c_int), pointer :: dRowPtr(:), dColInd(:)
real(c_double), pointer :: dVal(:)
type(c_ptr) :: dBuf
integer(c_size_t) :: bufSize
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dcsric0' (Fortran 2008 interfaces) - "
L11 = sqrt(4.0d0)
L21 = 1.0d0 / L11
L22 = sqrt(4.0d0 - L21*L21)
L32 = 1.0d0 / L22
L33 = sqrt(4.0d0 - L32*L32)
call hipCheck(hipMalloc(dRowPtr, source=hRowPtr))
call hipCheck(hipMalloc(dColInd, source=hColInd))
call hipCheck(hipMalloc(dVal, source=hVal))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_create_mat_descr(descr))
call rocsparseCheck(rocsparse_create_mat_info(info))
call rocsparseCheck(rocsparse_dcsric0_buffer_size(handle, m, nnz, descr, dVal, dRowPtr, &
dColInd, info, bufSize))
call hipCheck(hipMalloc(dBuf, max(bufSize, 1_c_size_t)))
call rocsparseCheck(rocsparse_dcsric0_analysis(handle, m, nnz, descr, dVal, dRowPtr, dColInd, &
info, rocsparse_analysis_policy_reuse, rocsparse_solve_policy_auto, dBuf))
call rocsparseCheck(rocsparse_dcsric0(handle, m, nnz, descr, dVal, dRowPtr, dColInd, &
info, rocsparse_solve_policy_auto, dBuf))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hOut, dVal, hipMemcpyDeviceToHost))
if (abs(hOut(1) - L11) > 1.0d-11 .or. abs(hOut(3) - L21) > 1.0d-11 .or. &
abs(hOut(4) - L22) > 1.0d-11 .or. abs(hOut(6) - L32) > 1.0d-11 .or. &
abs(hOut(7) - L33) > 1.0d-11) then
write(*,*) "FAILED! L = ", hOut(1), hOut(3), hOut(4), hOut(6), hOut(7), &
" expected ", L11, L21, L22, L32, L33
call exit(1)
end if
call rocsparseCheck(rocsparse_destroy_mat_info(info))
call rocsparseCheck(rocsparse_destroy_mat_descr(descr))
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(dRowPtr)); call hipCheck(hipFree(dColInd)); call hipCheck(hipFree(dVal))
call hipCheck(hipFree(dBuf))
write(*,*) "PASSED!"
end program dcsric0
spildlt0 computes an incomplete LDLH factorization through the
generic descriptor API, with descriptor-create, set-input, analysis and compute
stages, and a get-output query for the singularity status.
!!!!!!!!!!!!!!
! rocsparse spildlt0 example (incomplete LDL^H factorization, level 0)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Exercises the generic staged SpILDLT0 preconditioner API: create the config
! descriptor, set inputs (algorithm, compute datatype, analysis policy), then run
! the analysis and compute stages (each with its own workspace queried by
! spildlt0_buffer_size). The matrix to factorize (A) is also passed as the
! preconditioner (P), as required. Success is verified by the post-compute
! singularity output being rocsparse_singularity_none.
!
! The descriptor/spmat/config arguments are c_ptr-only, so device buffers are
! passed via c_loc and scalar host inputs via c_loc(...).
!!!!!!!!!!!!!!
!
program spildlt0
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer(c_int), parameter :: m = 4, nnz = 7
! SPD lower-CSR (from the rocSPARSE example): factorizable without singularity.
integer(c_int), target :: h_row_ptr(5) = (/0, 1, 3, 5, 7/)
integer(c_int), target :: h_col_ind(7) = (/0, 0, 1, 1, 2, 2, 3/)
real(c_double), target :: h_val(7) = (/4.0d0, 2.0d0, 8.0d0, 1.0d0, 8.0d0, 2.0d0, 4.0d0/)
integer(c_int), pointer :: d_row_ptr(:), d_col_ind(:)
real(c_double), pointer :: d_val(:)
type(c_ptr) :: handle, matA, descr, dBuf
integer(kind(rocsparse_spildlt0_alg_default)), target :: alg = rocsparse_spildlt0_alg_default
integer(kind(rocsparse_analysis_policy_reuse)), target :: apol = rocsparse_analysis_policy_reuse
integer(kind(rocsparse_datatype_f64_r)), target :: cdt = rocsparse_datatype_f64_r
integer(kind(rocsparse_singularity_none)), target :: sing
integer(c_size_t), target :: bufSize
write(*,"(a)",advance="no") "-- Running test 'rocsparse_spildlt0' (Fortran 2008 interfaces) - "
call hipCheck(hipMalloc(d_row_ptr, source=h_row_ptr))
call hipCheck(hipMalloc(d_col_ind, source=h_col_ind))
call hipCheck(hipMalloc(d_val, source=h_val))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_create_csr_descr(matA, int(m,c_int64_t), int(m,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_row_ptr), c_loc(d_col_ind), c_loc(d_val), &
rocsparse_indextype_i32, rocsparse_indextype_i32, rocsparse_index_base_zero, rocsparse_datatype_f64_r))
call rocsparseCheck(rocsparse_spildlt0_descr_create(handle, descr, c_null_ptr))
call rocsparseCheck(rocsparse_spildlt0_set_input(handle, descr, rocsparse_spildlt0_input_alg, &
c_loc(alg), int(c_sizeof(alg),c_size_t), c_null_ptr))
call rocsparseCheck(rocsparse_spildlt0_set_input(handle, descr, rocsparse_spildlt0_input_compute_datatype, &
c_loc(cdt), int(c_sizeof(cdt),c_size_t), c_null_ptr))
call rocsparseCheck(rocsparse_spildlt0_set_input(handle, descr, rocsparse_spildlt0_input_analysis_policy, &
c_loc(apol), int(c_sizeof(apol),c_size_t), c_null_ptr))
! Analysis stage.
call rocsparseCheck(rocsparse_spildlt0_buffer_size(handle, descr, matA, matA, &
rocsparse_spildlt0_stage_analysis, c_loc(bufSize), c_null_ptr))
call hipCheck(hipMalloc(dBuf, max(bufSize, 1_c_size_t)))
call rocsparseCheck(rocsparse_spildlt0(handle, descr, matA, matA, &
rocsparse_spildlt0_stage_analysis, bufSize, dBuf, c_null_ptr))
call hipCheck(hipFree(dBuf))
! Compute stage.
call rocsparseCheck(rocsparse_spildlt0_buffer_size(handle, descr, matA, matA, &
rocsparse_spildlt0_stage_compute, c_loc(bufSize), c_null_ptr))
call hipCheck(hipMalloc(dBuf, max(bufSize, 1_c_size_t)))
call rocsparseCheck(rocsparse_spildlt0(handle, descr, matA, matA, &
rocsparse_spildlt0_stage_compute, bufSize, dBuf, c_null_ptr))
call hipCheck(hipDeviceSynchronize())
call rocsparseCheck(rocsparse_set_pointer_mode(handle, rocsparse_pointer_mode_host))
call rocsparseCheck(rocsparse_spildlt0_get_output(handle, descr, rocsparse_spildlt0_output_singularity, &
c_loc(sing), int(c_sizeof(sing),c_size_t), c_null_ptr))
if (sing /= rocsparse_singularity_none) then
write(*,*) "FAILED! singularity = ", sing, " (expected none=0)"; call exit(1)
end if
call rocsparseCheck(rocsparse_spildlt0_descr_destroy(handle, descr, c_null_ptr))
call rocsparseCheck(rocsparse_destroy_spmat_descr(matA))
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_row_ptr)); call hipCheck(hipFree(d_col_ind)); call hipCheck(hipFree(d_val))
call hipCheck(hipFree(dBuf))
write(*,*) "PASSED!"
end program spildlt0
Tridiagonal and pentadiagonal solvers#
gtsv solves a tridiagonal system given its three diagonals. It is a direct
banded solver rather than an iterative one.
!!!!!!!!!!!!!!
! sgtsv example (single-precision tridiagonal solve, rocSPARSE)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/precond.html
!
! Solves A x = b for a tridiagonal system. Here A is the identity tridiagonal
! (dl = du = 0, d = 1), so the exact solution of A x = b is x = b. dl/d/du/B are
! passed as Fortran arrays, exercising the generic array form of the diagonal
! arguments (they used to be declared type(c_ptr), SWDEV-485451).
!!!!!!!!!!!!!!
!
program rocsparse_sgtsv_test
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer, parameter :: m = 512 ! tridiagonal system size
integer, parameter :: n = 1 ! single right-hand side
! Identity tridiagonal (dl = du = 0, d = 1) so the solution of A x = b is x = b.
real(c_float), allocatable, dimension(:) :: hdl, hd, hdu, hB
real(c_float), pointer, dimension(:) :: ddl => null(), dd => null(), &
ddu => null(), dB => null()
type(c_ptr) :: handle = c_null_ptr
type(c_ptr) :: dbuffer = c_null_ptr
integer(c_size_t) :: buffer_size
integer :: i
real(c_float) :: error
real(c_float), parameter :: error_max = 10*epsilon(error)
write(*,"(a)",advance="no") "-- Running test 'SGTSV' (Fortran 2008 interfaces) - "
call rocsparseCheck(rocsparse_create_handle(handle))
allocate(hdl(m), hd(m), hdu(m), hB(m))
hdl(:) = 0.0
hd(:) = 1.0
hdu(:) = 0.0
do i = 1, m
hB(i) = real(i) ! b(i) = i -> exact solution x(i) = i (identity system)
end do
call hipCheck(hipMalloc(ddl, source=hdl))
call hipCheck(hipMalloc(dd, source=hd))
call hipCheck(hipMalloc(ddu, source=hdu))
call hipCheck(hipMalloc(dB, source=hB))
! Query the temporary buffer size. dl/d/du/B are passed as Fortran arrays —
! this is the generic array form that used to fail to compile because the
! diagonals were declared type(c_ptr) (SWDEV-485451).
call rocsparseCheck(rocsparse_sgtsv_buffer_size(handle, m, n, ddl, dd, ddu, dB, m, buffer_size))
call hipCheck(hipMalloc(dbuffer, buffer_size))
! Tridiagonal solve A x = b, in place; dl/d/du/B passed as Fortran arrays.
call rocsparseCheck(rocsparse_sgtsv(handle, m, n, ddl, dd, ddu, dB, m, dbuffer))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hB, dB, hipMemcpyDeviceToHost))
do i = 1, m
error = abs((real(i) - hB(i))/real(i))
if( error > error_max )then
write(*,*) "FAILED! Error bigger than max! Error = ", error, " hB(", i, ") = ", hB(i)
call exit(1)
end if
end do
call hipCheck(hipFree(ddl))
call hipCheck(hipFree(dd))
call hipCheck(hipFree(ddu))
call hipCheck(hipFree(dB))
call hipCheck(hipFree(dbuffer))
call rocsparseCheck(rocsparse_destroy_handle(handle))
deallocate(hdl, hd, hdu, hB)
write(*,*) "PASSED!"
end program rocsparse_sgtsv_test
gpsv_interleaved_batch solves a batch of pentadiagonal systems whose data is
interleaved across the batch, a layout that lets the GPU coalesce memory access
across the independent systems.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Copyright (c) 2026 Advanced Micro Devices, Inc.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program rocsparse_zgpsv_interleaved_batch_test
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
! Pentadiagonal system size and (single) batch layout
integer(c_int), parameter :: M = 6
integer(c_int), parameter :: batch_count = 1
integer(c_int), parameter :: batch_stride = 1
! Five diagonals of the pentadiagonal system (interleaved, batch_count = 1)
! ds : lower diagonal at distance 2 (first two entries are zero)
! dl : lower diagonal (first entry is zero)
! d : main diagonal
! du : upper diagonal (last entry is zero)
! dw : upper diagonal at distance 2 (last two entries are zero)
complex(8) :: h_ds(M), h_dl(M), h_d(M), h_du(M), h_dw(M)
complex(8) :: h_x(M), h_x_exact(M)
complex(8) :: A(M,M)
complex(8), pointer :: d_ds(:), d_dl(:), d_d(:), d_du(:), d_dw(:), d_x(:)
type(c_ptr) :: d_buffer
integer(c_size_t) :: buffer_size
integer :: i
real(8) :: err
type(c_ptr) :: handle
write(*,"(a)",advance="no") "-- Running test 'zgpsv_interleaved_batch' (Fortran 2008 interfaces) - "
! Build a diagonally dominant pentadiagonal matrix
h_ds = (0.0d0, 0.0d0)
h_dl = (0.0d0, 0.0d0)
h_d = (0.0d0, 0.0d0)
h_du = (0.0d0, 0.0d0)
h_dw = (0.0d0, 0.0d0)
do i = 1, M
h_d(i) = (10.0d0, 1.0d0)
end do
do i = 2, M
h_dl(i) = (-1.0d0, 0.0d0)
end do
do i = 1, M - 1
h_du(i) = (-1.0d0, 0.0d0)
end do
do i = 3, M
h_ds(i) = (-2.0d0, 1.0d0)
end do
do i = 1, M - 2
h_dw(i) = (-2.0d0, -1.0d0)
end do
! Known solution and matching right-hand side (rhs = A * x_exact)
do i = 1, M
h_x_exact(i) = cmplx(real(i, 8), -real(i, 8), 8)
end do
A = (0.0d0, 0.0d0)
do i = 1, M
A(i,i) = h_d(i)
if (i >= 2) A(i,i-1) = h_dl(i)
if (i >= 3) A(i,i-2) = h_ds(i)
if (i <= M - 1) A(i,i+1) = h_du(i)
if (i <= M - 2) A(i,i+2) = h_dw(i)
end do
h_x = matmul(A, h_x_exact)
! Allocate device memory and copy host data to device
call hipCheck(hipMalloc(d_ds, source=h_ds))
call hipCheck(hipMalloc(d_dl, source=h_dl))
call hipCheck(hipMalloc(d_d, source=h_d))
call hipCheck(hipMalloc(d_du, source=h_du))
call hipCheck(hipMalloc(d_dw, source=h_dw))
call hipCheck(hipMalloc(d_x, source=h_x))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
! Query the required temporary buffer size
call rocsparseCheck(rocsparse_zgpsv_interleaved_batch_buffer_size(handle, &
rocsparse_gpsv_interleaved_alg_qr, &
M, &
d_ds, d_dl, d_d, &
d_du, d_dw, &
d_x, &
batch_count, &
batch_stride, &
buffer_size))
call hipCheck(hipMalloc(d_buffer, buffer_size))
! Solve the batched pentadiagonal system (x is overwritten with the solution)
call rocsparseCheck(rocsparse_zgpsv_interleaved_batch(handle, &
rocsparse_gpsv_interleaved_alg_qr, &
M, &
d_ds, d_dl, d_d, &
d_du, d_dw, &
d_x, &
batch_count, &
batch_stride, &
d_buffer))
! Copy the solution back to the host
call hipCheck(hipMemcpy(h_x, d_x, hipMemcpyDeviceToHost))
! Verification
err = 0.0d0
do i = 1, M
err = max(err, abs(h_x(i) - h_x_exact(i)))
end do
if (err > 1.0d-8) then
write(*,*) 'FAILED! max error =', err
call exit(1)
end if
! Clear rocSPARSE
call rocsparseCheck(rocsparse_destroy_handle(handle))
! Clear device memory
call hipCheck(hipFree(d_ds))
call hipCheck(hipFree(d_dl))
call hipCheck(hipFree(d_d))
call hipCheck(hipFree(d_du))
call hipCheck(hipFree(d_dw))
call hipCheck(hipFree(d_x))
call hipCheck(hipFree(d_buffer))
! Print success
write(*,*) 'PASSED!'
end program rocsparse_zgpsv_interleaved_batch_test
Sparse vector operations#
gthr gathers the entries of a dense vector y at a set of indices into a
compact sparse vector x_val, and sctr scatters a sparse vector back into
a dense one. They are the pack/unpack pair for the sparse-vector format.
!!!!!!!!!!!!!!
! rocsparse sgthr example (gather y[x_ind] -> x_val, single)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Gathers the entries of a dense vector y at the sparse index set x_ind into the
! packed vector x_val, then checks the gathered values.
!!!!!!!!!!!!!!
!
program dgthr
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: n = 5, nnz = 3
real(c_double) :: hY(5) = (/10.0d0, 20.0d0, 30.0d0, 40.0d0, 50.0d0/)
integer(c_int) :: hXind(3) = (/0, 2, 4/)
real(c_double) :: hXval(3)
real(c_double) :: hExp(3) = (/10.0d0, 30.0d0, 50.0d0/)
type(c_ptr) :: handle = c_null_ptr
real(c_double), pointer :: dY(:), dXval(:)
integer(c_int), pointer :: dXind(:)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dgthr' (Fortran 2008 interfaces) - "
call hipCheck(hipMalloc(dY, source=hY))
call hipCheck(hipMalloc(dXind, source=hXind))
call hipCheck(hipMalloc(dXval, mold=hXval))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_dgthr(handle, nnz, dY, dXval, dXind, rocsparse_index_base_zero))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hXval, dXval, hipMemcpyDeviceToHost))
do i = 1, nnz
if (abs(hXval(i) - hExp(i)) > 1.0d-12) then
write(*,*) "FAILED! x_val(", i, ") = ", hXval(i), " expected ", hExp(i); call exit(1)
end if
end do
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(dY)); call hipCheck(hipFree(dXind)); call hipCheck(hipFree(dXval))
write(*,*) "PASSED!"
end program dgthr
!!!!!!!!!!!!!!
! rocsparse ssctr example (scatter x_val -> y[x_ind], single)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Scatters the packed vector x_val into the dense vector y at the sparse index
! set x_ind, then checks the resulting dense vector.
!!!!!!!!!!!!!!
!
program dsctr
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: n = 5, nnz = 3
real(c_double) :: hXval(3) = (/100.0d0, 200.0d0, 300.0d0/)
integer(c_int) :: hXind(3) = (/0, 2, 4/)
real(c_double) :: hY(5) = (/0.0d0, 0.0d0, 0.0d0, 0.0d0, 0.0d0/)
real(c_double) :: hExp(5) = (/100.0d0, 0.0d0, 200.0d0, 0.0d0, 300.0d0/)
type(c_ptr) :: handle = c_null_ptr
real(c_double), pointer :: dXval(:), dY(:)
integer(c_int), pointer :: dXind(:)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dsctr' (Fortran 2008 interfaces) - "
call hipCheck(hipMalloc(dXval, source=hXval))
call hipCheck(hipMalloc(dXind, source=hXind))
call hipCheck(hipMalloc(dY, source=hY))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_dsctr(handle, nnz, dXval, dXind, dY, rocsparse_index_base_zero))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hY, dY, hipMemcpyDeviceToHost))
do i = 1, n
if (abs(hY(i) - hExp(i)) > 1.0d-12) then
write(*,*) "FAILED! y(", i, ") = ", hY(i), " expected ", hExp(i); call exit(1)
end if
end do
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(dXval)); call hipCheck(hipFree(dXind)); call hipCheck(hipFree(dY))
write(*,*) "PASSED!"
end program dsctr
doti computes the dot product of a sparse vector with a dense one, and
gemvi multiplies a dense matrix by a sparse vector,
y = alpha*A*x + beta*y.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Copyright (c) 2020-2022 Advanced Micro Devices, Inc.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program rocsparse_ddoti_test
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: h_xind(3)
real(8) :: h_xval(3), h_y(9)
real(8) :: h_dot
integer, pointer :: d_xind(:)
real(8), pointer :: d_xval(:), d_y(:)
real(8), pointer :: d_dot
integer :: i
integer(c_int) :: M, nnz
type(c_ptr) :: handle
write(*,"(a)",advance="no") "-- Running test 'ddoti' (Fortran 2008 interfaces) - "
! Input data
! Number of rows
M = 9
! Number of non-zero entries
nnz = 3
! Fill structures
h_xind = (/0, 3, 5/)
h_xval = (/1, 2, 3/)
h_y = (/1, 2, 3, 4, 5, 6, 7, 8, 9/)
! Allocate device memory and copy host data to device
call hipCheck(hipMalloc(d_xind, source=h_xind))
call hipCheck(hipMalloc(d_xval,source=h_xval))
call hipCheck(hipMalloc(d_y,source=h_y))
call hipCheck(hipMalloc(d_dot,source=h_dot))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_set_pointer_mode(handle, rocsparse_pointer_mode_device))
! Call ddoti
call rocsparseCheck(rocsparse_ddoti(handle, &
nnz, &
d_xval(1), &
d_xind(1), &
d_y(1), &
d_dot, &
rocsparse_index_base_zero))
! Copy result back to host
call hipCheck(hipMemcpy(h_dot, d_dot, hipMemcpyDeviceToHost))
! Verification
if(h_dot /= 27d0) then
write(*,*) 'FAILED!'
call exit
end if
! Clear rocSPARSE
call rocsparseCheck(rocsparse_destroy_handle(handle))
! Clear device memory
call hipCheck(hipFree(d_xind))
call hipCheck(hipFree(d_xval))
call hipCheck(hipFree(d_y))
call hipCheck(hipFree(d_dot))
! Print success
write(*,*) 'PASSED!'
end program rocsparse_ddoti_test
!!!!!!!!!!!!!!
! rocsparse sgemvi example (dense matrix * sparse vector, single)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Computes y = alpha * A * x + beta * y, where A is a dense m-by-n matrix and x
! is a sparse vector (nnz values x_val at indices x_ind). The result is checked
! against a dense host reference (alpha * matmul(A, x_dense) + beta * y).
!!!!!!!!!!!!!!
!
program dgemvi
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: m = 3, n = 4, lda = 3, nnz = 2
real(c_double) :: hA(3,4) = reshape((/ &
1.0d0, 2.0d0, 3.0d0, 4.0d0, 5.0d0, 6.0d0, 7.0d0, 8.0d0, 9.0d0, 10.0d0, 11.0d0, 12.0d0/), (/3,4/))
real(c_double) :: hXval(2) = (/2.0d0, 3.0d0/)
integer(c_int) :: hXind(2) = (/0, 2/)
real(c_double) :: hY(3) = (/1.0d0, 1.0d0, 1.0d0/)
real(c_double) :: alpha = 2.0d0, beta = 3.0d0
real(c_double) :: xDense(4), hRef(3)
type(c_ptr) :: handle = c_null_ptr
real(c_double), pointer :: dA(:,:), dXval(:), dY(:)
integer(c_int), pointer :: dXind(:)
type(c_ptr) :: dBuf
integer(c_size_t) :: bufSize
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dgemvi' (Fortran 2008 interfaces) - "
xDense = 0.0d0
do i = 1, nnz
xDense(hXind(i) + 1) = hXval(i)
end do
hRef = alpha * matmul(hA, xDense) + beta * hY
call hipCheck(hipMalloc(dA, source=hA))
call hipCheck(hipMalloc(dXval, source=hXval))
call hipCheck(hipMalloc(dXind, source=hXind))
call hipCheck(hipMalloc(dY, source=hY))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_dgemvi_buffer_size(handle, rocsparse_operation_none, m, n, nnz, bufSize))
call hipCheck(hipMalloc(dBuf, max(bufSize, 1_c_size_t)))
call rocsparseCheck(rocsparse_dgemvi(handle, rocsparse_operation_none, m, n, alpha, &
dA, lda, nnz, dXval, dXind, beta, dY, rocsparse_index_base_zero, dBuf))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hY, dY, hipMemcpyDeviceToHost))
do i = 1, m
if (abs(hY(i) - hRef(i)) > 1.0d-11) then
write(*,*) "FAILED! y(", i, ") = ", hY(i), " expected ", hRef(i); call exit(1)
end if
end do
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(dA)); call hipCheck(hipFree(dXval)); call hipCheck(hipFree(dXind))
call hipCheck(hipFree(dY)); call hipCheck(hipFree(dBuf))
write(*,*) "PASSED!"
end program dgemvi
Format conversions#
rocSPARSE converts between the sparse formats. csr2csc converts CSR to CSC,
which is equivalent to transposing the sparse matrix.
!!!!!!!!!!!!!/
! dcsr2csc example (double-precision CSR -> CSC conversion / sparse transpose)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/conversion.html
!
! Converting A from CSR to CSC is equivalent to producing the CSR of A**T.
! We check the resulting csc_col_ptr / csc_row_ind / csc_val against the known
! transpose. csr2csc needs a workspace buffer sized by csr2csc_buffer_size.
!!!!!!!!!!!!!!/
!
program dcsr2csc
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
! 3x3 sparse matrix in CSR (0-based):
! row 0: (0,0)=1, (0,2)=2
! row 1: (1,1)=3
! row 2: (2,0)=4, (2,2)=5
integer(c_int), parameter :: M = 3, N = 3, nnz = 5
integer(c_int) :: h_csr_row_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_csr_col_ind(5) = (/0, 2, 1, 0, 2/)
real(c_double) :: h_csr_val(5) = (/1, 2, 3, 4, 5/)
! Expected CSC (= CSR of the transpose)
integer(c_int) :: h_exp_col_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_exp_row_ind(5) = (/0, 2, 1, 0, 2/)
real(c_double) :: h_exp_val(5) = (/1, 4, 3, 2, 5/)
integer(c_int) :: h_csc_col_ptr(4)
integer(c_int) :: h_csc_row_ind(5)
real(c_double) :: h_csc_val(5)
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:)
integer(c_int), pointer :: d_csc_col_ptr(:), d_csc_row_ind(:)
real(c_double), pointer :: d_csc_val(:)
type(c_ptr) :: handle, d_buffer
integer(c_size_t) :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'rocsparse_dcsr2csc' (Fortran 2008 interfaces) - "
! Allocate device memory and copy the CSR matrix to device
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_csr_col_ind, source=h_csr_col_ind))
call hipCheck(hipMalloc(d_csr_val, source=h_csr_val))
call hipCheck(hipMalloc(d_csc_col_ptr, mold=h_csc_col_ptr))
call hipCheck(hipMalloc(d_csc_row_ind, mold=h_csc_row_ind))
call hipCheck(hipMalloc(d_csc_val, mold=h_csc_val))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
! Query and allocate the required workspace
call rocsparseCheck(rocsparse_csr2csc_buffer_size(handle, M, N, nnz, &
d_csr_row_ptr(1), d_csr_col_ind(1), rocsparse_action_numeric, buffer_size))
call hipCheck(hipMalloc(d_buffer, buffer_size))
! Convert CSR -> CSC (numeric: also permute values)
call rocsparseCheck(rocsparse_dcsr2csc(handle, M, N, nnz, &
d_csr_val(1), d_csr_row_ptr(1), d_csr_col_ind(1), &
d_csc_val(1), d_csc_row_ind(1), d_csc_col_ptr(1), &
rocsparse_action_numeric, rocsparse_index_base_zero, d_buffer))
! Copy the result back to host
call hipCheck(hipMemcpy(h_csc_col_ptr, d_csc_col_ptr, hipMemcpyDeviceToHost))
call hipCheck(hipMemcpy(h_csc_row_ind, d_csc_row_ind, hipMemcpyDeviceToHost))
call hipCheck(hipMemcpy(h_csc_val, d_csc_val, hipMemcpyDeviceToHost))
! Verify structure
do i = 1,N+1
if(h_csc_col_ptr(i) /= h_exp_col_ptr(i)) then
write(*,*) "FAILED! csc_col_ptr(", i, ") = ", h_csc_col_ptr(i), " expected ", h_exp_col_ptr(i)
call exit
end if
end do
do i = 1,nnz
if(h_csc_row_ind(i) /= h_exp_row_ind(i)) then
write(*,*) "FAILED! csc_row_ind(", i, ") = ", h_csc_row_ind(i), " expected ", h_exp_row_ind(i)
call exit
end if
error = abs(h_csc_val(i) - h_exp_val(i))
if(error .gt. error_max) then
write(*,*) "FAILED! csc_val(", i, ") = ", h_csc_val(i), " expected ", h_exp_val(i)
call exit
end if
end do
! Clean up
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_csr_col_ind))
call hipCheck(hipFree(d_csr_val))
call hipCheck(hipFree(d_csc_col_ptr))
call hipCheck(hipFree(d_csc_row_ind))
call hipCheck(hipFree(d_csc_val))
call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program dcsr2csc
csr2coo and coo2csr convert between the CSR row-pointer array and the
COO row-index array, the compressed and expanded forms of the same row
information.
!!!!!!!!!!!!!/
! csr2coo example (CSR -> COO row-index conversion)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/reference/conversion.html
!
! csr2coo expands the CSR row-pointer array into per-nonzero COO row indices.
! It is integer-only (no s/d/c/z variants).
!!!!!!!!!!!!!!/
!
program csr2coo
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: M = 3 ! number of rows
integer(c_int), parameter :: nnz = 4 ! number of nonzeros
! CSR row-pointer for a 3x3 matrix with 2, 1, 1 nonzeros per row (0-based)
integer(c_int) :: h_csr_row_ptr(4) = (/0, 2, 3, 4/)
integer(c_int) :: h_coo_row_ind(4)
integer(c_int) :: h_expected(4) = (/0, 0, 1, 2/)
integer(c_int), pointer :: d_csr_row_ptr(:)
integer(c_int), pointer :: d_coo_row_ind(:)
type(c_ptr) :: handle
write(*,"(a)",advance="no") "-- Running test 'rocsparse_csr2coo' (Fortran 2008 interfaces) - "
! Allocate device memory and copy the row-pointer to device
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_coo_row_ind, mold=h_coo_row_ind))
! Create rocSPARSE handle
call rocsparseCheck(rocsparse_create_handle(handle))
! Convert CSR row pointers to COO row indices
call rocsparseCheck(rocsparse_csr2coo(handle, d_csr_row_ptr(1), nnz, M, &
d_coo_row_ind(1), rocsparse_index_base_zero))
! Copy the result back to host
call hipCheck(hipMemcpy(h_coo_row_ind, d_coo_row_ind, hipMemcpyDeviceToHost))
! Verify against the expected COO row indices
do i = 1,nnz
if(h_coo_row_ind(i) /= h_expected(i)) then
write(*,*) "FAILED! coo_row_ind(", i, ") = ", h_coo_row_ind(i), " expected ", h_expected(i)
call exit
end if
end do
! Clean up
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_csr_row_ptr))
call hipCheck(hipFree(d_coo_row_ind))
write(*,*) "PASSED!"
end program csr2coo
!!!!!!!!!!!!!!
! rocsparse coo2csr example (COO row indices -> CSR row pointers)
! see: https:!rocm.docs.amd.com/projects/rocSPARSE/en/latest/
!
! Compresses a per-nonzero COO row-index array into a CSR row-pointer array and
! checks the offsets. Inverse of csr2coo.
!!!!!!!!!!!!!!
!
program coo2csr
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_rocsparse
implicit none
integer :: i
integer(c_int), parameter :: M = 3, nnz = 5
integer(c_int) :: h_coo_row(5) = (/0, 0, 1, 2, 2/)
integer(c_int) :: h_exp_csr_row_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_csr_row_ptr(4)
integer(c_int), pointer :: d_coo_row(:), d_csr_row_ptr(:)
type(c_ptr) :: handle = c_null_ptr
write(*,"(a)",advance="no") "-- Running test 'rocsparse_coo2csr' (Fortran 2008 interfaces) - "
call hipCheck(hipMalloc(d_coo_row, source=h_coo_row))
call hipCheck(hipMalloc(d_csr_row_ptr, mold=h_csr_row_ptr))
call rocsparseCheck(rocsparse_create_handle(handle))
call rocsparseCheck(rocsparse_coo2csr(handle, d_coo_row, nnz, M, d_csr_row_ptr, rocsparse_index_base_zero))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(h_csr_row_ptr, d_csr_row_ptr, hipMemcpyDeviceToHost))
do i = 1, M + 1
if (h_csr_row_ptr(i) /= h_exp_csr_row_ptr(i)) then
write(*,*) "FAILED! csr_row_ptr(", i, ") = ", h_csr_row_ptr(i), " expected ", h_exp_csr_row_ptr(i); call exit(1)
end if
end do
call rocsparseCheck(rocsparse_destroy_handle(handle))
call hipCheck(hipFree(d_coo_row)); call hipCheck(hipFree(d_csr_row_ptr))
write(*,*) "PASSED!"
end program coo2csr