hipBLAS examples

hipBLAS examples#

hipBLAS is a thin layer over rocBLAS whose API follows cuBLAS. hipFORT exposes it through the hipfort_hipblas module, which mirrors the hipBLAS C API one to one and re-exports the enumerations (HIPBLAS_OP_N, HIPBLAS_FILL_MODE_LOWER and so on) from hipfort_hipblas_enums.

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/hipblas 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/hipblas. The packed triangular solve, stpsv, is Fortran 2008 only.

If you want direct access to rocBLAS rather than a cuBLAS-style interface, see the rocBLAS examples, where the equivalent programs are written against the hipfort_rocblas module.

The examples are grouped the way the BLAS routines themselves are: Level 1 operates on vectors, Level 2 on a matrix and a vector, and Level 3 on two matrices. Most routines are provided in four precisions, identified by the usual BLAS prefix: s (single-precision real), d (double-precision real), c (single-precision complex) and z (double-precision complex). Each section below shows one precision and names the sibling test files that cover the others.

hipBLAS call pattern#

A hipBLAS program always follows the same sequence:

  1. Create a handle with hipblasCreate.

  2. Allocate device memory and copy the input data over, either with hipMalloc/hipMemcpy or, in the Fortran 2008 interfaces, with the hipMalloc(source=...) shortcut that allocates and copies in one call.

  3. Call the hipBLAS routine.

  4. Call hipDeviceSynchronize before reading a result, whether it was written to a host scalar or to device memory.

  5. Copy device results back to the host.

  6. Free the device memory and release the handle with hipblasDestroy.

Keep the following conventions in mind:

  • hipBLAS starts in host pointer mode, so alpha, beta and scalar results such as a dot product are read from host memory without any setup call. Only a program that wants those values to live on the device needs hipblasSetPointerMode. This is the main day-to-day difference from rocBLAS, whose examples set the pointer mode explicitly.

  • hipBLAS matrices are stored column-major, which matches Fortran’s native array layout directly, so no transpose trick is needed to call hipBLAS from Fortran.

  • The leading dimension of a device matrix is usually just its first dimension, size(dA,1).

  • Enumerations such as HIPBLAS_OP_N, HIPBLAS_FILL_MODE_LOWER, HIPBLAS_DIAG_NON_UNIT and HIPBLAS_SIDE_LEFT come from the hipfort_hipblas_enums module, which hipfort_hipblas re-exports, so use hipfort_hipblas on its own is enough.

  • Every hipBLAS call returns a status code. The examples wrap them in hipblasCheck from the hipfort_check module, which aborts on failure. (sgemv.f08 and sger.f08 route their hipBLAS calls through hipCheck instead; both abort on a non-zero status.)

Building an example#

The examples only need the hipblas and hip hipFORT components:

find_package(hipfort REQUIRED COMPONENTS hip hipblas)

add_executable(my_blas saxpy.f08)
target_link_libraries(my_blas PRIVATE hipfort::hipblas hipfort::hip)

See Using hipFORT in your application for the full set of build options.

Level 1: vector operations#

Scaled vector update (axpy)#

hipblas?axpy computes y := alpha * x + y. This example runs it in single precision and checks the result against the expected value.

program hip_saxpy
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas


  implicit none

  integer :: n = 6
  type(c_ptr) :: handle = c_null_ptr
  integer :: j
  real, allocatable, dimension(:) :: x, y, y_exact

  real, parameter :: alpha = 2.0
  real, pointer, dimension(:) :: dx, dy

  real :: error
  real, parameter :: error_max = 10*epsilon(error)

  allocate(x(n))
  allocate(y(n))
  allocate(y_exact(n))

  do j = 1,n
    x(j) = j
    y(j) = j
  end do

  do j = 1,n
    y_exact(j) = alpha*x(j) + y(j)
  end do
 
  write(*,"(a)",advance="no") "-- Running test 'SAXPY' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  call hipCheck(hipMalloc(dx,shape(x)))
  call hipCheck(hipMalloc(dy,shape(y)))

  call hipCheck(hipMemcpy(dx, x, hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(dy, y, hipMemcpyHostToDevice))

  call hipblasCheck(hipblasSaxpy(handle,n,alpha,dx,1,dy,1))

  call hipCheck(hipDeviceSynchronize())

  call hipCheck(hipMemcpy(y, dy, hipMemcpyDeviceToHost))

  do j = 1,n
    error = abs((y_exact(j) - y(j))/y_exact(j))
      if( error > error_max )then
        write(*,*) "FAILED! Error bigger than max! Error = ", error
        call exit(1)
      end if
  end do

  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))

  call hipblasCheck(hipblasDestroy(handle))

  deallocate(x,y)

  write(*,*) "PASSED!"

end program hip_saxpy

test/f2008/hipblas/daxpy.f08, caxpy.f08 and zaxpy.f08 run the same computation in double-precision real, single-precision complex and double-precision complex, respectively.

Vector scaling (scal)#

hipblas?scal computes x := alpha * x. The operation is in place, so the device vector is both the input and the output.

program hip_dscal

  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer, parameter :: N = 10240;

  double precision, parameter :: alpha = 10.d0

  double precision,pointer,dimension(:) :: dx

  double precision,allocatable,target,dimension(:) :: hx
  double precision,allocatable,dimension(:) ::hx_scaled

  double precision :: error
  double precision, parameter :: error_max = 10*epsilon(error_max)

  type(c_ptr) :: hip_blas_handle = c_null_ptr

  integer :: i

  write(*,"(a)",advance="no") "-- Running test 'dscal' (Fortran 2008 interfaces) - "
  allocate(hx(N))
  allocate(hx_scaled(N))

  hx(:) = 10.d0

  hx_scaled = alpha*hx

  call hipblasCheck(hipblasCreate(hip_blas_handle))

  call hipCheck(hipMalloc(dx,shape(hx)))

   ! Transfer data from host to device memory
  call hipCheck(hipMemcpy(dx, hx, hipMemcpyHostToDevice))

  call hipblasCheck(hipblasDscal(hip_blas_handle, N, alpha, dx, 1))

  call hipCheck(hipDeviceSynchronize())

  ! Transfer data back to host memory
  call hipCheck(hipMemcpy(hx, dx, hipMemcpyDeviceToHost))

  call hipCheck(hipFree(dx))

  ! Verification
  do i = 1,N
     error = abs(hx(i) - hx_scaled(i) )
     if( error .gt. error_max ) then
        write(*,*) "FAILED! Error bigger than max! Error = ", error, " hx(i) = ", hx(i)
        call exit
     endif
  end do

  call hipblasCheck(hipblasDestroy(hip_blas_handle))

  deallocate(hx_scaled,hx)

  write(*,*) "PASSED!"

end program hip_dscal

This is the only scal example in the test suite; there is no single-precision or complex counterpart. The rocBLAS examples cover the remaining precisions, including the mixed-precision forms that scale a complex vector by a real scalar.

Vector copy and swap#

hipblas?copy copies x into y on the device, leaving x untouched.

program hip_scopy
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer :: n = 6
  type(c_ptr) :: handle = c_null_ptr
  integer :: j
  real, allocatable, target, dimension(:) :: x, y
  integer, parameter :: bytes_per_element = 4 !float precision
 
  real, pointer, dimension(:) :: dx,dy

  real :: error
  real, parameter :: error_max = 10*epsilon(error)

  allocate(x(n))
  allocate(y(n))

  do j = 1,n
    x(j) = j
   ! write(*,*) "value of x(j)" , x(j)
  end do

  write(*,"(a)",advance="no") "-- Running test 'Scopy' (Fortran 2008 interfaces) - "

  do j = 1,n
    y(j) = x(j)
   ! write(*,*) "value of y(j)" , y(j)
  end do
  
  call hipblasCheck(hipblasCreate(handle))

  call hipCheck(hipMalloc(dx,source=x))
  call hipCheck(hipMalloc(dy,source=y))

  call hipblasCheck(hipblasScopy(handle,n,dx,1,dy,1))

  call hipCheck(hipDeviceSynchronize())

  call hipCheck(hipMemcpy(y, dy, hipMemcpyDeviceToHost))

  do j = 1,n
    error = abs(y(j) - x(j))
      if( error > error_max )then
        write(*,*) "FAILED! Error bigger than max! Error = ", error
        call exit(1)
      end if
  end do

  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))

  call hipblasCheck(hipblasDestroy(handle))

  deallocate(x,y)

  write(*,*) "PASSED!"
  
end program hip_scopy

hipblas?swap exchanges the contents of the two vectors instead, so both are modified. The program below checks that x and y have traded values.

program hip_sswap
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none
  
  integer :: n = 6
  integer :: i
  real(kind=4), allocatable, target, dimension(:) :: x, y, x_exact, y_exact
  real(kind=4), pointer, dimension(:) :: dx, dy
  type(c_ptr) :: handle = c_null_ptr

  real :: error
  real, parameter :: error_max = 10*epsilon(error)
  
  allocate(x(n))
  allocate(y(n))
  allocate(x_exact(n))
  allocate(y_exact(n))
  
  do i = 1,n
    x(i) = i
    x_exact(i) = x(i)
    y(i) = 2 * i
    y_exact(i) = y(i)
  end do
  
  write(*,"(a)",advance="no") "-- Running test 'SSWAP' (Fortran 2008 interfaces) - "
  
  call hipblasCheck(hipblasCreate(handle))
  
  call hipCheck(hipMalloc(dx,shape(x)))  
  call hipCheck(hipMalloc(dy,shape(y)))
  
  call hipCheck(hipMemcpy(dx, x, hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(dy, y, hipMemcpyHostToDevice))
  
  call hipblasCheck(hipblasSswap(handle,n,dx,1,dy,1))
  
  call hipCheck(hipMemcpy(x, dx, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(y, dy, hipMemcpyDeviceToHost))
  
  do i = 1,n
    error = MAX(abs((y_exact(i) - x(i))/y_exact(i)), abs((x_exact(i) - y(i))/x_exact(i)))
      if( error > error_max )then
        write(*,*) "FAILED! Error bigger than max! Error = ", error
        call exit(1)
      end if
  end do
  
  call hipblasCheck(hipblasDestroy(handle))
  
  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))
  
  deallocate(x,y)
  
  write(*,*) "PASSED!"
  
end program hip_sswap

Both are single-precision only in the test suite, and neither has a rocBLAS counterpart there.

Dot products#

hipblas?dot computes the dot product of two real vectors and returns the scalar result through a pointer. Because hipBLAS is in host pointer mode by default, that pointer is an ordinary host variable.

program hip_sdot
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer, parameter :: n = 10240
  real(kind=4), allocatable, dimension(:) :: hx, hy
  real(kind=4), pointer, dimension(:) :: dx, dy
  real(kind=4), target :: res
  real(kind=4) :: res_exact, error
  real(kind=4), parameter :: error_max = 10*epsilon(error)
  type(c_ptr) :: handle = c_null_ptr

  write(*,"(a)",advance="no") "-- Running test 'SDOT' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  allocate(hx(n), hy(n))
  hx = 1.0
  hy = 2.0
  res_exact = 2.0 * n

  call hipCheck(hipMalloc(dx, shape(hx)))
  call hipCheck(hipMalloc(dy, shape(hy)))
  call hipCheck(hipMemcpy(dx, hx, hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(dy, hy, hipMemcpyHostToDevice))

  res = 0.0
  call hipblasCheck(hipblasSdot(handle, n, dx, 1, dy, 1, c_loc(res)))
  call hipCheck(hipDeviceSynchronize())

  error = abs((res_exact - res) / res_exact)
  if (error > error_max) then
    write(*,*) "FAILED! error = ", error, " result = ", res
    call exit(1)
  end if

  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))
  call hipblasCheck(hipblasDestroy(handle))
  deallocate(hx, hy)
  write(*,*) "PASSED!"

end program hip_sdot

test/f2008/hipblas/ddot.f08 is the double-precision equivalent.

Complex vectors have two dot product variants: hipblas?dotc conjugates the first vector before multiplying, and hipblas?dotu does not. Compare test/f2008/hipblas/cdotc.f08 and cdotu.f08 (also available in double precision as zdotc.f08 and zdotu.f08) to see the different expected results for the same input data.

Euclidean norm (nrm2)#

hipblas?nrm2 computes the Euclidean norm of a vector and returns it through a pointer, which in host pointer mode is an ordinary host variable.

program hip_snrm2
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas
  use hipfort_hipblas_enums

  implicit none

  ! nrm2(x) with x = 1 over n elements, so the result is sqrt(n).
  integer, parameter :: n = 1024
  real(c_float), allocatable, dimension(:) :: hx
  real(c_float), target :: res
  real(c_float) :: res_exact, error
  real(c_float), parameter :: error_max = 10*epsilon(error)
  real(c_float), pointer, dimension(:) :: dx
  type(c_ptr) :: handle = c_null_ptr

  write(*,"(a)",advance="no") "-- Running test 'snrm2' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  allocate(hx(n))
  hx = 1.0
  res_exact = sqrt(real(n, kind=kind(res_exact)))

  call hipCheck(hipMalloc(dx, shape(hx)))
  call hipCheck(hipMemcpy(dx, hx, hipMemcpyHostToDevice))

  res = 0.0
  call hipblasCheck(hipblasSnrm2(handle, n, dx, 1, c_loc(res)))
  call hipCheck(hipDeviceSynchronize())

  error = abs((res_exact - res) / res_exact)
  if (error > error_max) then
    write(*,*) "FAILED! error = ", error, " result = ", res
    call exit(1)
  end if

  call hipCheck(hipFree(dx))
  call hipblasCheck(hipblasDestroy(handle))

  write(*,*) "PASSED!"

end program hip_snrm2

test/f2008/hipblas/dnrm2.f08 is the double-precision equivalent. The complex forms are named for both types involved, because the norm of a complex vector is real: hipblasScnrm2 takes a single-precision complex vector and returns a single-precision real result, and hipblasDznrm2 is its double-precision counterpart. See scnrm2.f08 and dznrm2.f08.

Sum of absolute values (asum)#

hipblas?asum sums the absolute values of a vector’s elements. For complex vectors it sums abs(real(x)) + abs(aimag(x)) per element rather than the complex modulus.

program hip_sasum
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas
  use hipfort_hipblas_enums

  implicit none

  ! asum(x) with x = 1 over n elements, so the result is n.
  integer, parameter :: n = 1024
  real(c_float), allocatable, dimension(:) :: hx
  real(c_float), target :: res
  real(c_float) :: res_exact, error
  real(c_float), parameter :: error_max = 10*epsilon(error)
  real(c_float), pointer, dimension(:) :: dx
  type(c_ptr) :: handle = c_null_ptr

  write(*,"(a)",advance="no") "-- Running test 'sasum' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  allocate(hx(n))
  hx = 1.0
  res_exact = real(n, kind=kind(res_exact))

  call hipCheck(hipMalloc(dx, shape(hx)))
  call hipCheck(hipMemcpy(dx, hx, hipMemcpyHostToDevice))

  res = 0.0
  call hipblasCheck(hipblasSasum(handle, n, dx, 1, c_loc(res)))
  call hipCheck(hipDeviceSynchronize())

  error = abs((res_exact - res) / res_exact)
  if (error > error_max) then
    write(*,*) "FAILED! error = ", error, " result = ", res
    call exit(1)
  end if

  call hipCheck(hipFree(dx))
  call hipblasCheck(hipblasDestroy(handle))

  write(*,*) "PASSED!"

end program hip_sasum

test/f2008/hipblas/dasum.f08 is the double-precision equivalent, and scasum.f08 and dzasum.f08 are the mixed real/complex forms named on the same convention as scnrm2.

Index of the largest element (iamax)#

hipblasI?amax returns the index of the element with the largest absolute value. The returned index is 1-based, so it can be used to subscript a Fortran array directly.

program hip_isamax
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas
  use hipfort_hipblas_enums

  implicit none

  ! The largest |x| sits at index 5, and hipBLAS returns a 1-based index.
  integer, parameter :: n = 8
  real(c_float), allocatable, dimension(:) :: hx
  integer(c_int), target :: res
  integer(c_int), parameter :: res_exact = 5
  real(c_float), pointer, dimension(:) :: dx
  type(c_ptr) :: handle = c_null_ptr

  write(*,"(a)",advance="no") "-- Running test 'isamax' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  allocate(hx(n))
  hx = 1.0
  hx(5) = 10.0

  call hipCheck(hipMalloc(dx, shape(hx)))
  call hipCheck(hipMemcpy(dx, hx, hipMemcpyHostToDevice))

  res = 0
  call hipblasCheck(hipblasIsamax(handle, n, dx, 1, c_loc(res)))
  call hipCheck(hipDeviceSynchronize())

  if (res /= res_exact) then
    write(*,*) "FAILED! wrong index: result = ", res, " expected ", res_exact
    call exit(1)
  end if

  call hipCheck(hipFree(dx))
  call hipblasCheck(hipblasDestroy(handle))

  write(*,*) "PASSED!"

end program hip_isamax

test/f2008/hipblas/idamax.f08, icamax.f08 and izamax.f08 cover the remaining precisions. The matching iamin routines have no hipBLAS example here; the rocBLAS examples cover them.

Level 2: matrix-vector operations#

Matrix-vector multiplication#

hipblas?gemv computes y := alpha * op(A) * x + beta * y. This example uses constant matrix and vector entries so the expected result is a constant vector and easy to check.

program hip_sgemv
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer :: m = 6
  integer :: n = 5
  integer :: i, j

  real, parameter :: alpha = 1.0
  real, parameter :: beta = 0.0

  type(c_ptr) :: handle = c_null_ptr

  real(kind=4), allocatable, target, dimension(:) :: a, x, y

  real(kind=4), pointer, dimension(:) :: da, dx, dy

  real :: error
  real, parameter :: error_max = 10*epsilon(error)

  allocate(x(n))
  allocate(y(m))
  allocate(a(m*n))

  a(:) = 1.0
  x(:) = 1.0
  y(:) = 1.0

  write(*,"(a)",advance="no") "-- Running test 'SGEMV' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  call hipCheck(hipMalloc(dx,source=x))
  call hipCheck(hipMalloc(dy,source=y))
  call hipCheck(hipMalloc(da,source=a))

  call hipCheck(hipblasSgemv(handle,HIPBLAS_OP_N,m,n,alpha,da,m,dx,1,beta,dy,1))

  call hipCheck(hipMemcpy(y, dy, hipMemcpyDeviceToHost))

  do i = 1,m
    error = abs(5.0 - y(i))
      if( error > 10*epsilon(error) )then
        write(*,*) "FAILED! Error bigger than max! Error = ", error, "y(i) = ", y(i)
        call exit(1)
     end if
  end do

  call hipblasCheck(hipblasDestroy(handle))

  call hipCheck(hipFree(da))
  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))

  deallocate(a,x,y)

  write(*,*) "PASSED!"

end program hip_sgemv

test/f2008/hipblas/dgemv.f08, cgemv.f08 and zgemv.f08 cover the remaining precisions.

Rank-1 update (ger)#

hipblas?ger computes A := alpha * x * y**T + A, adding the outer product of two vectors to a matrix in place. This program keeps the type(c_ptr) device pointers and explicit byte counts of the Fortran 2003 style even though it lives among the Fortran 2008 sources.

! HIPBLAS, i.e. cuBLAS and rocBLAS, assumes column-major matrix memory layouts.
! Hence, no matrix must be transposed when interfacing with Fortran.
program hip_sger
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none
  
  integer :: m = 6 
  integer :: n = 5
  real, parameter :: alpha = 2.0
  type(c_ptr) :: handle = c_null_ptr
  integer :: i, j
  
  real :: error
  real, parameter :: error_max = 10*epsilon(error)

  real(kind=4), allocatable, target, dimension(:,:) :: a
  real(kind=4), allocatable, target, dimension(:) :: x, y
 
  type(c_ptr) :: da = c_null_ptr, dx = c_null_ptr, dy = c_null_ptr
  
  integer(c_size_t) :: Nabytes, Nxbytes, Nybytes
  integer, parameter :: bytes_per_element = 4 !float precision
  
  Nxbytes = m * bytes_per_element
  Nybytes = n * bytes_per_element
  Nabytes = m * n * bytes_per_element
  
  allocate(x(m))
  allocate(y(n))
  allocate(a(m,n))
  
  do i = 1,m
    do j = 1,n
      a(i,j) = 1.0 
    end do
  end do
  
  do i = 1,m
    x(i) = 1.0
  end do

  do i = 1,n
    y(i) = 1.0
  end do
 

  write(*,"(a)",advance="no") "-- Running test 'SGER' (Fortran 2008 interfaces) - "
  
  call hipblasCheck(hipblasCreate(handle))
  
  call hipCheck(hipMalloc(dx,Nxbytes))  
  call hipCheck(hipMalloc(dy,Nybytes))
  call hipCheck(hipMalloc(da,Nabytes))
  
  !call hipCheck(hipblasSetMatrix(m,n,bytes_per_element,a,m,da,m))
  !call hipCheck(hipblasSetVector(m,bytes_per_element,x,1,dx,1))
  !call hipCheck(hipblasSetVector(n,bytes_per_element,y,1,dy,1))
  
  call hipCheck(hipMemcpy(da, c_loc(a), Nabytes, hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(dx, c_loc(x), Nxbytes, hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(dy, c_loc(y), Nybytes, hipMemcpyHostToDevice)) 

  call hipCheck(hipblasSger(handle,m,n,alpha,dx,1,dy,1,da,m))
  
  !call hipCheck(hipblasGetMatrix(m,n,bytes_per_element,da,m,a,m));
  call hipCheck(hipDeviceSynchronize())

  call hipCheck(hipMemcpy(c_loc(a), da, Nabytes, hipMemcpyDeviceToHost))

  !do i=1,m
  !  do j = 1,n
  !    write(*,*) a(i,j) 
  !  end do
  !end do

  do i = 1,m
  do j = 1,n
    error = abs(3.0 - a(i,j))
      if( error > 10*epsilon(error) )then
        write(*,*) "FAILED! Error bigger than max! Error = ", error, "a(i,j) = ", a(i,j)
        call exit(1)
     end if
  end do
  end do

  call hipblasCheck(hipblasDestroy(handle))

  call hipCheck(hipFree(da))
  call hipCheck(hipFree(dx))
  call hipCheck(hipFree(dy))
  
  deallocate(a)
  deallocate(x)
  deallocate(y)
  
  write(*,*) "PASSED!"
  
end program hip_sger

test/f2008/hipblas/dger.f08 is the double-precision equivalent. The rocBLAS suite additionally provides the conjugated and unconjugated complex forms, gerc and geru; see the rocBLAS examples.

Triangular solve#

hipblas?trsv solves A * x = b in place for a triangular matrix A, given the fill mode (HIPBLAS_FILL_MODE_LOWER or HIPBLAS_FILL_MODE_UPPER), the transpose operation and whether the diagonal is unit or not. dx holds b on entry and x on exit.

program hip_strsv
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer(kind(HIPBLAS_FILL_MODE_UPPER)), parameter :: uplo   = HIPBLAS_FILL_MODE_LOWER
  integer(kind(HIPBLAS_OP_N)),            parameter :: transA = HIPBLAS_OP_N
  integer(kind(HIPBLAS_DIAG_NON_UNIT)),   parameter :: diag   = HIPBLAS_DIAG_NON_UNIT
  integer, parameter :: m = 1024
  real(kind=4), allocatable, target, dimension(:,:) :: hA
  real(kind=4), allocatable, target, dimension(:)   :: hx
  real(kind=4), pointer, dimension(:,:) :: dA
  real(kind=4), pointer, dimension(:)   :: dx
  type(c_ptr) :: handle = c_null_ptr
  real(kind=4), parameter :: x_exact = 1.0
  real(kind=4) :: error
  real(kind=4), parameter :: error_max = 10*epsilon(error)
  integer :: i, j

  write(*,"(a)",advance="no") "-- Running test 'STRSV' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  allocate(hA(m,m), hx(m))

  hA = 0.0
  do j = 1, m
    do i = j, m
      hA(i,j) = 1.0
    end do
  end do
  do i = 1, m
    hx(i) = real(i)
  end do

  call hipCheck(hipMalloc(dA, source=hA))
  call hipCheck(hipMalloc(dx, source=hx))

  call hipblasCheck(hipblasStrsv(handle, uplo, transA, diag, m, dA, size(dA,1), dx, 1))
  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hx, dx, hipMemcpyDeviceToHost))

  do i = 1, m
    error = abs((x_exact - hx(i)) / x_exact)
    if (error > error_max) then
      write(*,*) "FAILED! Error bigger than max! Error = ", error, " hx(", i, ") = ", hx(i)
      call exit(1)
    end if
  end do

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dx))
  call hipblasCheck(hipblasDestroy(handle))
  deallocate(hA, hx)
  write(*,*) "PASSED!"

end program hip_strsv

test/f2008/hipblas/dtrsv.f08, ctrsv.f08 and ztrsv.f08 cover the remaining precisions.

Packed triangular solve#

hipblas?tpsv solves the same problem as trsv, but the triangular matrix is stored in packed form: only the referenced triangle is kept, in a single one-dimensional array, which halves the memory footprint for large matrices.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Copyright (c) 2020-2026 Advanced Micro Devices, Inc.
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in
! all copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


program hipblas_stpsv_test

  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer(kind(HIPBLAS_FILL_MODE_LOWER)), parameter :: uplo   = HIPBLAS_FILL_MODE_LOWER
  integer(kind(HIPBLAS_OP_N)),            parameter :: transA = HIPBLAS_OP_N
  integer(kind(HIPBLAS_DIAG_NON_UNIT)),   parameter :: diag   = HIPBLAS_DIAG_NON_UNIT

  integer, parameter :: n = 64

  ! Packed lower-triangular matrix (column-by-column). Every stored element of an
  ! all-ones lower-triangular matrix is 1, so the packing order is irrelevant here.
  real(c_float), allocatable, dimension(:) :: hAP, hx
  real(c_float), parameter :: x_exact = 1.0

  real(c_float), pointer, dimension(:) :: dAP => null(), dx => null()
  type(c_ptr) :: handle = c_null_ptr

  integer :: i
  real(c_float) :: error
  real(c_float), parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'STPSV' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  allocate(hAP(n*(n+1)/2), hx(n))

  hAP(:) = 1.0            ! all-ones lower-triangular L, packed
  do i = 1, n
    hx(i) = real(i)       ! b(i) = i  ->  exact solution of L x = b is x(i) = 1
  end do

  ! Allocate device memory (source= implies a blocking memcpy)
  call hipCheck(hipMalloc(dAP, source=hAP))
  call hipCheck(hipMalloc(dx,  source=hx))

  ! Packed triangular solve L x = b, in place. dAP is passed as a Fortran array
  ! (the generic array form) — this is the case that used to fail to compile
  ! because the packed matrix dummy was declared type(c_ptr) (SWDEV-485451).
  call hipblasCheck(hipblasStpsv(handle, uplo, transA, diag, n, dAP, dx, 1))

  call hipCheck(hipDeviceSynchronize())

  call hipCheck(hipMemcpy(hx, dx, hipMemcpyDeviceToHost))

  do i = 1, n
    error = abs((x_exact - hx(i))/x_exact)
    if( error > error_max )then
      write(*,*) "FAILED! Error bigger than max! Error = ", error, " hx(", i, ") = ", hx(i)
      call exit(1)
    end if
  end do

  call hipCheck(hipFree(dAP))
  call hipCheck(hipFree(dx))

  call hipblasCheck(hipblasDestroy(handle))

  deallocate(hAP, hx)

  write(*,*) "PASSED!"

end program hipblas_stpsv_test

This is the only packed-storage example; there is no double-precision or complex counterpart in the test suite.

Level 3: matrix-matrix operations#

Matrix-matrix multiplication#

hipblas?gemm computes C := alpha * op(A) * op(B) + beta * C. As with gemv, this example uses constant matrix entries so the exact result is a constant matrix.

program hip_dgemm

  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer(kind(HIPBLAS_OP_N)), parameter :: transa = HIPBLAS_OP_N, transb = HIPBLAS_OP_N;
  double precision, parameter ::  alpha = 1.1d0, beta = 0.9d0;

  integer, parameter ::  m = 1024, n = 1024, k = 1024;

  double precision, allocatable, dimension(:,:) :: ha, hb, hc, hc_exact

  double precision, pointer, dimension(:,:) :: da, db, dc
  type(c_ptr) :: handle = c_null_ptr

  integer :: i,j
  double precision :: error
  double precision, parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'DGEMM' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  allocate(ha(m,k))
  allocate(hb(k,n))
  allocate(hc(m,n))
  allocate(hc_exact(m,n))

  ! Use these constant matrices so the exact answer is also a
  ! constant matrix and therefore easy to check
  ha(:,:) = 1.d0
  hb(:,:) = 2.d0
  hc(:,:) = 3.d0
  hc_exact = alpha*k*2.d0 + beta*3.d0

  ! Allocate device memory
  call hipCheck(hipMalloc(da,source=ha)) ! implies (blocking) memcpy
  call hipCheck(hipMalloc(db,source=hb))
  call hipCheck(hipMalloc(dc,source=hc))

  call hipblasCheck(hipblasDgemm(handle,transa,transb,m,n,k,alpha,da,size(da,1),db,size(db,1),beta,dc,size(dc,1)))

  call hipCheck(hipDeviceSynchronize())

  ! Transfer data back to host memory
  call hipCheck(hipMemcpy(hc, dc, hipMemcpyDeviceToHost))

  do j = 1,n
    do i = 1,m
       error = abs((hc_exact(i,j) - hc(i,j))/hc_exact(i,j))
       if( error > error_max )then
          write(*,*) "FAILED! Error bigger than max! Error = ", error
          call exit(1)
       end if
    end do
  end do

  call hipCheck(hipFree(da))
  call hipCheck(hipFree(db))
  call hipCheck(hipFree(dc))

  call hipblasCheck(hipblasDestroy(handle))

  deallocate(ha,hb,hc,hc_exact)

  write(*,*) "PASSED!"

end program hip_dgemm

test/f2008/hipblas/sgemm.f08, cgemm.f08 and zgemm.f08 cover the remaining precisions.

Batched matrix multiplication#

hipblas?gemmBatched runs several independent gemm calls in one launch, with A, B and C each passed as a device array of device pointers. hipFORT’s Fortran 2008 pointer-convenience interfaces cannot be used for this argument shape: passing a plain array would only give hipBLAS a single pointer, not the array of device pointers the routine requires. This test therefore builds the pointer array with the same type(c_ptr) idiom used in the Fortran 2003 interfaces.

! Note: hipfort's Fortran 2008 pointer-convenience interfaces for
! hipblas?gemmBatched call c_loc() on a plain array and therefore pass a
! single pointer, not the array of device pointers the routine requires. This
! test consequently uses the type(c_ptr) array-of-pointers idiom (as in f2003);
! only the per-batch device allocation uses the f2008 hipMalloc(source=) style.

program hip_dgemm_batched

  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer(kind(HIPBLAS_OP_N)), parameter :: transa = HIPBLAS_OP_N, transb = HIPBLAS_OP_N
  double precision, parameter :: alpha = 1.1d0, beta = 0.9d0

  integer, parameter :: m = 512, n = 512, k = 512, batch_count = 4
  integer, parameter :: bytes_per_element = 8 ! double precision
  integer, parameter :: size_a = m*k, size_b = k*n, size_c = m*n
  integer(c_size_t), parameter :: Ncbytes = size_c*bytes_per_element
  integer :: lda, ldb, ldc, i, b

  double precision, allocatable, target, dimension(:,:) :: ha, hb, hc ! (elems, batch)
  double precision, allocatable, dimension(:) :: hc_exact ! one value per batch

  ! Per-batch device matrices as typed Fortran pointers (reused each iteration)
  double precision, pointer, dimension(:) :: da_b, db_b, dc_b
  ! Host arrays of device pointers (one device matrix per batch entry)
  type(c_ptr), target :: da(batch_count), db(batch_count), dc(batch_count)
  ! Device-resident pointer arrays passed to the batched routine
  type(c_ptr) :: da_p = c_null_ptr, db_p = c_null_ptr, dc_p = c_null_ptr
  type(c_ptr) :: handle = c_null_ptr

  double precision :: error
  double precision, parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'DGEMM_BATCHED' (Fortran 2008 interfaces) - "

  ! hipBLAS defaults to host pointer mode: no set-pointer-mode call needed
  call hipblasCheck(hipblasCreate(handle))

  lda = m; ldb = k; ldc = m

  allocate(ha(size_a,0:batch_count-1))
  allocate(hb(size_b,0:batch_count-1))
  allocate(hc(size_c,0:batch_count-1))
  allocate(hc_exact(0:batch_count-1))

  ! Constant matrices with a distinct per-batch value so the exact answer is
  ! a distinct constant per batch. A is held constant; B varies per batch.
  ha(:,:) = 1.d0
  do b = 0, batch_count-1
     hb(:,b) = dble(b+1)
     hc(:,b) = 3.d0
     hc_exact(b) = alpha*k*dble(b+1) + beta*3.d0
  end do

  ! Allocate one device matrix per batch entry (f2008 source= copy) and record
  ! its device address in the host pointer array
  do b = 1, batch_count
     call hipCheck(hipMalloc(da_b,source=ha(:,b-1)))
     call hipCheck(hipMalloc(db_b,source=hb(:,b-1)))
     call hipCheck(hipMalloc(dc_b,source=hc(:,b-1)))
     da(b) = c_loc(da_b(1))
     db(b) = c_loc(db_b(1))
     dc(b) = c_loc(dc_b(1))
  end do

  ! Allocate the device-resident pointer arrays and copy the host pointer
  ! arrays into them (hipBLAS requires the pointer array in device memory)
  call hipCheck(hipMalloc(da_p, int(batch_count,c_size_t)*c_sizeof(da(1))))
  call hipCheck(hipMalloc(db_p, int(batch_count,c_size_t)*c_sizeof(db(1))))
  call hipCheck(hipMalloc(dc_p, int(batch_count,c_size_t)*c_sizeof(dc(1))))
  call hipCheck(hipMemcpy(da_p, c_loc(da(1)), int(batch_count,c_size_t)*c_sizeof(da(1)), hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(db_p, c_loc(db(1)), int(batch_count,c_size_t)*c_sizeof(db(1)), hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(dc_p, c_loc(dc(1)), int(batch_count,c_size_t)*c_sizeof(dc(1)), hipMemcpyHostToDevice))

  call hipblasCheck(hipblasDgemmBatched(handle,transa,transb,m,n,k, &
       alpha,da_p,lda,db_p,ldb,beta,dc_p,ldc,batch_count))

  call hipCheck(hipDeviceSynchronize())

  ! Transfer each batch result back to host memory
  do b = 1, batch_count
     call hipCheck(hipMemcpy(c_loc(hc(1,b-1)), dc(b), Ncbytes, hipMemcpyDeviceToHost))
  end do

  do b = 0, batch_count-1
     do i = 1, size_c
        error = abs((hc_exact(b) - hc(i,b))/hc_exact(b))
        if( error > error_max )then
           write(*,*) "FAILED! Error bigger than max! batch = ", b, " error = ", error
           call exit(1)
        end if
     end do
  end do

  do b = 1, batch_count
     call hipCheck(hipFree(da(b)))
     call hipCheck(hipFree(db(b)))
     call hipCheck(hipFree(dc(b)))
  end do
  call hipCheck(hipFree(da_p))
  call hipCheck(hipFree(db_p))
  call hipCheck(hipFree(dc_p))

  call hipblasCheck(hipblasDestroy(handle))

  deallocate(ha)
  deallocate(hb)
  deallocate(hc)
  deallocate(hc_exact)

  write(*,*) "PASSED!"

end program hip_dgemm_batched

test/f2008/hipblas/sgemm_batched.f08, cgemm_batched.f08 and zgemm_batched.f08 cover the remaining precisions.

Strided-batched matrix multiplication#

hipblas?gemmStridedBatched also runs several gemm calls in one launch, but A, B and C are contiguous device arrays with a fixed stride between the start of each batch’s matrix, instead of an array of pointers. This is usually simpler to set up than the batched form when the batches are already laid out contiguously in memory.

program hip_dgemm_strided_batched

  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer(kind(HIPBLAS_OP_N)), parameter :: transa = HIPBLAS_OP_N, transb = HIPBLAS_OP_N
  double precision, parameter :: alpha = 1.1d0, beta = 0.9d0

  integer, parameter :: m = 512, n = 512, k = 512, batch_count = 4
  integer :: lda, ldb, ldc, i, b
  integer(c_int64_t) :: stride_a, stride_b, stride_c
  integer :: size_a, size_b, size_c

  double precision, allocatable, dimension(:) :: ha, hb, hc
  double precision, allocatable, dimension(:) :: hc_exact ! one value per batch

  double precision, pointer, dimension(:) :: da, db, dc
  type(c_ptr) :: handle = c_null_ptr

  double precision :: error
  double precision, parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'DGEMM_STRIDED_BATCHED' (Fortran 2008 interfaces) - "

  ! hipBLAS defaults to host pointer mode: no set-pointer-mode call needed
  call hipblasCheck(hipblasCreate(handle))

  lda = m; ldb = k; ldc = m
  stride_a = int(lda,c_int64_t)*k
  stride_b = int(ldb,c_int64_t)*n
  stride_c = int(ldc,c_int64_t)*n
  size_a = int(stride_a)*batch_count
  size_b = int(stride_b)*batch_count
  size_c = int(stride_c)*batch_count

  allocate(ha(size_a))
  allocate(hb(size_b))
  allocate(hc(size_c))
  allocate(hc_exact(0:batch_count-1))

  ! Constant matrices with a distinct per-batch value so the exact answer is
  ! a distinct constant per batch. A is held constant; B varies per batch.
  ha(:) = 1.d0
  do b = 0, batch_count-1
     hb(b*stride_b+1 : (b+1)*stride_b) = dble(b+1)
     hc(b*stride_c+1 : (b+1)*stride_c) = 3.d0
     hc_exact(b) = alpha*k*dble(b+1) + beta*3.d0
  end do

  ! Allocate device memory
  call hipCheck(hipMalloc(da,source=ha)) ! implies (blocking) memcpy
  call hipCheck(hipMalloc(db,source=hb))
  call hipCheck(hipMalloc(dc,source=hc))

  call hipblasCheck(hipblasDgemmStridedBatched(handle,transa,transb,m,n,k, &
       alpha,da,lda,stride_a,db,ldb,stride_b,beta,dc,ldc,stride_c,batch_count))

  call hipCheck(hipDeviceSynchronize())

  ! Transfer data back to host memory
  call hipCheck(hipMemcpy(hc, dc, hipMemcpyDeviceToHost))

  do b = 0, batch_count-1
     do i = 1, int(stride_c)
        error = abs((hc_exact(b) - hc(b*stride_c+i))/hc_exact(b))
        if( error > error_max )then
           write(*,*) "FAILED! Error bigger than max! batch = ", b, " error = ", error
           call exit(1)
        end if
     end do
  end do

  call hipCheck(hipFree(da))
  call hipCheck(hipFree(db))
  call hipCheck(hipFree(dc))

  call hipblasCheck(hipblasDestroy(handle))

  deallocate(ha)
  deallocate(hb)
  deallocate(hc)
  deallocate(hc_exact)

  write(*,*) "PASSED!"

end program hip_dgemm_strided_batched

test/f2008/hipblas/sgemm_strided_batched.f08, cgemm_strided_batched.f08 and zgemm_strided_batched.f08 cover the remaining precisions.

Triangular solve with multiple right-hand sides#

hipblas?trsm solves A * X = alpha * B in place for a triangular matrix A and a matrix of right-hand sides B, which is the gemm-like counterpart of trsv. dB holds B on entry and the solution X on exit.

program hip_dtrsm
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas

  implicit none

  integer(kind(HIPBLAS_SIDE_LEFT)),       parameter :: side   = HIPBLAS_SIDE_LEFT
  integer(kind(HIPBLAS_FILL_MODE_LOWER)), parameter :: uplo   = HIPBLAS_FILL_MODE_LOWER
  integer(kind(HIPBLAS_OP_N)),            parameter :: transA = HIPBLAS_OP_N
  integer(kind(HIPBLAS_DIAG_NON_UNIT)),   parameter :: diag   = HIPBLAS_DIAG_NON_UNIT
  integer, parameter :: m = 1024, n = 1024
  double precision, parameter :: alpha = 2.d0
  double precision, allocatable, target, dimension(:,:) :: hA, hB
  double precision, pointer, dimension(:,:) :: dA, dB
  type(c_ptr) :: handle = c_null_ptr
  double precision, parameter :: x_exact = 1.d0
  double precision :: error
  double precision, parameter :: error_max = 10*epsilon(error)
  integer :: i, j

  write(*,"(a)",advance="no") "-- Running test 'DTRSM' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  allocate(hA(m,m), hB(m,n))

  hA = 0.d0
  do j = 1, m
    do i = j, m
      hA(i,j) = 1.d0
    end do
  end do
  do j = 1, n
    do i = 1, m
      hB(i,j) = dble(i) / 2.d0
    end do
  end do

  call hipCheck(hipMalloc(dA, source=hA))
  call hipCheck(hipMalloc(dB, source=hB))

  call hipblasCheck(hipblasDtrsm(handle, side, uplo, transA, diag, m, n, alpha, dA, size(dA,1), dB, size(dB,1)))
  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hB, dB, hipMemcpyDeviceToHost))

  do j = 1, n
    do i = 1, m
      error = abs((x_exact - hB(i,j)) / x_exact)
      if (error > error_max) then
        write(*,*) "FAILED! Error bigger than max! Error = ", error, " hB(", i, ",", j, ") = ", hB(i,j)
        call exit(1)
      end if
    end do
  end do

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipblasCheck(hipblasDestroy(handle))
  deallocate(hA, hB)
  write(*,*) "PASSED!"

end program hip_dtrsm

test/f2008/hipblas/strsm.f08, ctrsm.f08 and ztrsm.f08 cover the remaining precisions.

Triangular matrix multiplication#

hipblas?trmm computes C := alpha * op(A) * B (or the mirrored right-hand form), where A is triangular and only the triangle chosen by the fill mode is referenced.

program hip_strmm
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas
  use hipfort_hipblas_enums

  implicit none

  ! C := alpha*op(A)*B with A upper triangular (out-of-place trmm).
  ! A = [1 2; 0 3] and B = I, so C = A.
  integer, parameter :: ld = 2
  real(c_float), parameter :: alpha = 1.0

  real(c_float) :: hA(ld,ld) = reshape([1.0, 0.0, 2.0, 3.0], [ld,ld])
  real(c_float) :: hB(ld,ld) = reshape([1.0, 0.0, 0.0, 1.0], [ld,ld])
  real(c_float) :: hC(ld,ld)
  real(c_float) :: expected(ld,ld) = reshape([1.0, 0.0, 2.0, 3.0], [ld,ld])
  real(c_float), pointer, dimension(:,:) :: dA, dB, dC
  type(c_ptr) :: handle = c_null_ptr
  integer :: i, j
  real(c_float) :: error
  real(c_float), parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'strmm' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  hC = 0.0

  call hipCheck(hipMalloc(dA, shape(hA)))
  call hipCheck(hipMemcpy(dA, hA, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dB, shape(hB)))
  call hipCheck(hipMemcpy(dB, hB, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dC, shape(hC)))
  call hipCheck(hipMemcpy(dC, hC, hipMemcpyHostToDevice))

  call hipblasCheck(hipblasStrmm(handle, HIPBLAS_SIDE_LEFT, HIPBLAS_FILL_MODE_UPPER, HIPBLAS_OP_N, &
       HIPBLAS_DIAG_NON_UNIT, ld, ld, alpha, dA, size(dA,1), dB, size(dB,1), dC, size(dA,1)))
  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hC, dC, hipMemcpyDeviceToHost))

  do j = 1, ld
    do i = 1, ld
      error = abs(expected(i,j) - hC(i,j))
      if (error > error_max) then
        write(*,*) "FAILED! error = ", error, " at ", i, j
        call exit(1)
      end if
    end do
  end do

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dC))
  call hipblasCheck(hipblasDestroy(handle))

  write(*,*) "PASSED!"

end program hip_strmm

test/f2008/hipblas/dtrmm.f08, ctrmm.f08 and ztrmm.f08 cover the remaining precisions.

Rank-k update (syrk)#

hipblas?syrk computes C := alpha * op(A) * op(A)**T + beta * C, where C is symmetric and only the triangle chosen by the fill mode is referenced.

program hip_ssyrk
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas
  use hipfort_hipblas_enums

  implicit none

  ! C := alpha*A*A**T + beta*C, upper triangle. A = [1 0; 2 3] gives
  ! A*A**T = [1 2; 2 13], so the referenced upper triangle is 1, 2, 13.
  integer, parameter :: ld = 2
  real(c_float), parameter :: alpha = 1.0, beta = 0.0

  real(c_float) :: hA(ld,ld) = reshape([1.0, 2.0, 0.0, 3.0], [ld,ld])
  real(c_float) :: hC(ld,ld)
  real(c_float) :: expected(ld,ld) = reshape([1.0, 0.0, 2.0, 13.0], [ld,ld])
  real(c_float), pointer, dimension(:,:) :: dA, dC
  type(c_ptr) :: handle = c_null_ptr
  integer :: i, j
  real(c_float) :: error
  real(c_float), parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'ssyrk' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  hC = 0.0

  call hipCheck(hipMalloc(dA, shape(hA)))
  call hipCheck(hipMemcpy(dA, hA, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dC, shape(hC)))
  call hipCheck(hipMemcpy(dC, hC, hipMemcpyHostToDevice))

  call hipblasCheck(hipblasSsyrk(handle, HIPBLAS_FILL_MODE_UPPER, HIPBLAS_OP_N, ld, ld, &
       alpha, dA, size(dA,1), beta, dC, size(dA,1)))
  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hC, dC, hipMemcpyDeviceToHost))

  do j = 1, ld
    do i = 1, j
      error = abs(expected(i,j) - hC(i,j))
      if (error > error_max) then
        write(*,*) "FAILED! error = ", error, " at ", i, j
        call exit(1)
      end if
    end do
  end do

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dC))
  call hipblasCheck(hipblasDestroy(handle))

  write(*,*) "PASSED!"

end program hip_ssyrk

test/f2008/hipblas/dsyrk.f08, csyrk.f08 and zsyrk.f08 cover the remaining precisions. The Hermitian form, herk, has no hipBLAS example here; the rocBLAS examples cover it.

Symmetric matrix product (symm)#

hipblas?symm computes C := alpha * A * B + beta * C with A symmetric, or the mirrored right-hand form selected by the side argument. As with syrk, only one triangle of A is referenced.

program hip_ssymm
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas
  use hipfort_hipblas_enums

  implicit none

  ! C := alpha*A*B + beta*C with A symmetric (upper triangle referenced).
  ! A = [1 2; 2 3] and B = I, so C = A.
  integer, parameter :: ld = 2
  real(c_float), parameter :: alpha = 1.0, beta = 0.0

  real(c_float) :: hA(ld,ld) = reshape([1.0, 2.0, 2.0, 3.0], [ld,ld])
  real(c_float) :: hB(ld,ld) = reshape([1.0, 0.0, 0.0, 1.0], [ld,ld])
  real(c_float) :: hC(ld,ld)
  real(c_float) :: expected(ld,ld) = reshape([1.0, 2.0, 2.0, 3.0], [ld,ld])
  real(c_float), pointer, dimension(:,:) :: dA, dB, dC
  type(c_ptr) :: handle = c_null_ptr
  integer :: i, j
  real(c_float) :: error
  real(c_float), parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'ssymm' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  hC = 0.0

  call hipCheck(hipMalloc(dA, shape(hA)))
  call hipCheck(hipMemcpy(dA, hA, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dB, shape(hB)))
  call hipCheck(hipMemcpy(dB, hB, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dC, shape(hC)))
  call hipCheck(hipMemcpy(dC, hC, hipMemcpyHostToDevice))

  call hipblasCheck(hipblasSsymm(handle, HIPBLAS_SIDE_LEFT, HIPBLAS_FILL_MODE_UPPER, ld, ld, &
       alpha, dA, size(dA,1), dB, size(dB,1), beta, dC, size(dA,1)))
  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hC, dC, hipMemcpyDeviceToHost))

  do j = 1, ld
    do i = 1, ld
      error = abs(expected(i,j) - hC(i,j))
      if (error > error_max) then
        write(*,*) "FAILED! error = ", error, " at ", i, j
        call exit(1)
      end if
    end do
  end do

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dC))
  call hipblasCheck(hipblasDestroy(handle))

  write(*,*) "PASSED!"

end program hip_ssymm

test/f2008/hipblas/dsymm.f08, csymm.f08 and zsymm.f08 cover the remaining precisions. The Hermitian form, hemm, likewise appears only among the rocBLAS programs.

Matrix addition and transposition (geam)#

hipblas?geam computes C := alpha * op(A) + beta * op(B). Because each operand has its own transpose flag and either scalar may be zero, the same routine also serves as an out-of-place transpose or a scaled copy.

program hip_sgeam
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas
  use hipfort_hipblas_enums

  implicit none

  ! C := alpha*op(A) + beta*op(B) with alpha = beta = 1 and no
  ! transposition, so C is the elementwise sum A + B.
  integer, parameter :: ld = 2
  real(c_float), parameter :: alpha = 1.0, beta = 1.0

  real(c_float) :: hA(ld,ld) = reshape([1.0, 2.0, 3.0, 4.0], [ld,ld])
  real(c_float) :: hB(ld,ld) = reshape([10.0, 20.0, 30.0, 40.0], [ld,ld])
  real(c_float) :: hC(ld,ld)
  real(c_float) :: expected(ld,ld) = reshape([11.0, 22.0, 33.0, 44.0], [ld,ld])
  real(c_float), pointer, dimension(:,:) :: dA, dB, dC
  type(c_ptr) :: handle = c_null_ptr
  integer :: i, j
  real(c_float) :: error
  real(c_float), parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'sgeam' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  hC = 0.0

  call hipCheck(hipMalloc(dA, shape(hA)))
  call hipCheck(hipMemcpy(dA, hA, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dB, shape(hB)))
  call hipCheck(hipMemcpy(dB, hB, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dC, shape(hC)))
  call hipCheck(hipMemcpy(dC, hC, hipMemcpyHostToDevice))

  call hipblasCheck(hipblasSgeam(handle, HIPBLAS_OP_N, HIPBLAS_OP_N, ld, ld, &
       alpha, dA, size(dA,1), beta, dB, size(dB,1), dC, size(dA,1)))
  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hC, dC, hipMemcpyDeviceToHost))

  do j = 1, ld
    do i = 1, ld
      error = abs(expected(i,j) - hC(i,j))
      if (error > error_max) then
        write(*,*) "FAILED! error = ", error, " at ", i, j
        call exit(1)
      end if
    end do
  end do

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dC))
  call hipblasCheck(hipblasDestroy(handle))

  write(*,*) "PASSED!"

end program hip_sgeam

test/f2008/hipblas/dgeam.f08, cgeam.f08 and zgeam.f08 cover the remaining precisions.

Extended-precision matrix multiplication (GemmEx)#

hipblasGemmEx computes D := alpha * op(A) * op(B) + beta * C with the type of every buffer, and the type used for the arithmetic, given explicitly as arguments rather than fixed by the routine name. That makes it the entry point for mixed precision work.

program hip_gemmex
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipblas
  use hipfort_hipblas_enums

  implicit none

  ! C := alpha*op(A)*op(B) + beta*C through the extended-precision entry point,
  ! with all buffers f32 and an f32 compute type. A = I and beta = 0, so C = B.
  integer, parameter :: ld = 2
  real(c_float), target :: alpha = 1.0, beta = 0.0

  real(c_float) :: hA(ld,ld) = reshape([1.0, 0.0, 0.0, 1.0], [ld,ld])
  real(c_float) :: hB(ld,ld) = reshape([1.0, 2.0, 3.0, 4.0], [ld,ld])
  real(c_float) :: hC(ld,ld)
  real(c_float) :: expected(ld,ld) = reshape([1.0, 2.0, 3.0, 4.0], [ld,ld])
  real(c_float), pointer, dimension(:,:) :: dA, dB, dC
  type(c_ptr) :: handle = c_null_ptr
  integer :: i, j
  real(c_float) :: error
  real(c_float), parameter :: error_max = 10*epsilon(error)

  write(*,"(a)",advance="no") "-- Running test 'gemmex' (Fortran 2008 interfaces) - "

  call hipblasCheck(hipblasCreate(handle))

  hC = 0.0

  call hipCheck(hipMalloc(dA, shape(hA)))
  call hipCheck(hipMemcpy(dA, hA, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dB, shape(hB)))
  call hipCheck(hipMemcpy(dB, hB, hipMemcpyHostToDevice))
  call hipCheck(hipMalloc(dC, shape(hC)))
  call hipCheck(hipMemcpy(dC, hC, hipMemcpyHostToDevice))

  call hipblasCheck(hipblasGemmEx(handle, HIPBLAS_OP_N, HIPBLAS_OP_N, ld, ld, ld, &
       c_loc(alpha), c_loc(dA(1,1)), HIP_R_32F, ld, &
       c_loc(dB(1,1)), HIP_R_32F, ld, c_loc(beta), &
       c_loc(dC(1,1)), HIP_R_32F, ld, HIPBLAS_COMPUTE_32F, HIPBLAS_GEMM_DEFAULT))
  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hC, dC, hipMemcpyDeviceToHost))

  do j = 1, ld
    do i = 1, ld
      error = abs(expected(i,j) - hC(i,j))
      if (error > error_max) then
        write(*,*) "FAILED! error = ", error, " at ", i, j
        call exit(1)
      end if
    end do
  end do

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dC))
  call hipblasCheck(hipblasDestroy(handle))

  write(*,*) "PASSED!"

end program hip_gemmex

Because the buffer types are runtime arguments, there is a single GemmEx program rather than one per precision. Note the source file is gemmex.f08, without the underscore used by the rocBLAS equivalent gemm_ex.f08.

Contents