hipSPARSE examples#
hipSPARSE is a thin
portability layer over rocSPARSE on AMD GPUs and cuSPARSE on NVIDIA GPUs. Its API
mirrors cuSPARSE, so the same source builds against either backend. hipFORT
exposes it through the hipfort_hipsparse 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/hipsparse 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/hipsparse.
If you want direct access to rocSPARSE rather than a portable interface, the
equivalent programs are written against the hipfort_rocsparse module.
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 hipsparse prefix
letter.
Conventions#
hipSPARSE 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. A few routines take COO (coordinate) row/column arrays.
Zero-based indexing. The examples use
HIPSPARSE_INDEX_BASE_ZERO, so CSR row pointers and column indices start at 0, matching the cuSPARSE samples. The Fortran host arrays that hold them are ordinary 1-based arrays whose values are 0-based.Two API generations. The generic API (SpMV, SpMM, SDDMM, SpSV, SpSM) wraps the operands in matrix/vector descriptors (
hipsparseCreateCsr,hipsparseCreateDnMat,hipsparseCreateDnVec) and runs in stages: query a workspace size, optionally preprocess/analyze, then compute. The older level-2/level-3 API (csrsv2,csrilu02,gemvi) uses an info handle and a matrix descriptor (hipsparseCreateMatDescr).Zero-size buffers. When a workspace query returns 0, pass a null pointer, not an allocated one: hipSPARSE returns
HIPSPARSE_STATUS_INVALID_VALUEif a non-null buffer is supplied for a zero-size workspace. The examples allocate the buffer only when the queried size is positive.Every call returns a status code. The examples wrap hipSPARSE calls in
hipsparseCheckand HIP calls inhipCheckfrom thehipfort_checkmodule, both of which abort on failure.
Building an example#
The examples only need the hipsparse and hip hipFORT components:
find_package(hipfort REQUIRED COMPONENTS hip hipsparse)
add_executable(my_sparse hipsparse_dspmv.f08)
target_link_libraries(my_sparse PRIVATE hipfort::hipsparse 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
SpMV_bufferSize and SpMV stages.
!!!!!!!!!!!!!!
! hipsparse SpMV example (double, y = alpha*A*x + beta*y)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Generic API: build a CSR descriptor for A and dense-vector descriptors for x
! and y, query the workspace with SpMV_bufferSize, then run SpMV. Result is
! checked against A*x. Uses the named algorithm enum HIPSPARSE_SPMV_ALG_DEFAULT.
!
! NOTE: the descriptor constructors are c_ptr-only (no array overloads), so
! device buffers are passed via c_loc(...).
!!!!!!!!!!!!!!
!
program hipsparse_dspmv
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
use hipfort_enums
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(:), 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 'hipsparse_dspmv' (Fortran 2008 interfaces) - "
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))
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseCreateCsr(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), &
HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_BASE_ZERO, HIP_R_64F))
call hipsparseCheck(hipsparseCreateDnVec(vecX, int(N,c_int64_t), c_loc(d_x), HIP_R_64F))
call hipsparseCheck(hipsparseCreateDnVec(vecY, int(M,c_int64_t), c_loc(d_y), HIP_R_64F))
call hipsparseCheck(hipsparseSpMV_bufferSize(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
c_loc(alpha), matA, vecX, c_loc(beta), vecY, HIP_R_64F, HIPSPARSE_SPMV_ALG_DEFAULT, buffer_size))
! hipSPARSE requires a null buffer when the queried size is 0; a non-null
! (dummy) pointer makes SpMV return HIPSPARSE_STATUS_INVALID_VALUE.
d_buffer = c_null_ptr
if (buffer_size > 0) call hipCheck(hipMalloc(d_buffer, buffer_size))
call hipsparseCheck(hipsparseSpMV(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
c_loc(alpha), matA, vecX, c_loc(beta), vecY, HIP_R_64F, HIPSPARSE_SPMV_ALG_DEFAULT, d_buffer))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(h_y, d_y, hipMemcpyDeviceToHost))
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! y(", i, ") = ", h_y(i), " expected ", h_expected(i); call exit(1)
end if
end do
call hipsparseCheck(hipsparseDestroyDnVec(vecX))
call hipsparseCheck(hipsparseDestroyDnVec(vecY))
call hipsparseCheck(hipsparseDestroySpMat(matA))
call hipsparseCheck(hipsparseDestroy(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))
if (c_associated(d_buffer)) call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program hipsparse_dspmv
SpMM multiplies a sparse matrix by a dense matrix,
C = alpha*A*B + beta*C, with dense-matrix descriptors for B and C.
!!!!!!!!!!!!!!
! hipsparse SpMM example (double, C = alpha*A*B + beta*C)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Generic API: build a CSR descriptor for the sparse A and dense-matrix
! descriptors for B and C, query the workspace with SpMM_bufferSize, then run
! SpMM. Result is checked against a dense host reference (matmul(A_dense, B)).
! Dense matrices are column-major (HIPSPARSE_ORDER_COL).
!
! NOTE: the descriptor constructors are c_ptr-only (no array overloads), so
! device buffers are passed via c_loc(...).
!!!!!!!!!!!!!!
!
program hipsparse_dspmm
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
use hipfort_enums
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) and C (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)
real(c_double) :: h_Adense(3,3), h_expected(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(:), 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 'hipsparse_dspmm' (Fortran 2008 interfaces) - "
h_Adense = 0.0_c_double
h_Adense(1,1) = 1; h_Adense(1,3) = 2
h_Adense(2,2) = 3
h_Adense(3,1) = 4; h_Adense(3,3) = 5
h_expected = matmul(h_Adense, h_B)
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))
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseCreateCsr(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), &
HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_BASE_ZERO, HIP_R_64F))
call hipsparseCheck(hipsparseCreateDnMat(matB, int(K,c_int64_t), int(Ncol,c_int64_t), int(K,c_int64_t), &
c_loc(d_B), HIP_R_64F, HIPSPARSE_ORDER_COL))
call hipsparseCheck(hipsparseCreateDnMat(matC, int(M,c_int64_t), int(Ncol,c_int64_t), int(M,c_int64_t), &
c_loc(d_C), HIP_R_64F, HIPSPARSE_ORDER_COL))
call hipsparseCheck(hipsparseSpMM_bufferSize(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
HIPSPARSE_OPERATION_NON_TRANSPOSE, c_loc(alpha), matA, matB, c_loc(beta), matC, &
HIP_R_64F, HIPSPARSE_SPMM_ALG_DEFAULT, buffer_size))
! hipSPARSE requires a null buffer when the queried size is 0; a non-null
! (dummy) pointer makes SpMM return HIPSPARSE_STATUS_INVALID_VALUE.
d_buffer = c_null_ptr
if (buffer_size > 0) call hipCheck(hipMalloc(d_buffer, buffer_size))
call hipsparseCheck(hipsparseSpMM(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
HIPSPARSE_OPERATION_NON_TRANSPOSE, c_loc(alpha), matA, matB, c_loc(beta), matC, &
HIP_R_64F, HIPSPARSE_SPMM_ALG_DEFAULT, d_buffer))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(h_C, d_C, hipMemcpyDeviceToHost))
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! C(", i, j, ") = ", h_C(i,j), " expected ", h_expected(i,j); call exit(1)
end if
end do
end do
call hipsparseCheck(hipsparseDestroyDnMat(matB))
call hipsparseCheck(hipsparseDestroyDnMat(matC))
call hipsparseCheck(hipsparseDestroySpMat(matA))
call hipsparseCheck(hipsparseDestroy(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))
if (c_associated(d_buffer)) call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program hipsparse_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. The example uses dense descriptors for
A and B, a CSR descriptor for C, and the three
SDDMM_bufferSize / SDDMM_preprocess / SDDMM stages.
!!!!!!!!!!!!!!
! hipsparse SDDMM example (d, sampled dense-dense matmul)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! SDDMM computes C = alpha * (A * B) .* spy(C) + beta * C: the dense product
! A*B is evaluated only at the nonzero positions of the sparse (CSR) C. Generic
! API with three stages (SDDMM_bufferSize -> SDDMM_preprocess -> SDDMM). The
! sampled values are checked against matmul(A,B) on the host.
!
! NOTE: the descriptor constructors are c_ptr-only (no array overloads), so
! device buffers are passed via c_loc(...).
!!!!!!!!!!!!!!
!
program hipsparse_dsddmm
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
use hipfort_enums
implicit none
integer :: i
integer(c_int), parameter :: M = 3, N = 2, K = 3, nnz = 4
real(c_double) :: h_A(3,3) = reshape((/1, 4, 7, 2, 5, 8, 3, 6, 10/), (/3,3/))
real(c_double) :: h_B(3,2) = reshape((/1, 3, 5, 2, 4, 6/), (/3,2/))
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 :: error
real, parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'hipsparse_dsddmm' (Fortran 2008 interfaces) - "
h_AB = matmul(h_A, h_B)
h_expected(1) = h_AB(1,1)
h_expected(2) = h_AB(2,2)
h_expected(3) = h_AB(3,1)
h_expected(4) = h_AB(3,2)
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 hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseCreateDnMat(matA, int(M,c_int64_t), int(K,c_int64_t), int(M,c_int64_t), &
c_loc(d_A), HIP_R_64F, HIPSPARSE_ORDER_COL))
call hipsparseCheck(hipsparseCreateDnMat(matB, int(K,c_int64_t), int(N,c_int64_t), int(K,c_int64_t), &
c_loc(d_B), HIP_R_64F, HIPSPARSE_ORDER_COL))
call hipsparseCheck(hipsparseCreateCsr(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), &
HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_BASE_ZERO, HIP_R_64F))
call hipsparseCheck(hipsparseSDDMM_bufferSize(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
HIPSPARSE_OPERATION_NON_TRANSPOSE, c_loc(alpha), matA, matB, c_loc(beta), matC, &
HIP_R_64F, HIPSPARSE_SDDMM_ALG_DEFAULT, buffer_size))
d_buffer = c_null_ptr
if (buffer_size > 0) call hipCheck(hipMalloc(d_buffer, buffer_size))
call hipsparseCheck(hipsparseSDDMM_preprocess(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
HIPSPARSE_OPERATION_NON_TRANSPOSE, c_loc(alpha), matA, matB, c_loc(beta), matC, &
HIP_R_64F, HIPSPARSE_SDDMM_ALG_DEFAULT, d_buffer))
call hipsparseCheck(hipsparseSDDMM(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
HIPSPARSE_OPERATION_NON_TRANSPOSE, c_loc(alpha), matA, matB, c_loc(beta), matC, &
HIP_R_64F, HIPSPARSE_SDDMM_ALG_DEFAULT, d_buffer))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(h_csr_val, d_csr_val, hipMemcpyDeviceToHost))
do i = 1, nnz
error = abs(h_csr_val(i) - h_expected(i)) / max(abs(h_expected(i)), 1.0)
if(error .gt. error_max) then
write(*,*) "FAILED! val(", i, ") = ", h_csr_val(i), " expected ", h_expected(i); call exit(1)
end if
end do
call hipsparseCheck(hipsparseDestroyDnMat(matA))
call hipsparseCheck(hipsparseDestroyDnMat(matB))
call hipsparseCheck(hipsparseDestroySpMat(matC))
call hipsparseCheck(hipsparseDestroy(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))
if (c_associated(d_buffer)) call hipCheck(hipFree(d_buffer))
write(*,*) "PASSED!"
end program hipsparse_dsddmm
Sparse triangular solves#
The generic SpSV solves a sparse triangular system for a single right-hand
side, and SpSM solves it for several right-hand sides at once. Both add an
analysis stage between the buffer-size query and the solve.
!!!!!!!!!!!!!/
! dsptrsv example (double-precision sparse triangular solve, op(A)*y = alpha*x)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! 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 SpMatSetAttribute (fill_mode + diag_type).
! hipSPARSE flow: createDescr -> bufferSize -> analysis -> solve.
!
! NOTE: descriptor/array arguments are c_ptr-only, so device buffers and the
! attribute values are passed via c_loc(...). bufferSize is a c_ptr to a host
! size_t.
!!!!!!!!!!!!!!/
!
program dsptrsv
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
use hipfort_hipsparse_enums
use hipfort_enums
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(HIPSPARSE_FILL_MODE_LOWER)), target :: fill = HIPSPARSE_FILL_MODE_LOWER
integer(kind(HIPSPARSE_DIAG_TYPE_NON_UNIT)), target :: diag = HIPSPARSE_DIAG_TYPE_NON_UNIT
integer(c_int), pointer :: d_csr_row_ptr(:), d_csr_col_ind(:)
real(c_double), pointer :: d_csr_val(:), d_x(:), d_y(:)
type(c_ptr) :: handle = c_null_ptr
type(c_ptr) :: matL, vecX, vecY, spsvDescr, d_buffer
integer(c_size_t), target :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 100 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'hipsparse_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 handle, CSR descriptor for L (lower / non-unit diag), dense vectors
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseCreateCsr(matL, int(M,c_int64_t), int(N,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_csr_row_ptr(1)), c_loc(d_csr_col_ind(1)), c_loc(d_csr_val(1)), &
HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_BASE_ZERO, HIP_R_64F))
call hipsparseCheck(hipsparseSpMatSetAttribute(matL, HIPSPARSE_SPMAT_FILL_MODE, c_loc(fill), int(4,c_size_t)))
call hipsparseCheck(hipsparseSpMatSetAttribute(matL, HIPSPARSE_SPMAT_DIAG_TYPE, c_loc(diag), int(4,c_size_t)))
call hipsparseCheck(hipsparseCreateDnVec(vecX, int(M,c_int64_t), c_loc(d_x(1)), HIP_R_64F))
call hipsparseCheck(hipsparseCreateDnVec(vecY, int(M,c_int64_t), c_loc(d_y(1)), HIP_R_64F))
call hipsparseCheck(hipsparseSpSV_createDescr(spsvDescr))
! Stage 1: workspace size
call hipsparseCheck(hipsparseSpSV_bufferSize(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, c_loc(alpha), matL, vecX, vecY, &
HIP_R_64F, HIPSPARSE_SPSV_ALG_DEFAULT, spsvDescr, buffer_size))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Stage 2: analysis
call hipsparseCheck(hipsparseSpSV_analysis(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, c_loc(alpha), matL, vecX, vecY, &
HIP_R_64F, HIPSPARSE_SPSV_ALG_DEFAULT, spsvDescr, d_buffer))
! Stage 3: solve
call hipsparseCheck(hipsparseSpSV_solve(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, c_loc(alpha), matL, vecX, vecY, &
HIP_R_64F, HIPSPARSE_SPSV_ALG_DEFAULT, spsvDescr))
! 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 hipsparseCheck(hipsparseSpSV_destroyDescr(spsvDescr))
call hipsparseCheck(hipsparseDestroyDnVec(vecX))
call hipsparseCheck(hipsparseDestroyDnVec(vecY))
call hipsparseCheck(hipsparseDestroySpMat(matL))
call hipsparseCheck(hipsparseDestroy(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))
call hipCheck(hipDeviceReset())
write(*,*) "PASSED!"
end program dsptrsv
!!!!!!!!!!!!!/
! dsptrsm example (double-precision sparse triangular solve with multiple rhs)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! 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 SpMatSetAttribute.
! hipSPARSE flow: createDescr -> bufferSize -> analysis -> solve.
!
! NOTE: descriptor/array arguments are c_ptr-only, so device buffers and the
! attribute values are passed via c_loc(...). bufferSize is a c_ptr to a host
! size_t.
!!!!!!!!!!!!!!/
!
program dsptrsm
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
use hipfort_hipsparse_enums
use hipfort_enums
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(HIPSPARSE_FILL_MODE_LOWER)), target :: fill = HIPSPARSE_FILL_MODE_LOWER
integer(kind(HIPSPARSE_DIAG_TYPE_NON_UNIT)), target :: diag = HIPSPARSE_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 = c_null_ptr
type(c_ptr) :: matL, matB, matC, spsmDescr, d_buffer
integer(c_size_t), target :: buffer_size
real(c_double) :: error
real(c_double), parameter :: error_max = 100 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'hipsparse_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 handle, CSR descriptor for L (lower / non-unit diag), dense matrices
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseCreateCsr(matL, int(M,c_int64_t), int(M,c_int64_t), int(nnz,c_int64_t), &
c_loc(d_csr_row_ptr(1)), c_loc(d_csr_col_ind(1)), c_loc(d_csr_val(1)), &
HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_32I, HIPSPARSE_INDEX_BASE_ZERO, HIP_R_64F))
call hipsparseCheck(hipsparseSpMatSetAttribute(matL, HIPSPARSE_SPMAT_FILL_MODE, c_loc(fill), int(4,c_size_t)))
call hipsparseCheck(hipsparseSpMatSetAttribute(matL, HIPSPARSE_SPMAT_DIAG_TYPE, c_loc(diag), int(4,c_size_t)))
call hipsparseCheck(hipsparseCreateDnMat(matB, int(M,c_int64_t), int(nrhs,c_int64_t), int(M,c_int64_t), &
c_loc(d_X(1,1)), HIP_R_64F, HIPSPARSE_ORDER_COLUMN))
call hipsparseCheck(hipsparseCreateDnMat(matC, int(M,c_int64_t), int(nrhs,c_int64_t), int(M,c_int64_t), &
c_loc(d_C(1,1)), HIP_R_64F, HIPSPARSE_ORDER_COLUMN))
call hipsparseCheck(hipsparseSpSM_createDescr(spsmDescr))
! Stage 1: workspace size
call hipsparseCheck(hipsparseSpSM_bufferSize(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
c_loc(alpha), matL, matB, matC, HIP_R_64F, HIPSPARSE_SPSM_ALG_DEFAULT, spsmDescr, buffer_size))
call hipCheck(hipMalloc(d_buffer, max(buffer_size, 1_c_size_t)))
! Stage 2: analysis
call hipsparseCheck(hipsparseSpSM_analysis(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
c_loc(alpha), matL, matB, matC, HIP_R_64F, HIPSPARSE_SPSM_ALG_DEFAULT, spsmDescr, d_buffer))
! Stage 3: solve
call hipsparseCheck(hipsparseSpSM_solve(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
c_loc(alpha), matL, matB, matC, HIP_R_64F, HIPSPARSE_SPSM_ALG_DEFAULT, spsmDescr, 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, " C(", i, ",", j, ") = ", h_C(i,j)
call exit
end if
end do
end do
! Clean up
call hipsparseCheck(hipsparseSpSM_destroyDescr(spsmDescr))
call hipsparseCheck(hipsparseDestroyDnMat(matB))
call hipsparseCheck(hipsparseDestroyDnMat(matC))
call hipsparseCheck(hipsparseDestroySpMat(matL))
call hipsparseCheck(hipsparseDestroy(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))
call hipCheck(hipDeviceReset())
write(*,*) "PASSED!"
end program dsptrsm
The older csrsv2 triangular solve uses the info-handle API instead: create
a matrix descriptor and a csrsv2 info object, query the buffer size, run the
analysis phase, then solve.
!!!!!!!!!!!!!!
! hipsparse Dcsrsv2 example (sparse triangular solve v2, single)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Solves the lower-triangular system L*x = alpha*f for x using the legacy
! csrsv2 API (bufferSize -> analysis -> solve, with a mat descriptor and a
! csrsv2Info object). The right-hand side f is built from a known solution so
! the recovered x can be checked directly.
!
! L = [ 2 0 0 ] x = [ 1 ] f = L*x = [ 2 ]
! [ 1 2 0 ] [ 2 ] [ 5 ]
! [ 3 1 2 ] [ 3 ] [ 11 ]
!!!!!!!!!!!!!!
!
program hipsparse_dcsrsv2
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
implicit none
integer :: i
integer(c_int), parameter :: m = 3, nnz = 6
integer(c_int) :: hRowPtr(4) = (/0, 1, 3, 6/)
integer(c_int) :: hColInd(6) = (/0, 0, 1, 0, 1, 2/)
real(c_double) :: hVal(6) = (/2.0d0, 1.0d0, 2.0d0, 3.0d0, 1.0d0, 2.0d0/)
real(c_double) :: hF(3) = (/2.0d0, 5.0d0, 11.0d0/)
real(c_double) :: hX(3)
real(c_double) :: hExp(3) = (/1.0d0, 2.0d0, 3.0d0/)
real(c_double) :: alpha = 1.0d0
type(c_ptr) :: handle = c_null_ptr
type(c_ptr) :: descrA = c_null_ptr
type(c_ptr) :: info = c_null_ptr
integer(c_int), pointer :: dRowPtr(:), dColInd(:)
real(c_double), pointer :: dVal(:), dF(:), dX(:)
type(c_ptr) :: dBuf
integer(c_int) :: bufSize
write(*,"(a)",advance="no") "-- Running test 'hipsparse_dcsrsv2' (Fortran 2008 interfaces) - "
call hipCheck(hipMalloc(dRowPtr, source=hRowPtr))
call hipCheck(hipMalloc(dColInd, source=hColInd))
call hipCheck(hipMalloc(dVal, source=hVal))
call hipCheck(hipMalloc(dF, source=hF))
call hipCheck(hipMalloc(dX, mold=hX))
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseCreateMatDescr(descrA))
call hipsparseCheck(hipsparseSetMatType(descrA, HIPSPARSE_MATRIX_TYPE_GENERAL))
call hipsparseCheck(hipsparseSetMatIndexBase(descrA, HIPSPARSE_INDEX_BASE_ZERO))
call hipsparseCheck(hipsparseSetMatFillMode(descrA, HIPSPARSE_FILL_MODE_LOWER))
call hipsparseCheck(hipsparseSetMatDiagType(descrA, HIPSPARSE_DIAG_TYPE_NON_UNIT))
call hipsparseCheck(hipsparseCreateCsrsv2Info(info))
call hipsparseCheck(hipsparseDcsrsv2_bufferSize(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
m, nnz, descrA, dVal, dRowPtr, dColInd, info, bufSize))
call hipCheck(hipMalloc(dBuf, int(max(bufSize,1),c_size_t)))
call hipsparseCheck(hipsparseDcsrsv2_analysis(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
m, nnz, descrA, dVal, dRowPtr, dColInd, info, HIPSPARSE_SOLVE_POLICY_NO_LEVEL, dBuf))
call hipsparseCheck(hipsparseDcsrsv2_solve(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
m, nnz, alpha, descrA, dVal, dRowPtr, dColInd, info, dF, dX, &
HIPSPARSE_SOLVE_POLICY_NO_LEVEL, dBuf))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hX, dX, hipMemcpyDeviceToHost))
do i = 1, m
if (abs(hX(i) - hExp(i)) > 1.0d-12) then
write(*,*) "FAILED! x(", i, ") = ", hX(i), " expected ", hExp(i); call exit(1)
end if
end do
call hipsparseCheck(hipsparseDestroyCsrsv2Info(info))
call hipsparseCheck(hipsparseDestroyMatDescr(descrA))
call hipsparseCheck(hipsparseDestroy(handle))
call hipCheck(hipFree(dRowPtr)); call hipCheck(hipFree(dColInd)); call hipCheck(hipFree(dVal))
call hipCheck(hipFree(dF)); call hipCheck(hipFree(dX)); call hipCheck(hipFree(dBuf))
write(*,*) "PASSED!"
end program hipsparse_dcsrsv2
Sparse matrix-matrix multiplication#
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 = A*B)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Two-phase flow with the classic hipSPARSE API: XcsrgemmNnz fills row_ptr_C and
! the total nnz, then Xcsrgemm computes the values. Here B = A, so C = A*A,
! checked against the known product.
!
! NOTE: the mat-descr/array arguments are c_ptr-only, so device buffers are
! passed via c_loc(...). nnzTotal is a host int (hipSPARSE default pointer mode).
!!!!!!!!!!!!!!/
!
program dcsrgemm
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
use hipfort_hipsparse_enums
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/)
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(:)
type(c_ptr) :: handle = c_null_ptr
type(c_ptr) :: 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 'hipsparse_dcsrgemm' (Fortran 2008 interfaces) - "
! Create handle and matrix descriptors
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseCreateMatDescr(descr_A))
call hipsparseCheck(hipsparseCreateMatDescr(descr_B))
call hipsparseCheck(hipsparseCreateMatDescr(descr_C))
! Allocate device memory and copy A (B aliases A)
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))
! Phase 1: compute the sparsity of C (row_ptr_C + total nnz_C)
call hipsparseCheck(hipsparseXcsrgemmNnz(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
M, N, K, descr_A, nnz_A, c_loc(d_csr_row_ptr(1)), c_loc(d_csr_col_ind(1)), &
descr_B, nnz_A, c_loc(d_csr_row_ptr(1)), c_loc(d_csr_col_ind(1)), &
descr_C, c_loc(d_row_ptr_C(1)), nnz_C))
! 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 hipsparseCheck(hipsparseDcsrgemm(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
M, N, K, descr_A, nnz_A, c_loc(d_csr_val(1)), c_loc(d_csr_row_ptr(1)), c_loc(d_csr_col_ind(1)), &
descr_B, nnz_A, c_loc(d_csr_val(1)), c_loc(d_csr_row_ptr(1)), c_loc(d_csr_col_ind(1)), &
descr_C, c_loc(d_val_C(1)), c_loc(d_row_ptr_C(1)), c_loc(d_col_ind_C(1))))
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 hipsparseCheck(hipsparseDestroyMatDescr(descr_A))
call hipsparseCheck(hipsparseDestroyMatDescr(descr_B))
call hipsparseCheck(hipsparseDestroyMatDescr(descr_C))
call hipsparseCheck(hipsparseDestroy(handle))
call hipCheck(hipDeviceReset())
write(*,*) "PASSED!"
end program dcsrgemm
Incomplete-LU preconditioner#
csrilu02 computes an incomplete LU factorization with zero fill-in, used as
a preconditioner. It follows the info-handle pattern: a buffer-size query, an
analysis phase that inspects the pattern, and the factorization itself, with a
zero-pivot query to detect breakdown.
!!!!!!!!!!!!!!
! hipsparse Dcsrilu02 example (incomplete LU, single)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Computes the ILU(0) factorization of a sparse matrix in place using the legacy
! csrilu02 API (bufferSize -> analysis -> compute, with a mat descriptor and a
! csrilu02Info object). For a tridiagonal matrix there is no fill-in, so ILU(0)
! equals the exact LU factorization and the overwritten CSR values can be 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 ] (l21=1/4, l32=4/15,
! [ 1 4 1 ] [ 1/4 15/4 1 ] u22=15/4, u33=56/15)
! [ 0 1 4 ] [ 0 4/15 56/15]
!!!!!!!!!!!!!!
!
program hipsparse_dcsrilu02
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
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) :: descrA = c_null_ptr
type(c_ptr) :: info = c_null_ptr
integer(c_int), pointer :: dRowPtr(:), dColInd(:)
real(c_double), pointer :: dVal(:)
type(c_ptr) :: dBuf
integer(c_int) :: bufSize
write(*,"(a)",advance="no") "-- Running test 'hipsparse_dcsrilu02' (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 hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseCreateMatDescr(descrA))
call hipsparseCheck(hipsparseSetMatType(descrA, HIPSPARSE_MATRIX_TYPE_GENERAL))
call hipsparseCheck(hipsparseSetMatIndexBase(descrA, HIPSPARSE_INDEX_BASE_ZERO))
call hipsparseCheck(hipsparseCreateCsrilu02Info(info))
call hipsparseCheck(hipsparseDcsrilu02_bufferSize(handle, m, nnz, descrA, dVal, dRowPtr, &
dColInd, info, bufSize))
call hipCheck(hipMalloc(dBuf, int(max(bufSize,1),c_size_t)))
call hipsparseCheck(hipsparseDcsrilu02_analysis(handle, m, nnz, descrA, dVal, dRowPtr, &
dColInd, info, HIPSPARSE_SOLVE_POLICY_NO_LEVEL, dBuf))
call hipsparseCheck(hipsparseDcsrilu02(handle, m, nnz, descrA, dVal, dRowPtr, dColInd, &
info, HIPSPARSE_SOLVE_POLICY_NO_LEVEL, dBuf))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(hOut, dVal, hipMemcpyDeviceToHost))
do i = 1, nnz
if (abs(hOut(i) - hExp(i)) > 1.0d-12) then
write(*,*) "FAILED! val(", i, ") = ", hOut(i), " expected ", hExp(i); call exit(1)
end if
end do
call hipsparseCheck(hipsparseDestroyCsrilu02Info(info))
call hipsparseCheck(hipsparseDestroyMatDescr(descrA))
call hipsparseCheck(hipsparseDestroy(handle))
call hipCheck(hipFree(dRowPtr)); call hipCheck(hipFree(dColInd)); call hipCheck(hipFree(dVal))
call hipCheck(hipFree(dBuf))
write(*,*) "PASSED!"
end program hipsparse_dcsrilu02
Tridiagonal solver#
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, hipSPARSE)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! 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 hipsparse_sgtsv_test
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
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 hipsparseCheck(hipsparseCreate(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 hipsparseCheck(hipsparseSgtsv2_bufferSizeExt(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 hipsparseCheck(hipsparseSgtsv2(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 hipsparseCheck(hipsparseDestroy(handle))
deallocate(hdl, hd, hdu, hB)
write(*,*) "PASSED!"
end program hipsparse_sgtsv_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.
!!!!!!!!!!!!!!
! hipsparse Dgthr example (gather y[xInd] -> xVal, double)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Gathers the entries of a dense vector y at the sparse index set xInd into the
! packed vector xVal, then checks the gathered values.
!!!!!!!!!!!!!!
!
program hipsparse_dgthr
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
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 'hipsparse_dgthr' (Fortran 2008 interfaces) - "
call hipCheck(hipMalloc(dY, source=hY))
call hipCheck(hipMalloc(dXind, source=hXind))
call hipCheck(hipMalloc(dXval, mold=hXval))
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseDgthr(handle, nnz, dY, dXval, dXind, HIPSPARSE_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! xVal(", i, ") = ", hXval(i), " expected ", hExp(i); call exit(1)
end if
end do
call hipsparseCheck(hipsparseDestroy(handle))
call hipCheck(hipFree(dY)); call hipCheck(hipFree(dXind)); call hipCheck(hipFree(dXval))
write(*,*) "PASSED!"
end program hipsparse_dgthr
!!!!!!!!!!!!!!
! hipsparse Dsctr example (scatter xVal -> y[xInd], double)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Scatters the packed vector xVal into the dense vector y at the sparse index
! set xInd, then checks the resulting dense vector.
!!!!!!!!!!!!!!
!
program hipsparse_dsctr
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
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 'hipsparse_dsctr' (Fortran 2008 interfaces) - "
call hipCheck(hipMalloc(dXval, source=hXval))
call hipCheck(hipMalloc(dXind, source=hXind))
call hipCheck(hipMalloc(dY, source=hY))
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseDsctr(handle, nnz, dXval, dXind, dY, HIPSPARSE_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 hipsparseCheck(hipsparseDestroy(handle))
call hipCheck(hipFree(dXval)); call hipCheck(hipFree(dXind)); call hipCheck(hipFree(dY))
write(*,*) "PASSED!"
end program hipsparse_dsctr
gemvi multiplies a dense matrix by a sparse vector,
y = alpha*A*x + beta*y, sizing its workspace with a gemvi_bufferSize
query.
!!!!!!!!!!!!!!
! hipsparse Dgemvi example (dense matrix * sparse vector, double)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/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 xVal at indices xInd). The result is checked
! against a dense host reference (alpha * matmul(A, x_dense) + beta * y).
!!!!!!!!!!!!!!
!
program hipsparse_dgemvi
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
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_int) :: bufSize
write(*,"(a)",advance="no") "-- Running test 'hipsparse_dgemvi' (Fortran 2008 interfaces) - "
! Dense host reference.
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 hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseDgemvi_bufferSize(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, &
m, n, nnz, bufSize))
call hipCheck(hipMalloc(dBuf, int(max(bufSize,1),c_size_t)))
call hipsparseCheck(hipsparseDgemvi(handle, HIPSPARSE_OPERATION_NON_TRANSPOSE, m, n, alpha, &
dA, lda, nnz, dXval, dXind, beta, dY, HIPSPARSE_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 hipsparseCheck(hipsparseDestroy(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 hipsparse_dgemvi
Format conversions#
hipSPARSE 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/hipSPARSE/en/latest/
!
! 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.
!
! NOTE: csr2csc args are c_ptr-only, so device buffers are passed via c_loc(...).
!!!!!!!!!!!!!!/
!
program dcsr2csc
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
use hipfort_hipsparse_enums
implicit none
integer :: i
! 3x3 sparse matrix 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/)
! 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 = c_null_ptr
real(c_double) :: error
real(c_double), parameter :: error_max = 10 * epsilon(error_max)
write(*,"(a)",advance="no") "-- Running test 'hipsparse_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 handle and convert CSR -> CSC (numeric: also permute values)
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseDcsr2csc(handle, M, N, nnz, &
c_loc(d_csr_val(1)), c_loc(d_csr_row_ptr(1)), c_loc(d_csr_col_ind(1)), &
c_loc(d_csc_val(1)), c_loc(d_csc_row_ind(1)), c_loc(d_csc_col_ptr(1)), &
HIPSPARSE_ACTION_NUMERIC, HIPSPARSE_INDEX_BASE_ZERO))
call hipCheck(hipDeviceSynchronize())
! 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 and values
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 hipsparseCheck(hipsparseDestroy(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(hipDeviceReset())
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. These are index-only conversions, so they have a single X
(type-agnostic) entry point rather than one per precision.
!!!!!!!!!!!!!!
! hipsparse Xcsr2coo example (CSR row pointers -> COO row indices)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Expands the CSR row-pointer array into one row index per nonzero and checks
! it against the expected COO row indices.
!!!!!!!!!!!!!!
!
program hipsparse_xcsr2coo
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
implicit none
integer :: i
integer(c_int), parameter :: M = 3, nnz = 5
integer(c_int) :: h_csr_row_ptr(4) = (/0, 2, 3, 5/)
integer(c_int) :: h_exp_coo_row(5) = (/0, 0, 1, 2, 2/)
integer(c_int) :: h_coo_row(5)
integer(c_int), pointer :: d_csr_row_ptr(:), d_coo_row(:)
type(c_ptr) :: handle = c_null_ptr
write(*,"(a)",advance="no") "-- Running test 'hipsparse_xcsr2coo' (Fortran 2008 interfaces) - "
call hipCheck(hipMalloc(d_csr_row_ptr, source=h_csr_row_ptr))
call hipCheck(hipMalloc(d_coo_row, mold=h_coo_row))
call hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseXcsr2coo(handle, d_csr_row_ptr, nnz, M, d_coo_row, HIPSPARSE_INDEX_BASE_ZERO))
call hipCheck(hipDeviceSynchronize())
call hipCheck(hipMemcpy(h_coo_row, d_coo_row, hipMemcpyDeviceToHost))
do i = 1, nnz
if (h_coo_row(i) /= h_exp_coo_row(i)) then
write(*,*) "FAILED! coo_row(", i, ") = ", h_coo_row(i), " expected ", h_exp_coo_row(i); call exit(1)
end if
end do
call hipsparseCheck(hipsparseDestroy(handle))
call hipCheck(hipFree(d_csr_row_ptr)); call hipCheck(hipFree(d_coo_row))
write(*,*) "PASSED!"
end program hipsparse_xcsr2coo
!!!!!!!!!!!!!!
! hipsparse Xcoo2csr example (COO row indices -> CSR row pointers)
! see: https:!rocm.docs.amd.com/projects/hipSPARSE/en/latest/
!
! Compresses the per-nonzero COO row-index array into a CSR row-pointer array
! and checks it against the expected offsets. Inverse of Xcsr2coo.
!!!!!!!!!!!!!!
!
program hipsparse_xcoo2csr
use iso_c_binding
use hipfort
use hipfort_check
use hipfort_hipsparse
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 'hipsparse_xcoo2csr' (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 hipsparseCheck(hipsparseCreate(handle))
call hipsparseCheck(hipsparseXcoo2csr(handle, d_coo_row, nnz, M, d_csr_row_ptr, HIPSPARSE_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 hipsparseCheck(hipsparseDestroy(handle))
call hipCheck(hipFree(d_coo_row)); call hipCheck(hipFree(d_csr_row_ptr))
write(*,*) "PASSED!"
end program hipsparse_xcoo2csr