rocSOLVER examples

rocSOLVER examples#

rocSOLVER is the AMD implementation of LAPACK for AMD GPUs. hipFORT exposes it through the hipfort_rocsolver module, which mirrors the rocSOLVER C API one to one. rocSOLVER is built on rocBLAS and reuses its handle type, so every program also uses the hipfort_rocblas module for rocblas_create_handle and the rocblas_* enumerators.

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/rocsolver 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/rocsolver.

Each example is provided in the four LAPACK precisions where the routine has them: 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 rocsolver_ prefix letter.

Conventions#

rocSOLVER follows a small number of conventions that recur in every example:

  • Column-major storage. rocSOLVER matrices are column-major, which matches Fortran’s native array layout, so a Fortran 2-D array maps directly onto a rocSOLVER matrix with leading dimension lda = size(A, 1).

  • The info output lives in device memory. rocSOLVER writes the factorization status info to a device pointer, so the hipFORT binding types that argument as type(c_ptr). Back it with a device allocation and pass c_loc(dInfo); passing a host scalar faults on the GPU. For the batched routines info is an array of batch_count integers on the device.

  • Pivots and scalar factors are device arrays. Arguments such as ipiv (pivot indices) and tau (Householder scalars) are outputs written on the device and are passed as device buffers.

  • rocBLAS enumerators select variants. rocblas_operation_none / rocblas_operation_transpose choose whether a routine works on A or A**T; rocblas_fill_upper / rocblas_fill_lower choose the stored triangle; rocblas_evect_* and rocblas_svect_* choose whether vectors are computed.

  • Every call returns a status code. The examples wrap rocSOLVER calls in rocsolverCheck and HIP calls in hipCheck from the hipfort_check module, both of which abort on failure.

Building an example#

The examples need the rocsolver, rocblas, and hip hipFORT components:

find_package(hipfort REQUIRED COMPONENTS hip rocblas rocsolver)

add_executable(my_solver rocsolver_dgetrf.f08)
target_link_libraries(my_solver PRIVATE hipfort::rocsolver hipfort::rocblas hipfort::hip)

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

LU factorization and solve#

getrf computes the LU factorization A = P*L*U with partial pivoting, writing the factors in place over A and the pivot indices into ipiv. The example factorizes a matrix and reconstructs L*U to confirm the result.

!!!!!!!!!!!!!/
! dgetrf example (double-precision LU factorization)
! see: https:!www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga0019443faea08275ca60a734d0593e60.html
!
! NOTE: rocSOLVER writes the `info` output to DEVICE memory. The hipfort
! binding types the info argument as a device pointer (type(c_ptr)), so it must be backed
! by a device allocation (dInfo below) and passed as c_loc(dInfo); passing a
! host scalar faults on the GPU.
!!!!!!!!!!!!!!/
!
program dgetrf
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocsolver

  implicit none
  integer :: i, j ! indices for iterating over results

  ! Define our input data (column-major)
  real(c_double) :: hA(3,3) = reshape((/12, 6, -4, -51, 167, 24, 4, -68, -41/), (/3, 3/))
  ! Reference: packed LU factors (L below diagonal, U on/above) from rocSOLVER
  real(c_double) :: hResult(3,3) = reshape((/&
    12.0000000000000000d0,   0.500000000000000000d0,  -0.333333333333333315d0,&
   -51.0000000000000000d0, 192.500000000000000d0,      0.363636363636363688d-01,&
     4.00000000000000000d0, -70.0000000000000000d0,   -37.1212121212121176d0/), shape(hResult), order=(/2,1/))
  integer(c_int) :: hIpiv_ref(3) = (/1, 2, 3/)
  integer(c_int), parameter :: M = 3
  integer(c_int), parameter :: N = 3
  integer(c_int), parameter :: lda = 3

  integer(c_int) :: hIpiv(3) ! CPU buffer for pivot indices

  real(c_double), pointer :: dA(:,:)    ! GPU buffer for A
  integer(c_int), pointer :: dIpiv(:)   ! GPU buffer for pivots
  integer(c_int), pointer :: dInfo(:)   ! GPU buffer for info (rocSOLVER writes to device)

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: error
  real(c_double), parameter :: error_max = 10 * epsilon(error_max)
  !
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dgetrf' (Fortran 2008 interfaces) - "

  ! Allocate device-side memory & copy memory from host to device
  call hipCheck(hipMalloc(dA,    source=hA))
  call hipCheck(hipMalloc(dIpiv, mold=hIpiv))
  call hipCheck(hipMalloc(dInfo, 1))

  ! Create rocBLAS handle
  call hipCheck(rocblas_create_handle(handle))

  ! Compute the LU factorization on the device
  call hipCheck(rocsolver_dgetrf(handle, M, N, dA, lda, dIpiv, c_loc(dInfo)))

  ! Copy result from device to host
  call hipCheck(hipMemcpy(hA,    dA,    hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hIpiv, dIpiv, hipMemcpyDeviceToHost))

  ! Check factor values
  do j = 1,size(hA,2)
    do i = 1,size(hA,1)
        error = abs(hA(i,j) - hResult(i,j)) / max(abs(hResult(i,j)), 1.0_c_double)
        if(error .gt. error_max) then
            write(*,*) "FAILED! Error bigger than max! Error = ", error, " hA(", i, ",", j, ") = ", hA(i,j)
            call exit
        end if
    end do
  end do

  ! Check pivots
  do i = 1,3
    if(hIpiv(i) .ne. hIpiv_ref(i)) then
        write(*,*) "FAILED! Pivot mismatch at ", i, " got ", hIpiv(i), " expected ", hIpiv_ref(i)
        call exit
    end if
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dIpiv))
  call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dgetrf

getrs uses the factors and pivots from getrf to solve A*X = B. The example picks a known solution x, forms b = A*x, factorizes, solves, and checks that the recovered X matches x.

!!!!!!!!!!!!!/
! dgetrs example (double-precision LU solve)
! see: https:!www.netlib.org/lapack/explore-html/d6/d49/group__double_g_ecomputational_gafa35ce1d7865b80563bbed6317050ad7.html
!
! Self-verifying: pick a known solution x, form b = A*x, factorize A with
! getrf, solve A*X = b with getrs, and confirm X recovers x.
!
! NOTE: getrf writes `info` to DEVICE memory, so info is backed by a device
! allocation and passed as c_loc(dInfo); passing a host scalar faults on the GPU.
!!!!!!!!!!!!!!/
!
program dgetrs
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocblas_enums
  use hipfort_rocsolver

  implicit none
  integer :: i ! index for iterating over results

  integer(c_int), parameter :: N = 3
  integer(c_int), parameter :: lda = 3
  integer(c_int), parameter :: ldb = 3
  integer(c_int), parameter :: nrhs = 1

  ! Nonsingular input matrix (column-major)
  real(c_double) :: hA(3,3) = reshape((/4, 6, -4, -51, 167, 24, 4, -68, -41/), (/3, 3/))
  real(c_double) :: hX(3,1) = reshape((/1, 2, 3/), (/3, 1/))  ! known solution
  real(c_double) :: hB(3,1)                                    ! RHS = A*x

  integer(c_int) :: hIpiv(3)

  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dB(:,:)
  integer(c_int), pointer :: dIpiv(:)
  integer(c_int), pointer :: dInfo(:)

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: error
  real(c_double), parameter :: error_max = 100 * epsilon(error_max)
  !
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dgetrs' (Fortran 2008 interfaces) - "

  ! Build a consistent RHS so that A*x = b
  hB = matmul(hA, hX)

  ! Allocate device-side memory & copy memory from host to device
  call hipCheck(hipMalloc(dA,    source=hA))
  call hipCheck(hipMalloc(dB,    source=hB))
  call hipCheck(hipMalloc(dIpiv, mold=hIpiv))
  call hipCheck(hipMalloc(dInfo, 1))

  ! Create rocBLAS handle
  call hipCheck(rocblas_create_handle(handle))

  ! Factorize, then solve A*X = B in place (B is overwritten with the solution)
  call hipCheck(rocsolver_dgetrf(handle, N, N, dA, lda, dIpiv, c_loc(dInfo)))
  call hipCheck(rocsolver_dgetrs(handle, rocblas_operation_none, N, nrhs, dA, lda, dIpiv, dB, ldb))

  ! Copy result from device to host
  call hipCheck(hipMemcpy(hB, dB, hipMemcpyDeviceToHost))

  ! Verify the recovered solution matches x
  do i = 1,N
    error = abs(hB(i,1) - hX(i,1)) / max(abs(hX(i,1)), 1.0_c_double)
    if(error .gt. error_max) then
        write(*,*) "FAILED! Error bigger than max! Error = ", error, " X(", i, ") = ", hB(i,1)
        call exit
    end if
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dIpiv))
  call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dgetrs

The getrf_npvt variant factorizes without pivoting (valid when no row interchanges are needed, as for a diagonally dominant matrix), and getf2 is the unblocked kernel with the same interface. See test/f2008/rocsolver/rocsolver_dgetrf_npvt.f08 and rocsolver_dgetf2.f08.

Batched LU#

rocSOLVER factorizes many matrices with one call in two forms. The array-of-pointers form, getrf_batched, takes A as a device array of per-matrix device pointers, and ipiv and info as device arrays indexed by batch.

!!!!!!!!!!!!!!
! dgetrf_batched example (batched LU, array-of-pointers form)
! Exercises the "array of device pointers" argument class: A is a device array
! of per-batch matrix pointers (not the strided form). The input is diagonally
! dominant (no pivoting), so each batch's L*U == A is checked by reconstruction,
! plus info == 0 per batch. rocSOLVER writes info to DEVICE memory.
!!!!!!!!!!!!!!
!
program dgetrf_batched
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocsolver
  implicit none
  integer :: i, j, k, b
  integer(c_int), parameter :: M = 3, N = 3, lda = 3, batch = 2
  real(c_double), target :: hA(3,3)  = reshape((/4, 1, 1,  1, 4, 1,  1, 1, 4/), (/3,3/))
  real(c_double), target :: hLU(3,3)
  integer(c_int) :: hInfo(batch)
  integer(c_int), target :: hInfo_t(batch)
  type(c_ptr) :: d1 = c_null_ptr, d2 = c_null_ptr
  type(c_ptr), target :: hptrs(batch)
  type(c_ptr) :: dA_ptrs = c_null_ptr, dIpiv = c_null_ptr, dInfo = c_null_ptr, handle = c_null_ptr
  integer(c_size_t) :: matbytes, psize
  real(c_double) :: lij, recon, error
  real(c_double), parameter :: error_max = 1.0d-10

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

  matbytes = int(M,c_size_t) * int(N,c_size_t) * 8
  psize    = c_sizeof(d1)

  call hipCheck(hipMalloc(d1, matbytes))
  call hipCheck(hipMalloc(d2, matbytes))
  call hipCheck(hipMemcpy(d1, c_loc(hA(1,1)), matbytes, hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(d2, c_loc(hA(1,1)), matbytes, hipMemcpyHostToDevice))

  hptrs(1) = d1
  hptrs(2) = d2
  call hipCheck(hipMalloc(dA_ptrs, int(batch,c_size_t) * psize))
  call hipCheck(hipMemcpy(dA_ptrs, c_loc(hptrs(1)), int(batch,c_size_t) * psize, hipMemcpyHostToDevice))

  call hipCheck(hipMalloc(dIpiv, int(N,c_size_t) * int(batch,c_size_t) * 4))
  call hipCheck(hipMalloc(dInfo, int(batch,c_size_t) * 4))

  call hipCheck(rocblas_create_handle(handle))
  call hipCheck(rocsolver_dgetrf_batched(handle, M, N, dA_ptrs, lda, dIpiv, int(N,c_int64_t), dInfo, batch))

  call hipCheck(hipMemcpy(c_loc(hInfo_t(1)), dInfo, int(batch,c_size_t)*4, hipMemcpyDeviceToHost))
  hInfo = hInfo_t
  do b = 1, batch
     if (hInfo(b) /= 0) then
        write(*,*) "FAILED! info(", b, ") = ", hInfo(b)
        call exit(1)
     end if
  end do

  ! Reconstruct L*U == A for each batch matrix.
  do b = 1, batch
     if (b == 1) then
        call hipCheck(hipMemcpy(c_loc(hLU(1,1)), d1, matbytes, hipMemcpyDeviceToHost))
     else
        call hipCheck(hipMemcpy(c_loc(hLU(1,1)), d2, matbytes, hipMemcpyDeviceToHost))
     end if
     do j = 1, N
        do i = 1, M
           recon = 0.0d0
           do k = 1, min(i,j)
              if (k == i) then; lij = 1.0d0; else; lij = hLU(i,k); end if
              recon = recon + lij * hLU(k,j)
           end do
           error = abs(recon - hA(i,j))
           if (error > error_max) then
              write(*,*) "FAILED! batch ", b, " (L*U)(", i, ",", j, ") = ", recon
              call exit(1)
           end if
        end do
     end do
  end do

  call hipCheck(hipFree(d1)); call hipCheck(hipFree(d2)); call hipCheck(hipFree(dA_ptrs))
  call hipCheck(hipFree(dIpiv)); call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle)); call hipCheck(hipDeviceReset())

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

The strided-batched form, getrf_strided_batched, instead stores the matrices contiguously in one device buffer and locates each by a fixed stride, which avoids building a pointer array. See test/f2008/rocsolver/rocsolver_dgetrf_strided_batched.f08.

The 64-bit integer API#

The _64 routines accept 64-bit problem dimensions, pivots, and info for problems that exceed the 32-bit range. getrf_64 performs the same factorization as getrf with integer(c_int64_t) dimensions; the device buffers, including the int64 ipiv and info, are passed as type(c_ptr) because the _64 routines have no native-array overloads.

!!!!!!!!!!!!!!
! dgetrf_64 example (LU with the 64-bit integer API)
! Same math as dgetrf; dimensions, pivots and info are 64-bit. Device buffers
! (int64 ipiv/info included) are passed as type(c_ptr) - the _64 routines have
! no native-array overloads. info lives on the DEVICE. The input is diagonally
! dominant so no pivoting occurs and L*U == A is checked by reconstruction
! (layout-agnostic).
!!!!!!!!!!!!!!
!
program dgetrf_64
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocsolver
  implicit none
  integer(c_int64_t), parameter :: M = 3, N = 3, lda = 3
  real(c_double), target :: hA(3,3)  = reshape((/4, 1, 1,  1, 4, 1,  1, 1, 4/), (/3, 3/))
  real(c_double), target :: hLU(3,3)
  integer(c_size_t) :: size_A = 9
  type(c_ptr) :: dA, dIpiv, dInfo, handle
  integer :: i, j, k
  real(c_double) :: lij, recon, error
  real(c_double), parameter :: error_max = 1.0d-10
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dgetrf_64' (Fortran 2008 interfaces) - "
  call hipCheck(hipMalloc(dA,    size_A * 8))
  call hipCheck(hipMalloc(dIpiv, 3_c_size_t * 8))   ! int64 pivots
  call hipCheck(hipMalloc(dInfo, 8_c_size_t))       ! int64 info
  call hipCheck(rocblas_create_handle(handle))
  call hipCheck(hipMemcpy(dA, c_loc(hA(1,1)), size_A * 8, hipMemcpyHostToDevice))
  call hipCheck(rocsolver_dgetrf_64(handle, M, N, dA, lda, dIpiv, dInfo))
  call hipCheck(hipMemcpy(c_loc(hLU(1,1)), dA, size_A * 8, hipMemcpyDeviceToHost))
  ! Reconstruct A = L*U (unit-lower L below the diagonal, U on/above).
  do j = 1, 3
     do i = 1, 3
        recon = 0.0d0
        do k = 1, min(i,j)
           if (k == i) then; lij = 1.0d0; else; lij = hLU(i,k); end if
           recon = recon + lij * hLU(k,j)
        end do
        error = abs(recon - hA(i,j))
        if (error > error_max) then
           write(*,*) "FAILED! (L*U)(", i, ",", j, ") = ", recon, " expected ", hA(i,j)
           call exit(1)
        end if
     end do
  end do
  call hipCheck(hipFree(dA)); call hipCheck(hipFree(dIpiv)); call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle)); call hipCheck(hipDeviceReset())
  write(*,*) "PASSED!"
end program dgetrf_64

Cholesky factorization and solve#

potrf computes the Cholesky factorization of a symmetric (or Hermitian) positive-definite matrix, writing the factor into the triangle chosen by the fill mode. The example uses rocblas_fill_upper and checks the factor against the known Cholesky root.

!!!!!!!!!!!!!/
! dpotrf example (double-precision Cholesky factorization)
! see: https:!www.netlib.org/lapack/explore-html/d1/d7a/group__double_p_ocomputational_ga2f55f604a6003d03b5cd4a0adcfb9e07.html
!
! NOTE: rocSOLVER writes the `info` output to DEVICE memory. The hipfort
! binding types the info argument as a device pointer (type(c_ptr)), so it must be backed
! by a device allocation (dInfo below) and passed as c_loc(dInfo); passing a
! host scalar faults on the GPU.
!
! With rocblas_fill_upper, the Cholesky factor is written to the upper
! triangle; the lower triangle keeps the original input values.
!!!!!!!!!!!!!!/
!
program dpotrf
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocblas_enums
  use hipfort_rocsolver

  implicit none
  integer :: i, j ! indices for iterating over results

  ! Symmetric positive-definite input (column-major)
  real(c_double) :: hA(3,3) = reshape((/4, 12, -16, 12, 37, -43, -16, -43, 98/), (/3, 3/))
  ! Reference from rocSOLVER (upper = Cholesky factor, lower = original), column-major
  real(c_double) :: hResult(3,3) = reshape((/&
     2.0000000000000000d0, 12.0000000000000000d0, -16.0000000000000000d0, &
     6.0000000000000000d0,  1.0000000000000000d0, -43.0000000000000000d0, &
    -8.0000000000000000d0,  5.0000000000000000d0,   3.0000000000000000d0/), (/3, 3/))
  integer(c_int), parameter :: N = 3
  integer(c_int), parameter :: lda = 3

  real(c_double), pointer :: dA(:,:)    ! GPU buffer for A
  integer(c_int), pointer :: dInfo(:)   ! GPU buffer for info (rocSOLVER writes to device)

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: error
  real(c_double), parameter :: error_max = 10 * epsilon(error_max)
  !
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dpotrf' (Fortran 2008 interfaces) - "

  ! Allocate device-side memory & copy memory from host to device
  call hipCheck(hipMalloc(dA,    source=hA))
  call hipCheck(hipMalloc(dInfo, 1))

  ! Create rocBLAS handle
  call hipCheck(rocblas_create_handle(handle))

  ! Compute the Cholesky factorization on the device (upper triangle)
  call hipCheck(rocsolver_dpotrf(handle, rocblas_fill_upper, N, dA, lda, c_loc(dInfo)))

  ! Copy result from device to host
  call hipCheck(hipMemcpy(hA, dA, hipMemcpyDeviceToHost))

  ! Check factor values
  do j = 1,size(hA,2)
    do i = 1,size(hA,1)
        error = abs(hA(i,j) - hResult(i,j)) / max(abs(hResult(i,j)), 1.0_c_double)
        if(error .gt. error_max) then
            write(*,*) "FAILED! Error bigger than max! Error = ", error, " hA(", i, ",", j, ") = ", hA(i,j)
            call exit
        end if
    end do
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dpotrf

potrs solves A*X = B from a potrf factorization, and posv combines the factorization and the solve in a single call. See test/f2008/rocsolver/rocsolver_dpotrs.f08 and rocsolver_dposv.f08.

QR factorization#

geqrf computes A = Q*R, storing R in the upper triangle of A and the Householder vectors that represent Q below it, with their scalar factors in ipiv/tau.

!!!!!!!!!!!!!/
! dgeqrf example
! see: http:!www.netlib.org/lapack/explore-html/df/dc5/group__variants_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html
!!!!!!!!!!!!!!/
!
program dgeqrf
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocsolver

  implicit none
  integer :: i, j ! indices for iterating over results

  ! Define our input data
  real(c_double) :: hA(3,3) = reshape((/12, 6, -4, -51, 167, 24, 4, -68, -41/), (/3, 3/))
  real(c_double) :: hResult(3,3) = reshape((/&
  -14.000000000000000,       -21.000000000000000,        14.000000000000002,&
  0.23076923076923078,       -175.00000000000000,        70.000000000000000,&
 -0.15384615384615385,        5.5555555555555559E-002,  -35.000000000000000/), shape(hResult), order=(/2,1/))
  integer(c_int), parameter :: M = 3
  integer(c_int), parameter :: N = 3
  integer(c_int), parameter :: lda = 3

  real(c_double) :: hIpiv(3) ! CPU buffer for Householder scalars

  real(c_double), pointer :: dA(:,:)   ! GPU buffer for A
  real(c_double), pointer :: dIpiv(:)  ! GPU buffer for Householder scalars

  type(c_ptr) :: handle ! rocblas_handle
    
  real :: error
  real, parameter :: error_max = 10 * epsilon(error_max)
  !
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dgeqrf' (Fortran 2008 interfaces) - "

  ! Allocate device-side memory & copy memory from host to device
  call hipCheck(hipMalloc(dA,    source=hA))
  call hipCheck(hipMalloc(dIpiv, source=hIpiv))

  ! Create rocBLAS handle
  call hipCheck(rocblas_create_handle(handle))

  ! Compute the QR factorization on the devi ce
  call hipCheck(rocsolver_dgeqrf(handle, M, N, dA, lda, dIpiv))

  ! Copy result from device to host
  call hipCheck(hipMemcpy(hA,    dA,    hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hIpiv, dIpiv, hipMemcpyDeviceToHost))

  ! Output results
  do j = 1,size(hA,2)
    do i = 1,size(hA,1)
        error = abs(hA(i,j) - hResult(i,j))
        if(error .gt. error_max) then
            write(*,*) "FAILED! Error bigger than max! Error = ", error, " hA(", i, ",", j, ") = ", hA(i,j)
            call exit
        end if
      ! print *, (hA(i,j), j=1,size(hA,2))
      ! print *, (hResult(i,j), j=1,size(hA,2))
    end do
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dIpiv))
  call hipCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())
    
  write(*,*) "PASSED!"

end program dgeqrf

Q is never formed explicitly by geqrf. Two follow-on routines use its compact representation: orgqr (ungqr for complex) generates the explicit orthogonal matrix Q, and ormqr (unmqr for complex) multiplies a given matrix by Q or Q**T without forming it. See test/f2008/rocsolver/rocsolver_dorgqr.f08 and rocsolver_dormqr.f08.

Linear least squares#

gels solves the least-squares problem min || A*X - B || (or the minimum-norm problem for underdetermined systems) using a QR or LQ factorization. The example solves a square nonsingular system and confirms the overwritten B recovers the known solution.

!!!!!!!!!!!!!/
! dgels example (double-precision least-squares solve)
! see: https:!rocm.docs.amd.com/projects/rocSOLVER/en/latest/reference/lapack.html
!
! Self-verifying: pick a known solution x, form b = A*x for a square
! nonsingular A, solve the least-squares problem with gels, and confirm the
! overwritten B recovers x.
!
! NOTE: gels writes `info` to DEVICE memory, so info is backed by a device
! allocation and passed as c_loc(dInfo); passing a host scalar faults on the GPU.
!!!!!!!!!!!!!!/
!
program dgels
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocblas_enums
  use hipfort_rocsolver

  implicit none
  integer :: i ! index for iterating over results

  integer(c_int), parameter :: M = 3
  integer(c_int), parameter :: N = 3
  integer(c_int), parameter :: nrhs = 1
  integer(c_int), parameter :: lda = 3
  integer(c_int), parameter :: ldb = 3

  ! Nonsingular input matrix (column-major)
  real(c_double) :: hA(3,3) = reshape((/1, 2, 1, 2, 5, 0, 3, 3, 8/), (/3, 3/))
  real(c_double) :: hX(3,1) = reshape((/1, 2, 3/), (/3, 1/))  ! known solution
  real(c_double) :: hB(3,1)                                    ! RHS = A*x

  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dB(:,:)
  integer(c_int), pointer :: dInfo(:)

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: error
  real(c_double), parameter :: error_max = 100 * epsilon(error_max)
  !
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dgels' (Fortran 2008 interfaces) - "

  ! Build a consistent RHS so that A*x = b
  hB = matmul(hA, hX)

  ! Allocate device-side memory & copy memory from host to device
  call hipCheck(hipMalloc(dA,    source=hA))
  call hipCheck(hipMalloc(dB,    source=hB))
  call hipCheck(hipMalloc(dInfo, 1))

  ! Create rocBLAS handle
  call hipCheck(rocblas_create_handle(handle))

  ! Solve min || A*X - B || in place (B is overwritten with the solution)
  call hipCheck(rocsolver_dgels(handle, rocblas_operation_none, M, N, nrhs, dA, lda, dB, ldb, c_loc(dInfo)))

  ! Copy result from device to host
  call hipCheck(hipMemcpy(hB, dB, hipMemcpyDeviceToHost))

  ! Verify the recovered solution matches x
  do i = 1,N
    error = abs(hB(i,1) - hX(i,1)) / max(abs(hX(i,1)), 1.0_c_double)
    if(error .gt. error_max) then
        write(*,*) "FAILED! Error bigger than max! Error = ", error, " X(", i, ") = ", hB(i,1)
        call exit
    end if
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dgels

Symmetric eigenvalues#

syev (heev for Hermitian matrices) computes the eigenvalues, and optionally the eigenvectors, of a symmetric matrix. The rocblas_evect_none/rocblas_evect_original argument selects whether eigenvectors are produced; the example requests eigenvalues only.

!!!!!!!!!!!!!!
! dsyev example (rocSOLVER)
! see: https:!rocm.docs.amd.com/projects/rocSOLVER/en/latest/reference/lapack.html
!
! Computes the eigenvalues of a real symmetric matrix. rocSOLVER writes `info`
! to DEVICE memory, so it is backed by a device allocation and passed as
! c_loc(dInfo).
!!!!!!!!!!!!!!
!
program dsyev
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocblas_enums
  use hipfort_rocsolver
  use hipfort_rocsolver_enums

  implicit none

  integer(c_int), parameter :: n = 4, lda = 4

  ! Symmetric 4x4 matrix (column-major); its trace is 10+11+12+13 = 46.
  real(c_double) :: hA(n,n) = reshape([ &
      10.0d0,  2.0d0,  3.0d0,  6.0d0, &
       2.0d0, 11.0d0,  1.0d0,  0.0d0, &
       3.0d0,  1.0d0, 12.0d0,  2.0d0, &
       6.0d0,  0.0d0,  2.0d0, 13.0d0], [n,n])
  real(c_double) :: hD(n) = 0.0d0   ! eigenvalues
  real(c_double) :: hE(n) = 0.0d0   ! workspace

  real(c_double), pointer :: dA(:,:)   ! GPU buffer for A
  real(c_double), pointer :: dD(:)     ! GPU buffer for the eigenvalues
  real(c_double), pointer :: dE(:)     ! GPU workspace
  integer(c_int), pointer :: dInfo(:)  ! GPU buffer for info

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: trace_A, error
  real(c_double), parameter :: rtol = 1.0d-9

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

  trace_A = hA(1,1) + hA(2,2) + hA(3,3) + hA(4,4)   ! = 46

  call hipCheck(hipMalloc(dA,    source=hA))
  call hipCheck(hipMalloc(dD,    source=hD))
  call hipCheck(hipMalloc(dE,    source=hE))
  call hipCheck(hipMalloc(dInfo, 1))

  call rocblasCheck(rocblas_create_handle(handle))

  ! Eigenvalues only (rocblas_evect_none); A/D/E as native arrays, info device.
  call rocsolverCheck(rocsolver_dsyev(handle, rocblas_evect_none, rocblas_fill_lower, n, dA, lda, &
                                      dD, dE, c_loc(dInfo)))

  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hD, dD, hipMemcpyDeviceToHost))

  ! An orthogonal diagonalization preserves the trace: sum of eigenvalues =
  ! trace(A). This is convention-independent (order of eigenvalues irrelevant).
  error = abs(sum(hD) - trace_A) / abs(trace_A)
  if (error > rtol) then
     write(*,*) "FAILED! sum(eigenvalues) = ", sum(hD), " expected trace(A) = ", trace_A
     call exit(1)
  end if

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dD))
  call hipCheck(hipFree(dE))
  call hipCheck(hipFree(dInfo))
  call rocblasCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dsyev

syevd/heevd solve the same problem with a divide-and-conquer algorithm, and syevj/heevj with a Jacobi algorithm. See test/f2008/rocsolver/rocsolver_dsyevd.f08 and rocsolver_dsyevj.f08.

Singular value decomposition#

gesvd computes the singular value decomposition A = U*S*V**T. The rocblas_svect_* arguments choose which of the singular-vector matrices are computed. The example requests all vectors and reconstructs A from the factors (rocSOLVER returns V**T).

!!!!!!!!!!!!!/
! dgesvd example (double-precision singular value decomposition)
! see: https:!www.netlib.org/lapack/explore-html/d1/d7e/group__double_g_esing_ga84fdf22a62b12ff364621e4713ce02f2.html
!
! Self-verifying: compute A = U * S * V (rocSOLVER returns V as V**T), then
! confirm the reconstruction matches the original A.
!
! NOTE: gesvd writes `info` to DEVICE memory, so info is backed by a device
! allocation and passed as c_loc(dInfo); passing a host scalar faults on the GPU.
! For real matrices S and E are real; A, U, V are real.
!!!!!!!!!!!!!!/
!
program dgesvd
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocblas_enums
  use hipfort_rocsolver
  use hipfort_rocsolver_enums

  implicit none
  integer :: i, j ! indices for iterating over results

  integer(c_int), parameter :: M = 3
  integer(c_int), parameter :: N = 3
  integer(c_int), parameter :: lda = 3
  integer(c_int), parameter :: ldu = 3
  integer(c_int), parameter :: ldv = 3

  ! Input matrix (column-major) and a copy kept for verification
  real(c_double) :: hA(3,3) = reshape((/4, 6, -4, -51, 167, 24, 4, -68, -41/), (/3, 3/))
  real(c_double) :: hA0(3,3)
  real(c_double) :: hS(3)      ! singular values
  real(c_double) :: hU(3,3)    ! left singular vectors
  real(c_double) :: hV(3,3)    ! right singular vectors (stored as V**T)
  real(c_double) :: hE(2)      ! superdiagonal workspace (min(M,N)-1)
  real(c_double) :: recon(3,3), Sd(3,3)

  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dS(:)
  real(c_double), pointer :: dU(:,:)
  real(c_double), pointer :: dV(:,:)
  real(c_double), pointer :: dE(:)
  integer(c_int), pointer :: dInfo(:)

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: error
  real(c_double), parameter :: error_max = 1.0d-9
  !
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dgesvd' (Fortran 2008 interfaces) - "

  hA0 = hA ! keep the original for the reconstruction check

  ! Allocate device-side memory & copy memory from host to device
  call hipCheck(hipMalloc(dA, source=hA))
  call hipCheck(hipMalloc(dS, mold=hS))
  call hipCheck(hipMalloc(dU, mold=hU))
  call hipCheck(hipMalloc(dV, mold=hV))
  call hipCheck(hipMalloc(dE, mold=hE))
  call hipCheck(hipMalloc(dInfo, 1))

  ! Create rocBLAS handle
  call hipCheck(rocblas_create_handle(handle))

  ! Compute the full SVD on the device
  call hipCheck(rocsolver_dgesvd(handle, rocblas_svect_all, rocblas_svect_all, M, N, dA, lda, &
                                 dS, dU, ldu, dV, ldv, dE, rocblas_outofplace, c_loc(dInfo)))

  ! Copy factors back to host
  call hipCheck(hipMemcpy(hS, dS, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hU, dU, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hV, dV, hipMemcpyDeviceToHost))

  ! Reconstruct A = U * diag(S) * V (V is already V**T from rocSOLVER)
  Sd = 0.0_c_double
  do i = 1,N
    Sd(i,i) = hS(i)
  end do
  recon = matmul(hU, matmul(Sd, hV))

  ! Verify reconstruction matches the original
  do j = 1,N
    do i = 1,M
        error = abs(recon(i,j) - hA0(i,j)) / max(abs(hA0(i,j)), 1.0_c_double)
        if(error .gt. error_max) then
            write(*,*) "FAILED! Error bigger than max! Error = ", error, " at (", i, ",", j, ")"
            call exit
        end if
    end do
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dS))
  call hipCheck(hipFree(dU))
  call hipCheck(hipFree(dV))
  call hipCheck(hipFree(dE))
  call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dgesvd

gesvdj computes the same decomposition with a Jacobi algorithm, which is often faster for small matrices. See test/f2008/rocsolver/rocsolver_dgesvdj.f08.

Symmetric indefinite factorization#

sytrf computes the Bunch-Kaufman factorization of a symmetric indefinite matrix, and sytrs uses that factorization to solve A*X = B. The example factorizes with rocblas_fill_upper and checks the recovered solution.

!!!!!!!!!!!!!!
! dsytrs example (solve A*X = B for a symmetric matrix, double)
! see: https:!rocm.docs.amd.com/projects/rocSOLVER/en/latest/
!
! Self-verifying: factorize a symmetric A with sytrf (Bunch-Kaufman), then solve
! A*x = b with sytrs for a right-hand side built from a known solution, and check
! that the recovered x matches. rocSOLVER writes info to DEVICE memory.
!
! sytrf has array overloads (exercised here via array pointers); sytrs is
! c_ptr-only, so the same device arrays are passed to it with c_loc.
!!!!!!!!!!!!!!
!
program dsytrs
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocblas_enums
  use hipfort_rocsolver
  implicit none
  integer :: i
  integer(c_int), parameter :: N = 4, nrhs = 1, lda = 4, ldb = 4
  real(c_double) :: hA(N,N) = reshape((/ &
      10.0d0,  2.0d0,  3.0d0,  6.0d0, &
       2.0d0, 11.0d0,  1.0d0,  0.0d0, &
       3.0d0,  1.0d0, 12.0d0,  2.0d0, &
       6.0d0,  0.0d0,  2.0d0, 13.0d0/), (/N,N/))
  real(c_double) :: hX(N) = (/1.0d0, 2.0d0, 3.0d0, 4.0d0/)   ! known solution
  real(c_double) :: hB(N)                                     ! RHS = A*x
  integer(c_int), target :: hInfo(1)
  integer(c_int) :: hIpiv(N)
  real(c_double), pointer :: dA(:,:), dB(:)
  integer(c_int), pointer :: dIpiv(:)
  type(c_ptr) :: dInfo
  type(c_ptr) :: handle
  real(c_double), parameter :: error_max = 1.0d-9
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dsytrs' (Fortran 2008 interfaces) - "

  hB = matmul(hA, hX)

  call hipCheck(hipMalloc(dA, source=hA))
  call hipCheck(hipMalloc(dB, source=hB))
  call hipCheck(hipMalloc(dIpiv, mold=hIpiv))
  call hipCheck(hipMalloc(dInfo, 4_c_size_t))

  call hipCheck(rocblas_create_handle(handle))
  ! sytrf: array overload (dA, dIpiv passed directly)
  call hipCheck(rocsolver_dsytrf(handle, rocblas_fill_upper, N, dA, lda, dIpiv, dInfo))
  ! sytrs: c_ptr-only, pass the same device arrays via c_loc
  call hipCheck(rocsolver_dsytrs(handle, rocblas_fill_upper, N, nrhs, &
       c_loc(dA), lda, c_loc(dIpiv), c_loc(dB), ldb))
  call hipCheck(hipMemcpy(c_loc(hInfo(1)), dInfo, 4_c_size_t, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hB, dB, hipMemcpyDeviceToHost))

  if (hInfo(1) /= 0) then
     write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if
  do i = 1, N
     if (abs(hB(i) - hX(i)) > error_max) then
        write(*,*) "FAILED! x(", i, ") = ", hB(i), " expected ", hX(i); call exit(1)
     end if
  end do

  call hipCheck(hipFree(dA)); call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dIpiv)); call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle)); call hipCheck(hipDeviceReset())
  write(*,*) "PASSED!"
end program dsytrs

Triangular inverse#

trtri inverts a triangular matrix in place. The example inverts an upper-triangular matrix and checks U * U^-1 == I.

!!!!!!!!!!!!!!
! dtrtri example (rocSOLVER)
! see: https:!rocm.docs.amd.com/projects/rocSOLVER/en/latest/reference/lapack.html
!
! Inverts an upper-triangular matrix in place and checks U * U^-1 == I.
! rocSOLVER writes `info` to DEVICE memory, so it is backed by a device
! allocation and passed as c_loc(dInfo).
!!!!!!!!!!!!!!
!
program dtrtri
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocblas_enums
  use hipfort_rocsolver

  implicit none

  integer(c_int), parameter :: n = 3, lda = 3

  ! Upper-triangular matrix (lower part zero), column-major:
  !   [2 1 1; 0 2 1; 0 0 2]
  real(c_double) :: hU(n,n) = reshape([ &
      2.0d0, 0.0d0, 0.0d0, &
      1.0d0, 2.0d0, 0.0d0, &
      1.0d0, 1.0d0, 2.0d0], [n,n])
  real(c_double) :: hUinv(n,n) = 0.0d0

  real(c_double), pointer :: dA(:,:)   ! GPU buffer for A (holds U^-1 on output)
  integer(c_int), pointer :: dInfo(:)  ! GPU buffer for info

  type(c_ptr) :: handle ! rocblas_handle

  integer :: i, j, l
  real(c_double) :: prod, error
  real(c_double), parameter :: error_max = 1.0d-10

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

  call hipCheck(hipMalloc(dA,    source=hU))
  call hipCheck(hipMalloc(dInfo, 1))

  call rocblasCheck(rocblas_create_handle(handle))

  call rocsolverCheck(rocsolver_dtrtri(handle, rocblas_fill_upper, rocblas_diagonal_non_unit, &
                                       n, dA, lda, c_loc(dInfo)))

  call hipCheck(hipDeviceSynchronize())
  call hipCheck(hipMemcpy(hUinv, dA, hipMemcpyDeviceToHost))

  ! U and U^-1 are both upper triangular (zero below), so U * U^-1 == I exactly.
  do j = 1, n
     do i = 1, n
        prod = 0.0d0
        do l = 1, n
           prod = prod + hU(i,l) * hUinv(l,j)
        end do
        error = abs(prod - merge(1.0d0, 0.0d0, i == j))
        if (error > error_max) then
           write(*,*) "FAILED! (U*Uinv)(", i, ",", j, ") = ", prod
           call exit(1)
        end if
     end do
  end do

  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dInfo))
  call rocblasCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dtrtri

Reductions to condensed form#

Several rocSOLVER routines reduce a matrix to a condensed form used inside the eigenvalue and SVD algorithms. getrf aside, these are lower-level building blocks:

  • gebrd reduces a general matrix to bidiagonal form.

  • sytrd reduces a symmetric matrix to tridiagonal form, and latrd reduces a leading block of it.

  • sterf computes the eigenvalues of a symmetric tridiagonal matrix, and steqr/stedc compute its eigenvalues and eigenvectors.

  • larft forms the triangular factor of a block of Householder reflectors.

!!!!!!!!!!!!!/
! dgebrd example (double-precision bidiagonal reduction)
! see: https:!rocm.docs.amd.com/projects/rocSOLVER/en/latest/reference/lapack.html
!
! gebrd reduces A to bidiagonal form B = Q**T * A * P. We check the diagonal D
! and superdiagonal E of B against reference values, and confirm the reduction
! preserves the Frobenius norm (sum(D**2)+sum(E**2) = ||A||_F**2), which holds
! because Q and P are orthogonal.
!!!!!!!!!!!!!!/
!
program dgebrd
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocsolver

  implicit none
  integer :: i ! index for iterating over results

  integer(c_int), parameter :: M = 3
  integer(c_int), parameter :: N = 3
  integer(c_int), parameter :: lda = 3

  ! Input matrix (column-major)
  real(c_double) :: hA(3,3) = reshape((/1, 4, 7, 2, 5, 8, 3, 6, 10/), (/3, 3/))
  real(c_double) :: hD(3)      ! diagonal of B
  real(c_double) :: hE(2)      ! superdiagonal of B
  real(c_double) :: hTauq(3)   ! Householder scalars (Q)
  real(c_double) :: hTaup(3)   ! Householder scalars (P)

  ! Reference bidiagonal entries (from rocSOLVER)
  real(c_double), parameter :: refD(3) = &
    (/-8.1240384046359608_c_double, -1.7704849623785914_c_double, 0.2085724989394371_c_double/)
  real(c_double), parameter :: refE(2) = &
    (/15.3213062185449527_c_double, 0.2818798826684673_c_double/)
  real(c_double), parameter :: normA2 = 304.0_c_double

  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dD(:)
  real(c_double), pointer :: dE(:)
  real(c_double), pointer :: dTauq(:)
  real(c_double), pointer :: dTaup(:)

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: error
  real(c_double), parameter :: error_max = 1.0d-9
  !
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dgebrd' (Fortran 2008 interfaces) - "

  ! Allocate device-side memory & copy memory from host to device
  call hipCheck(hipMalloc(dA,    source=hA))
  call hipCheck(hipMalloc(dD,    mold=hD))
  call hipCheck(hipMalloc(dE,    mold=hE))
  call hipCheck(hipMalloc(dTauq, mold=hTauq))
  call hipCheck(hipMalloc(dTaup, mold=hTaup))

  ! Create rocBLAS handle
  call hipCheck(rocblas_create_handle(handle))

  ! Reduce A to bidiagonal form
  call hipCheck(rocsolver_dgebrd(handle, M, N, dA, lda, dD, dE, dTauq, dTaup))

  ! Copy results back to host
  call hipCheck(hipMemcpy(hD, dD, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hE, dE, hipMemcpyDeviceToHost))

  ! Verify the bidiagonal entries against the reference
  do i = 1,3
    error = abs(hD(i) - refD(i))
    if(error .gt. error_max) then
        write(*,*) "FAILED! Error bigger than max! Error = ", error, " D(", i, ") = ", hD(i)
        call exit
    end if
  end do
  do i = 1,2
    error = abs(hE(i) - refE(i))
    if(error .gt. error_max) then
        write(*,*) "FAILED! Error bigger than max! Error = ", error, " E(", i, ") = ", hE(i)
        call exit
    end if
  end do

  ! Cross-check the norm-preserving invariant
  error = abs(sum(hD*hD) + sum(hE*hE) - normA2)
  if(error .gt. 1.0d-6) then
      write(*,*) "FAILED! Norm not preserved! Error = ", error
      call exit
  end if

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dD))
  call hipCheck(hipFree(dE))
  call hipCheck(hipFree(dTauq))
  call hipCheck(hipFree(dTaup))
  call hipCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dgebrd

The tridiagonal eigenvalue solver sterf takes only the diagonal and off-diagonal of the tridiagonal matrix:

!!!!!!!!!!!!!/
! dsterf example
!!!!!!!!!!!!!!/
!
program dsterf
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_rocblas
  use hipfort_rocsolver

  implicit none
  integer :: i ! index for iterating over results

  ! Define our input data
  real(c_double) :: hD(3) = (/ 2, 2, 2 /)
  real(c_double) :: hE(2) = (/ -1, -1 /)
  real(c_double) :: hResult(3) = (/ 0.58578643762690485d0, 2.0d0, 3.41421356237309515d0 /)
  integer(c_int) :: hInfo = -1

  integer(c_int), parameter :: n = 3

  real(c_double), pointer :: dD(:)   ! GPU buffer for D
  real(c_double), pointer :: dE(:)   ! GPU buffer for E
  integer(c_int), pointer :: dInfo   ! GPU buffer for myInfo

  type(c_ptr) :: handle ! rocblas_handle

  real :: error
  real, parameter :: error_max = 10 * epsilon(error_max)
  !
  write(*,"(a)",advance="no") "-- Running test 'rocsolver_dsterf' (Fortran 2008 interfaces) - "

  ! Allocate device-side memory & copy memory from host to device
  call hipCheck(hipMalloc(dD, source=hD))
  call hipCheck(hipMalloc(dE, source=hE))
  call hipCheck(hipMalloc(dInfo))

  ! Create rocBLAS handle
  call hipCheck(rocblas_create_handle(handle))

  ! Compute eigenvalues.
  ! `dInfo` is passed as a device pointer (c_loc), which is what the
  ! myInfo -> c_ptr binding fix enables.
  call hipCheck(rocsolver_dsterf(handle, n, dD, dE, c_loc(dInfo)))

  ! Copy result from device to host
  call hipCheck(hipMemcpy(hD,    dD,    hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))

  ! Check and output results
  if(hInfo .gt. 0) then
    write(*,*) "FAILED! ", n, " elements of E did not converge to 0."
    call exit(1)
  else
    do i = 1,n
      error = abs(hD(i) - hResult(i))
        if(error .gt. error_max) then
            write(*,*) "FAILED! Error bigger than max! Error = ", error, " hD(", i, ") = ", hD(i)
            call exit(1)
        end if
    end do
  end if

  ! Clean up
  call hipCheck(hipFree(dD))
  call hipCheck(hipFree(dE))
  call hipCheck(hipFree(dInfo))
  call hipCheck(rocblas_destroy_handle(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"
  write(*,*) "Eigenvalues in increasing order:"
  do i = 1,n
    write(*,*) hD(i)
  end do

end program dsterf

See test/f2008/rocsolver for the sytrd, latrd, steqr, stedc, and larft examples.

Contents