hipSPARSE examples#
hipSPARSE is a thin
layer over rocSPARSE whose API follows cuSPARSE. hipFORT exposes it through the
hipfort_hipsparse module.
Every program on this page is complete and self-contained, and 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 cuSPARSE-style interface,
see the rocSPARSE examples, where the equivalent
programs are written against the hipfort_rocsparse module.
Where a routine has the four precisions, a program 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 program:
Sparse matrix formats. Most programs 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 programs 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 programs allocate the buffer only when the queried size is positive.Every call returns a status code. The programs wrap hipSPARSE calls in
hipsparseCheckand HIP calls inhipCheckfrom thehipfort_checkmodule, both of which abort on failure.
Building and running#
The programs 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 program 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