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 complete and self-contained, and is built and run as part of the hipFORT test suite. The Fortran 2008 version of each program lives in test/f2008/rocsolver. Most programs also have an equivalent Fortran 2003 version, which uses type(c_ptr) device pointers and explicit byte counts instead of Fortran array pointers, in test/f2003/rocsolver.

hipSOLVER offers the same functionality through an API that follows cuSOLVER; see the hipSOLVER examples.

Each routine is provided in the four LAPACK precisions where it 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 program:

• 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, which holds the pivot indices or, for the QR routines, the Householder scalars that LAPACK calls tau, are outputs written on the device and are passed as device buffers.

• Enumerators select variants. rocblas_operation_none / rocblas_operation_transpose choose whether a routine works on A or A**T, and rocblas_fill_upper / rocblas_fill_lower choose the stored triangle; both come from hipfort_rocblas_enums. The rocblas_evect_* and rocblas_svect_* selectors, which choose whether vectors are computed, are rocSOLVER enumerators from hipfort_rocsolver_enums despite their prefix.

• Every call returns a status code. The hipfort_check module provides hipCheck for HIP calls, rocblasCheck for the rocBLAS handle calls, and rocsolverCheck for rocSOLVER calls; all abort on failure.

Building and running

The programs 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 program factorizes a matrix and checks the packed factors and the pivot indices against reference values.

!!!!!!!!!!!!!/
! 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))
  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)
  integer(c_int) :: hInfo(1)

  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 rocblasCheck(rocblas_create_handle(handle))

  ! Compute the LU factorization on the device
  call rocsolverCheck(rocsolver_dgetrf(handle, M, N, dA, lda, dIpiv, c_loc(dInfo)))
  call hipCheck(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))
  if (hInfo(1) /= 0) then
    write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if

  ! 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(1)
        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(1)
    end if
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dIpiv))
  call hipCheck(hipFree(dInfo))
  call rocblasCheck(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 program 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(:)
  integer(c_int) :: hInfo(1)

  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 rocblasCheck(rocblas_create_handle(handle))

  ! Factorize, then solve A*X = B in place (B is overwritten with the solution)
  call rocsolverCheck(rocsolver_dgetrf(handle, N, N, dA, lda, dIpiv, c_loc(dInfo)))
  call hipCheck(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))
  if (hInfo(1) /= 0) then
    write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if
  call rocsolverCheck(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(1)
    end if
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dIpiv))
  call hipCheck(hipFree(dInfo))
  call rocblasCheck(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, ipiv as one device buffer with a per-matrix stride (strideP), and info as a device array 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 rocblasCheck(rocblas_create_handle(handle))
  call rocsolverCheck(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 rocblasCheck(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 _64 routines have the same Fortran array overloads as the 32-bit ones; this example passes all its device buffers, including the int64 ipiv and info, as type(c_ptr), so it calls the C binding.

!!!!!!!!!!!!!!
! 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(c_int64_t), target :: hInfo(1)
  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 rocblasCheck(rocblas_create_handle(handle))
  call hipCheck(hipMemcpy(dA, c_loc(hA(1,1)), size_A * 8, hipMemcpyHostToDevice))
  call rocsolverCheck(rocsolver_dgetrf_64(handle, M, N, dA, lda, dIpiv, dInfo))
  call hipCheck(hipMemcpy(c_loc(hInfo(1)), dInfo, 8_c_size_t, hipMemcpyDeviceToHost))
  if (hInfo(1) /= 0) then
    write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if
  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 rocblasCheck(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 program 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)
  integer(c_int) :: hInfo(1)

  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 rocblasCheck(rocblas_create_handle(handle))

  ! Compute the Cholesky factorization on the device (upper triangle)
  call rocsolverCheck(rocsolver_dpotrf(handle, rocblas_fill_upper, N, dA, lda, c_loc(dInfo)))
  call hipCheck(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))
  if (hInfo(1) /= 0) then
    write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if

  ! 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(1)
        end if
    end do
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dInfo))
  call rocblasCheck(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 the argument rocSOLVER names ipiv (tau in LAPACK).

!!!!!!!!!!!!!/
! 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.000000000000000d0,       -21.000000000000000d0,        14.000000000000002d0,&
  0.23076923076923078d0,       -175.00000000000000d0,        70.000000000000000d0,&
 -0.15384615384615385d0,        5.5555555555555559D-002,  -35.000000000000000d0/), 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(c_double) :: error
  real(c_double), parameter :: error_max = 100 * 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 rocblasCheck(rocblas_create_handle(handle))

  ! Compute the QR factorization on the device
  call rocsolverCheck(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)) / 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(1)
        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 rocblasCheck(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 program 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(:)
  integer(c_int) :: hInfo(1)

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: error
  real(c_double), parameter :: error_max = 10000 * epsilon(error_max)   ! kappa_2(A) ~ 4.5e2; QR forward error ~ c*kappa_2*eps
  !
  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 rocblasCheck(rocblas_create_handle(handle))

  ! Solve min || A*X - B || in place (B is overwritten with the solution)
  call rocsolverCheck(rocsolver_dgels(handle, rocblas_operation_none, M, N, nrhs, dA, lda, dB, ldb, c_loc(dInfo)))
  call hipCheck(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))
  if (hInfo(1) /= 0) then
    write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if

  ! 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(1)
    end if
  end do

  ! Clean up
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dInfo))
  call rocblasCheck(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 program 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
  integer(c_int) :: hInfo(1)

  type(c_ptr) :: handle ! rocblas_handle

  real(c_double) :: trace_A, frob2_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
  frob2_A = sum(hA**2)   ! ||A||_F**2, from the original A

  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(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))
  if (hInfo(1) /= 0) then
    write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if

  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

  ! The squared Frobenius norm is also invariant: sum of squared eigenvalues = ||A||_F**2.
  error = abs(sum(hD**2) - frob2_A) / abs(frob2_A)
  if (error > rtol) then
     write(*,*) "FAILED! sum(eigenvalues**2) = ", sum(hD**2), " expected ||A||_F**2 = ", frob2_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 program 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(:)
  integer(c_int) :: hInfo(1)

  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 rocblasCheck(rocblas_create_handle(handle))

  ! Compute the full SVD on the device
  call rocsolverCheck(rocsolver_dgesvd(handle, rocblas_svect_all, rocblas_svect_all, M, N, dA, lda, &
                                 dS, dU, ldu, dV, ldv, dE, rocblas_outofplace, c_loc(dInfo)))
  call hipCheck(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))
  if (hInfo(1) /= 0) then
    write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if

  ! 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(1)
        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 rocblasCheck(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 program 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 rocblasCheck(rocblas_create_handle(handle))
  ! sytrf: array overload (dA, dIpiv passed directly)
  call rocsolverCheck(rocsolver_dsytrf(handle, rocblas_fill_upper, N, dA, lda, dIpiv, dInfo))
  ! sytrs: c_ptr-only, pass the same device arrays via c_loc
  call rocsolverCheck(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 rocblasCheck(rocblas_destroy_handle(handle)); call hipCheck(hipDeviceReset())
  write(*,*) "PASSED!"
end program dsytrs

Triangular inverse

trtri inverts a triangular matrix in place. The program 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
  integer(c_int) :: hInfo(1)

  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(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))
  if (hInfo(1) /= 0) then
    write(*,*) "FAILED! info = ", hInfo(1), " (expected 0)"; call exit(1)
  end if

  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, or solve the condensed problem itself. 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 rocblasCheck(rocblas_create_handle(handle))

  ! Reduce A to bidiagonal form
  call rocsolverCheck(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(1)
    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(1)
    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(1)
  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 rocblasCheck(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(c_double) :: error
  real(c_double), parameter :: error_max = 100 * 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 rocblasCheck(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 rocsolverCheck(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! ", hInfo, " elements of E did not converge to 0."
    call exit(1)
  else
    do i = 1,n
      error = abs(hD(i) - hResult(i)) / abs(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 rocblasCheck(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 programs.