hipSOLVER examples

hipSOLVER examples#

hipSOLVER is a thin layer over rocSOLVER whose API follows cuSOLVER. hipFORT exposes it through the hipfort_hipsolver module.

Every program on this page is complete and self-contained, and is built and run as part of the hipFORT test suite. The Fortran 2008 programs live in test/f2008/hipsolver, and the Fortran 2003 versions, which use type(c_ptr) device pointers and explicit byte counts instead of Fortran array pointers, live in test/f2003/hipsolver.

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

Many routines come in the four LAPACK precisions: s (real single), d (real double), c (complex single), and z (complex double). Where several are provided, this page shows the double-precision program of the group; the other precisions differ only in the host data type and the hipsolver prefix letter.

Solver workflow#

hipSOLVER routines need an explicit GPU workspace, sized by a separate query rather than by an lwork = -1 call. A typical call follows the same sequence as cuSOLVER:

  1. Create a handle with hipsolverCreate.

  2. Query the workspace size with the routine’s _bufferSize companion (for example hipsolverDgetrf_bufferSize). The query returns a number of elements, so allocate lwork * sizeof(element) bytes on the device — lwork * 8 for a double precision routine.

  3. Run the routine, passing the workspace and its size.

  4. Read back the device info output to check for success.

  5. Release the handle with hipsolverDestroy.

Keep the following conventions in mind:

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

  • An explicit workspace. Most routines take a device work buffer and its length. Size it with the matching _bufferSize query rather than guessing; the buffer stays a bare type(c_ptr) in both dialects.

  • The info output lives in device memory. hipSOLVER writes the factorization status to a device pointer, so it must be backed by a device allocation, not a host scalar. For the batched routines info is an array of batch_count integers on the device.

  • Enumerators select variants. HIPSOLVER_FILL_MODE_UPPER / HIPSOLVER_FILL_MODE_LOWER choose the stored triangle, and HIPSOLVER_EIG_MODE_NOVECTOR / HIPSOLVER_EIG_MODE_VECTOR choose whether eigenvectors are computed. The SVD job arguments are character(c_char) job codes ('N', 'A', 'S', 'O') passed by value.

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

Building and running#

The programs only need the hipsolver and hip hipFORT components:

find_package(hipfort REQUIRED COMPONENTS hip hipsolver)

add_executable(my_solver hipsolver_dgetrf.f08)
target_link_libraries(my_solver PRIVATE hipfort::hipsolver 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 queries the workspace with hipsolverDgetrf_bufferSize, factorizes, and checks the packed factors and the pivot indices against reference values.

!!!!!!!!!!!!!!
! hipsolver dgetrf example (double-precision LU factorization, Fortran 2008 interfaces)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! f2008 style: device buffers are native Fortran array pointers allocated with
! hipMalloc(source=/mold=); A and the pivots are passed as typed arrays. The
! workspace stays a bare type(c_ptr); devInfo is a device-backed scalar pointer.
!!!!!!!!!!!!!!
!
program hipsolver_dgetrf
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver

  implicit none
  integer :: i, j

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

  ! Input matrix (column-major) and the expected packed LU (same reference as
  ! the rocSOLVER dgetrf test; hipSOLVER wraps rocSOLVER)
  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((/&
    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) :: hIpiv(3), hInfo

  real(c_double),  pointer :: dA(:,:)   ! GPU buffer for A (holds packed LU on output)
  integer(c_int),  pointer :: dIpiv(:)  ! GPU buffer for pivots
  integer(c_int),  pointer :: dInfo     ! GPU scalar for devInfo (written on device)
  type(c_ptr) :: dWork                  ! opaque workspace
  type(c_ptr) :: handle = c_null_ptr
  integer(c_int) :: lwork

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

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

  call hipsolverCheck(hipsolverCreate(handle))

  ! Allocate device memory as native arrays and copy A to device
  call hipCheck(hipMalloc(dA,    source=hA))
  call hipCheck(hipMalloc(dIpiv, mold=hIpiv))
  call hipCheck(hipMalloc(dInfo))

  ! Query and allocate the workspace
  call hipsolverCheck(hipsolverDgetrf_bufferSize(handle, M, N, dA, lda, lwork))
  call hipCheck(hipMalloc(dWork, max(int(lwork,c_size_t) * 8, 1_c_size_t)))

  ! Compute the LU factorization (A/pivots as native arrays, devInfo by reference)
  call hipsolverCheck(hipsolverDgetrf(handle, M, N, dA, lda, dWork, lwork, dIpiv, dInfo))

  ! Copy results back to host
  call hipCheck(hipMemcpy(hA,    dA,    hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hIpiv, dIpiv, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))

  ! Check info
  if(hInfo /= 0) then
    write(*,*) "FAILED! info = ", hInfo, " (expected 0)"
    call exit(1)
  end if

  ! 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 hipCheck(hipFree(dWork))
  call hipsolverCheck(hipsolverDestroy(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program hipsolver_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.

!!!!!!!!!!!!!!
! hipsolver dgetrs example (solve A*X=B after LU, Fortran 2008 interfaces)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Factorizes A with getrf, then solves A*X=B with getrs and checks X.
! Native-array f2008 form; workspaces via *_bufferSize; devInfo device-backed.
!!!!!!!!!!!!!!
!
program hipsolver_dgetrs
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver
  implicit none
  integer :: i
  integer(c_int), parameter :: N = 3, nrhs = 1, lda = 3, ldb = 3
  real(c_double) :: hA(N,N)     = reshape((/2, 1, 0,  1, 2, 1,  0, 1, 2/), (/N,N/))
  real(c_double) :: hB(N,nrhs)  = reshape((/3, 4, 3/), (/N,nrhs/))
  real(c_double) :: hX_ref(N)   = (/1, 1, 1/)
  real(c_double), pointer :: dA(:,:), dB(:,:)
  integer(c_int), pointer :: dIpiv(:)
  integer(c_int), pointer :: dInfo
  integer(c_int) :: hIpiv(N)
  type(c_ptr) :: dWork1, dWork2, handle = c_null_ptr
  integer(c_int) :: lwork1, lwork2
  real(c_double) :: error
  real(c_double), parameter :: error_max = 100 * epsilon(error_max)
  write(*,"(a)",advance="no") "-- Running test 'hipsolver_dgetrs' (Fortran 2008 interfaces) - "
  call hipsolverCheck(hipsolverCreate(handle))
  call hipCheck(hipMalloc(dA, source=hA))
  call hipCheck(hipMalloc(dB, source=hB))
  call hipCheck(hipMalloc(dIpiv, mold=hIpiv))
  call hipCheck(hipMalloc(dInfo))
  ! LU factorization
  call hipsolverCheck(hipsolverDgetrf_bufferSize(handle, N, N, dA, lda, lwork1))
  call hipCheck(hipMalloc(dWork1, max(int(lwork1,c_size_t) * 8, 1_c_size_t)))
  call hipsolverCheck(hipsolverDgetrf(handle, N, N, dA, lda, dWork1, lwork1, dIpiv, dInfo))
  ! Solve
  call hipsolverCheck(hipsolverDgetrs_bufferSize(handle, HIPSOLVER_OP_N, N, nrhs, dA, lda, dIpiv, dB, ldb, lwork2))
  call hipCheck(hipMalloc(dWork2, max(int(lwork2,c_size_t) * 8, 1_c_size_t)))
  call hipsolverCheck(hipsolverDgetrs(handle, HIPSOLVER_OP_N, N, nrhs, dA, lda, dIpiv, dB, ldb, dWork2, lwork2, dInfo))
  call hipCheck(hipMemcpy(hB, dB, hipMemcpyDeviceToHost))
  do i = 1, N
     error = abs(hB(i,1) - hX_ref(i)) / max(abs(hX_ref(i)), 1.0d0)
     if (error > error_max) then
        write(*,*) "FAILED! X(", i, ") = ", hB(i,1), " expected ", hX_ref(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(hipFree(dWork1)); call hipCheck(hipFree(dWork2))
  call hipsolverCheck(hipsolverDestroy(handle))
  write(*,*) "PASSED!"
end program hipsolver_dgetrs

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 HIPSOLVER_FILL_MODE_UPPER and checks the factor against the known Cholesky root.

!!!!!!!!!!!!!!
! hipsolver dpotrf example (Cholesky factorization, Fortran 2008 interfaces)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Native-array f2008 form: A is a Fortran array pointer; the workspace is a
! type(c_ptr) sized by hipsolverDpotrf_bufferSize; devInfo is device-backed.
!!!!!!!!!!!!!!
!
program hipsolver_dpotrf
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver
  implicit none
  integer :: i, j
  integer(c_int), parameter :: N = 3, lda = 3
  real(c_double) :: hA(3,3) = reshape((/4, 12, -16, 12, 37, -43, -16, -43, 98/), (/3, 3/))
  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/))
  real(c_double), pointer :: dA(:,:)
  integer(c_int), pointer :: dInfo
  type(c_ptr) :: dWork, handle = c_null_ptr
  integer(c_int) :: lwork
  real(c_double) :: error
  real(c_double), parameter :: error_max = 10 * epsilon(error_max)
  write(*,"(a)",advance="no") "-- Running test 'hipsolver_dpotrf' (Fortran 2008 interfaces) - "
  call hipsolverCheck(hipsolverCreate(handle))
  call hipCheck(hipMalloc(dA, source=hA))
  call hipCheck(hipMalloc(dInfo))
  call hipsolverCheck(hipsolverDpotrf_bufferSize(handle, HIPSOLVER_FILL_MODE_UPPER, N, dA, lda, lwork))
  call hipCheck(hipMalloc(dWork, max(int(lwork,c_size_t) * 8, 1_c_size_t)))
  call hipsolverCheck(hipsolverDpotrf(handle, HIPSOLVER_FILL_MODE_UPPER, N, dA, lda, dWork, lwork, dInfo))
  call hipCheck(hipMemcpy(hA, dA, hipMemcpyDeviceToHost))
  do j = 1,3
    do i = 1,3
      error = abs(hA(i,j) - hResult(i,j)) / max(abs(hResult(i,j)), 1.0_c_double)
      if (error > error_max) then
        write(*,*) "FAILED! hA(", i, ",", j, ") = ", hA(i,j)
        call exit(1)
      end if
    end do
  end do
  call hipCheck(hipFree(dA)); call hipCheck(hipFree(dInfo)); call hipCheck(hipFree(dWork))
  call hipsolverCheck(hipsolverDestroy(handle))
  write(*,*) "PASSED!"
end program hipsolver_dpotrf

potrs solves A*X = B from a potrf factorization. The program forms b = A*x for a known x and confirms the solve recovers it.

!!!!!!!!!!!!!/
! hipsolverDpotrs example (double-precision Cholesky solve)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Self-verifying: pick a known solution x, form b = A*x for a symmetric positive
! definite A, factorize with potrf, solve A*X = b with potrs, and confirm X
! recovers x.
!!!!!!!!!!!!!!/
!
program dpotrs
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver
  use hipfort_hipsolver_enums

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

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

  ! Symmetric positive definite input (column-major)
  real(c_double) :: hA(3,3) = reshape((/4, 2, 2, 2, 5, 3, 2, 3, 6/), (/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

  type(c_ptr) :: handle = c_null_ptr
  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dB(:,:)
  integer(c_int), pointer :: dInfo(:)
  type(c_ptr) :: dWork
  integer(c_int) :: lwork_f, lwork_s, lwork

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

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

  call hipsolverCheck(hipsolverCreate(handle))

  ! 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))

  ! Workspace big enough for both potrf and potrs
  call hipsolverCheck(hipsolverDpotrf_bufferSize(handle, HIPSOLVER_FILL_MODE_LOWER, N, dA, lda, lwork_f))
  call hipsolverCheck(hipsolverDpotrs_bufferSize(handle, HIPSOLVER_FILL_MODE_LOWER, N, nrhs, dA, lda, dB, ldb, lwork_s))
  lwork = max(lwork_f, lwork_s)
  call hipCheck(hipMalloc(dWork, int(lwork,c_size_t) * 8))

  ! Factorize A = L*L**T, then solve A*X = B in place
  call hipsolverCheck(hipsolverDpotrf(handle, HIPSOLVER_FILL_MODE_LOWER, N, dA, lda, dWork, lwork, dInfo(1)))
  call hipsolverCheck(hipsolverDpotrs(handle, HIPSOLVER_FILL_MODE_LOWER, N, nrhs, dA, lda, dB, ldb, dWork, lwork, dInfo(1)))

  ! 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(dWork))
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dB))
  call hipCheck(hipFree(dInfo))
  call hipsolverCheck(hipsolverDestroy(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dpotrs

Batched Cholesky#

potrfBatched factorizes many matrices with one call. The batched API takes A as an array of device pointers that itself lives in device memory: each matrix is allocated on the device, their device addresses are collected in a host array, and that array is copied to a device buffer whose address is passed as A. info is a device array indexed by batch.

!!!!!!!!!!!!!!
! hipsolver dpotrfBatched example (batched Cholesky factorization)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Factorizes a batch of SPD matrices. The batched API takes A as an array of
! device pointers that itself lives in DEVICE memory: each matrix is allocated on
! the device, their device addresses are collected in a host array, and that
! array is copied to a device buffer whose address (by value) is passed as A.
! Checks info == 0 for every batch entry and the Cholesky diagonal.
!!!!!!!!!!!!!!
!
program dpotrfbatched
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver
  implicit none
  integer :: b
  integer(c_int), parameter :: N = 3, lda = 3, batch = 2
  ! SPD A = diag(4, 9, 16) -> upper Cholesky diag = 2, 3, 4
  real(c_double), target :: hA(N,N) = reshape((/ &
      4.0d0, 0.0d0, 0.0d0, &
      0.0d0, 9.0d0, 0.0d0, &
      0.0d0, 0.0d0,16.0d0/), (/N,N/))
  real(c_double), target :: hOut(N,N)
  integer(c_int) :: hInfo(batch)
  type(c_ptr), target :: hostPtrs(batch)
  type(c_ptr) :: dA1, dA2, dPtrArray, dWork
  integer(c_int), pointer :: dInfo(:)
  type(c_ptr) :: handle
  integer(c_int) :: lwork
  integer(c_size_t) :: ptrbytes
  write(*,"(a)",advance="no") "-- Running test 'hipsolver_dpotrfbatched' (Fortran 2008 interfaces) - "

  call hipsolverCheck(hipsolverCreate(handle))
  call hipCheck(hipMalloc(dA1, int(N*N,c_size_t) * 8))
  call hipCheck(hipMalloc(dA2, int(N*N,c_size_t) * 8))
  call hipCheck(hipMemcpy(dA1, c_loc(hA(1,1)), int(N*N,c_size_t) * 8, hipMemcpyHostToDevice))
  call hipCheck(hipMemcpy(dA2, c_loc(hA(1,1)), int(N*N,c_size_t) * 8, hipMemcpyHostToDevice))

  ! Host array of device addresses, copied to a device-resident pointer array.
  hostPtrs(1) = dA1
  hostPtrs(2) = dA2
  ptrbytes = int(batch,c_size_t) * c_sizeof(c_null_ptr)
  call hipCheck(hipMalloc(dPtrArray, ptrbytes))
  call hipCheck(hipMemcpy(dPtrArray, c_loc(hostPtrs), ptrbytes, hipMemcpyHostToDevice))

  call hipCheck(hipMalloc(dInfo, batch))

  call hipsolverCheck(hipsolverDpotrfBatched_bufferSize(handle, HIPSOLVER_FILL_MODE_UPPER, &
       N, dPtrArray, lda, lwork, batch))
  call hipCheck(hipMalloc(dWork, int(max(lwork,1),c_size_t) * 8))
  call hipsolverCheck(hipsolverDpotrfBatched(handle, HIPSOLVER_FILL_MODE_UPPER, &
       N, dPtrArray, lda, dWork, lwork, dInfo(1), batch))
  call hipCheck(hipDeviceSynchronize())

  call hipCheck(hipMemcpy(hInfo, dInfo, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(c_loc(hOut(1,1)), dA1, int(N*N,c_size_t) * 8, hipMemcpyDeviceToHost))

  do b = 1, batch
     if (hInfo(b) /= 0) then
        write(*,*) "FAILED! info(", b, ") = ", hInfo(b), " (expected 0)"; call exit(1)
     end if
  end do
  if (abs(hOut(1,1)-2.0d0) > 1.0d-9 .or. abs(hOut(2,2)-3.0d0) > 1.0d-9 .or. &
      abs(hOut(3,3)-4.0d0) > 1.0d-9) then
     write(*,*) "FAILED! chol diag = ", hOut(1,1), hOut(2,2), hOut(3,3), " expected 2 3 4"; call exit(1)
  end if

  call hipCheck(hipFree(dA1)); call hipCheck(hipFree(dA2)); call hipCheck(hipFree(dPtrArray))
  call hipCheck(hipFree(dWork)); call hipCheck(hipFree(dInfo))
  call hipsolverCheck(hipsolverDestroy(handle)); call hipCheck(hipDeviceReset())
  write(*,*) "PASSED!"
end program dpotrfbatched

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 tau.

!!!!!!!!!!!!!/
! hipsolverDgeqrf example (double-precision QR factorization)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Self-verifying: geqrf overwrites the upper triangle of A with R. Because
! A = Q*R with Q orthogonal, A**T * A = R**T * R, so we recover R from the
! output, form R**T * R, and compare against A0**T * A0 computed on the host.
!
! hipSOLVER uses an explicit workspace sized by the matching _bufferSize query,
! and devInfo lives in device memory.
!!!!!!!!!!!!!!/
!
program dgeqrf
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver

  implicit none
  integer :: i, j, l ! indices 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) :: hA0(3,3)      ! original kept for verification
  real(c_double) :: hTau(3)
  real(c_double) :: R(3,3), lhs, rhs

  type(c_ptr) :: handle = c_null_ptr
  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dTau(:)
  integer(c_int), pointer :: dInfo(:)
  type(c_ptr) :: dWork
  integer(c_int) :: lwork

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

  hA0 = hA ! keep original for the A**T*A = R**T*R check

  call hipsolverCheck(hipsolverCreate(handle))

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

  ! Query workspace size and allocate it
  call hipsolverCheck(hipsolverDgeqrf_bufferSize(handle, M, N, dA, lda, lwork))
  call hipCheck(hipMalloc(dWork, int(lwork,c_size_t) * 8))

  ! Compute the QR factorization
  call hipsolverCheck(hipsolverDgeqrf(handle, M, N, dA, lda, dTau(1), dWork, lwork, dInfo(1)))

  ! Copy the factorized matrix back to host
  call hipCheck(hipMemcpy(hA, dA, hipMemcpyDeviceToHost))

  ! Extract R (upper triangle of the geqrf output)
  R = 0.0_c_double
  do j = 1,N
    do i = 1,j
      R(i,j) = hA(i,j)
    end do
  end do

  ! Verify A0**T * A0 = R**T * R
  do j = 1,N
    do i = 1,N
      lhs = 0.0_c_double
      rhs = 0.0_c_double
      do l = 1,M
        lhs = lhs + hA0(l,i) * hA0(l,j)
        rhs = rhs + R(l,i) * R(l,j)
      end do
      error = abs(lhs - rhs)
      if(error .gt. error_max) then
          write(*,*) "FAILED! Error bigger than max! Error = ", error, " (", i, ",", j, ")"
          call exit(1)
      end if
    end do
  end do

  ! Clean up
  call hipCheck(hipFree(dWork))
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dTau))
  call hipCheck(hipFree(dInfo))
  call hipsolverCheck(hipsolverDestroy(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.

!!!!!!!!!!!!!/
! hipsolverDorgqr example (double-precision generation of Q from a QR
! factorization)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Self-verifying: factorize A with geqrf, generate the orthogonal factor Q with
! orgqr, and confirm Q**T * Q = I. Orthogonality is sign-convention independent,
! so no reference matrix is needed.
!!!!!!!!!!!!!!/
!
program dorgqr
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver

  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 :: K = 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) :: hTau(3)
  real(c_double) :: gram

  type(c_ptr) :: handle = c_null_ptr
  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dTau(:)
  integer(c_int), pointer :: dInfo(:)
  type(c_ptr) :: dWork
  integer(c_int) :: lwork_qr, lwork_or, lwork

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

  call hipsolverCheck(hipsolverCreate(handle))

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

  ! Workspace big enough for both geqrf and orgqr
  call hipsolverCheck(hipsolverDgeqrf_bufferSize(handle, M, N, dA, lda, lwork_qr))
  call hipsolverCheck(hipsolverDorgqr_bufferSize(handle, M, N, K, dA, lda, dTau(1), lwork_or))
  lwork = max(lwork_qr, lwork_or)
  call hipCheck(hipMalloc(dWork, int(lwork,c_size_t) * 8))

  ! Factorize A = Q*R, then form the explicit Q in place
  call hipsolverCheck(hipsolverDgeqrf(handle, M, N, dA, lda, dTau(1), dWork, lwork, dInfo(1)))
  call hipsolverCheck(hipsolverDorgqr(handle, M, N, K, dA, lda, dTau(1), dWork, lwork, dInfo(1)))

  ! Copy Q back to host
  call hipCheck(hipMemcpy(hA, dA, hipMemcpyDeviceToHost))

  ! Verify Q**T * Q = I
  do j = 1,N
    do i = 1,N
      gram = sum(hA(:,i) * hA(:,j))
      if(i .eq. j) then
        error = abs(gram - 1.0_c_double)
      else
        error = abs(gram)
      end if
      if(error .gt. error_max) then
          write(*,*) "FAILED! Q not orthogonal! Error = ", error, " (", i, ",", j, ")"
          call exit(1)
      end if
    end do
  end do

  ! Clean up
  call hipCheck(hipFree(dWork))
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dTau))
  call hipCheck(hipFree(dInfo))
  call hipsolverCheck(hipsolverDestroy(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dorgqr

!!!!!!!!!!!!!/
! hipsolverDormqr example (double-precision multiply by Q from a QR
! factorization)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Self-verifying: factorize A with geqrf to obtain Q (as Householder vectors),
! then form Q*C with ormqr. Q is orthogonal, so it preserves the Frobenius norm:
! ||Q*C||_F = ||C||_F. Norm preservation is sign-convention independent, so no
! reference matrix is needed.
!!!!!!!!!!!!!!/
!
program dormqr
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver
  use hipfort_hipsolver_enums

  implicit none
  integer :: i, j

  integer(c_int), parameter :: M = 3
  integer(c_int), parameter :: N = 2   ! number of columns of C
  integer(c_int), parameter :: K = 3   ! number of reflectors
  integer(c_int), parameter :: lda = 3
  integer(c_int), parameter :: ldc = 3

  ! Matrix to factorize (column-major) and a separate C to multiply
  real(c_double) :: hA(3,3) = reshape((/1, 4, 7, 2, 5, 8, 3, 6, 10/), (/3, 3/))
  real(c_double) :: hC(3,2) = reshape((/1, 2, 3, 4, 5, 6/), (/3, 2/))
  real(c_double) :: hTau(3)
  real(c_double) :: norm_in, norm_out

  type(c_ptr) :: handle = c_null_ptr
  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dC(:,:)
  real(c_double), pointer :: dTau(:)
  integer(c_int), pointer :: dInfo(:)
  type(c_ptr) :: dWork
  integer(c_int) :: lwork_qr, lwork_mq, lwork

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

  ! Norm of the input C (Frobenius)
  norm_in = sqrt(sum(hC*hC))

  call hipsolverCheck(hipsolverCreate(handle))

  ! Allocate device-side memory & copy inputs to device
  call hipCheck(hipMalloc(dA, source=hA))
  call hipCheck(hipMalloc(dC, source=hC))
  call hipCheck(hipMalloc(dTau, mold=hTau))
  call hipCheck(hipMalloc(dInfo, 1))

  ! Workspace big enough for both geqrf and ormqr
  call hipsolverCheck(hipsolverDgeqrf_bufferSize(handle, M, K, dA, lda, lwork_qr))
  call hipsolverCheck(hipsolverDormqr_bufferSize(handle, HIPSOLVER_SIDE_LEFT, HIPSOLVER_OP_N, &
       M, N, K, dA, lda, dTau(1), dC, ldc, lwork_mq))
  lwork = max(lwork_qr, lwork_mq)
  call hipCheck(hipMalloc(dWork, int(lwork,c_size_t) * 8))

  ! Factorize A = Q*R (Q stored as reflectors), then form C <- Q*C in place
  call hipsolverCheck(hipsolverDgeqrf(handle, M, K, dA, lda, dTau(1), dWork, lwork, dInfo(1)))
  call hipsolverCheck(hipsolverDormqr(handle, HIPSOLVER_SIDE_LEFT, HIPSOLVER_OP_N, &
       M, N, K, dA, lda, dTau(1), dC, ldc, dWork, lwork, dInfo(1)))

  ! Copy the transformed C back to host
  call hipCheck(hipMemcpy(hC, dC, hipMemcpyDeviceToHost))

  ! Verify ||Q*C||_F = ||C||_F (Q is orthogonal)
  norm_out = sqrt(sum(hC*hC))
  error = abs(norm_out - norm_in) / max(norm_in, 1.0_c_double)
  if(error .gt. error_max) then
      write(*,*) "FAILED! Norm not preserved! ||C|| = ", norm_in, " ||Q*C|| = ", norm_out
      call exit(1)
  end if

  ! Clean up
  call hipCheck(hipFree(dWork))
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dC))
  call hipCheck(hipFree(dTau))
  call hipCheck(hipFree(dInfo))
  call hipsolverCheck(hipsolverDestroy(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dormqr

Symmetric eigenvalues#

syevd (heevd for Hermitian matrices) computes the eigenvalues, and optionally the eigenvectors, of a symmetric matrix with a divide-and-conquer algorithm. The HIPSOLVER_EIG_MODE_* argument selects whether eigenvectors are produced; the program requests eigenvalues only and checks their sum against the trace.

!!!!!!!!!!!!!!
! hipsolver dsyevd example (symmetric eigenvalues, Fortran 2008 interfaces)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Checks sum(eigenvalues) == trace(A). Native-array f2008 form; workspace via
! hipsolverDsyevd_bufferSize; devInfo device-backed.
!!!!!!!!!!!!!!
!
program hipsolver_dsyevd
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver
  implicit none
  integer(c_int), parameter :: N = 4, lda = 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) :: hD(N) = 0.0d0
  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dD(:)
  integer(c_int), pointer :: dInfo
  type(c_ptr) :: dWork, handle = c_null_ptr
  integer(c_int) :: lwork
  real(c_double) :: trace_A, error
  real(c_double), parameter :: rtol = 1.0d-9
  write(*,"(a)",advance="no") "-- Running test 'hipsolver_dsyevd' (Fortran 2008 interfaces) - "
  trace_A = hA(1,1) + hA(2,2) + hA(3,3) + hA(4,4)
  call hipsolverCheck(hipsolverCreate(handle))
  call hipCheck(hipMalloc(dA, source=hA))
  call hipCheck(hipMalloc(dD, source=hD))
  call hipCheck(hipMalloc(dInfo))
  call hipsolverCheck(hipsolverDsyevd_bufferSize(handle, HIPSOLVER_EIG_MODE_NOVECTOR, &
                                                 HIPSOLVER_FILL_MODE_UPPER, N, dA, lda, dD, lwork))
  call hipCheck(hipMalloc(dWork, max(int(lwork,c_size_t) * 8, 1_c_size_t)))
  call hipsolverCheck(hipsolverDsyevd(handle, HIPSOLVER_EIG_MODE_NOVECTOR, HIPSOLVER_FILL_MODE_UPPER, &
                                      N, dA, lda, dD, dWork, lwork, dInfo))
  call hipCheck(hipMemcpy(hD, dD, hipMemcpyDeviceToHost))
  error = abs(sum(hD) - trace_A) / abs(trace_A)
  if (error > rtol) then
     write(*,*) "FAILED! sum(eigenvalues) = ", sum(hD), " expected trace = ", trace_A
     call exit(1)
  end if
  call hipCheck(hipFree(dA)); call hipCheck(hipFree(dD)); call hipCheck(hipFree(dInfo)); call hipCheck(hipFree(dWork))
  call hipsolverCheck(hipsolverDestroy(handle))
  write(*,*) "PASSED!"
end program hipsolver_dsyevd

syevj/heevj solve the same problem with a Jacobi algorithm, which is often faster for small matrices. With HIPSOLVER_EIG_MODE_VECTOR the matrix is overwritten with the eigenvectors; the program confirms each eigenpair satisfies A*v = lambda*v.

!!!!!!!!!!!!!/
! hipsolverDsyevj example (double-precision Jacobi symmetric eigensolver)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Self-verifying: with jobz=vector, A is overwritten with the eigenvectors (as
! columns) and W holds the eigenvalues. We confirm A0*v_k = lambda_k*v_k.
!
! syevj has no generated array-pointer overload, so A and W are passed to the
! generic interface via c_loc of their device pointers.
!!!!!!!!!!!!!!/
!
program dsyevj
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver
  use hipfort_hipsolver_enums

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

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

  ! Symmetric input (column-major); eigenvalues are 2-sqrt(2), 2, 2+sqrt(2)
  real(c_double) :: hA(3,3) = reshape((/2, -1, 0, -1, 2, -1, 0, -1, 2/), (/3, 3/))
  real(c_double) :: hA0(3,3)   ! original kept for verification
  real(c_double) :: hW(3)      ! eigenvalues
  real(c_double) :: lhs(3), rhs(3)

  type(c_ptr) :: handle = c_null_ptr
  type(c_ptr) :: params = c_null_ptr
  real(c_double), pointer :: dA(:,:)
  real(c_double), pointer :: dW(:)
  integer(c_int), pointer :: dInfo(:)
  type(c_ptr) :: dWork
  integer(c_int) :: lwork

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

  hA0 = hA ! keep original for the A*v = lambda*v check

  call hipsolverCheck(hipsolverCreate(handle))
  call hipsolverCheck(hipsolverCreateSyevjInfo(params))

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

  ! Query workspace size and allocate it
  call hipsolverCheck(hipsolverDsyevj_bufferSize(handle, HIPSOLVER_EIG_MODE_VECTOR, &
       HIPSOLVER_FILL_MODE_UPPER, N, c_loc(dA(1,1)), lda, c_loc(dW(1)), lwork, params))
  call hipCheck(hipMalloc(dWork, int(lwork,c_size_t) * 8))

  ! Compute eigenvalues and eigenvectors (A overwritten with eigenvectors)
  call hipsolverCheck(hipsolverDsyevj(handle, HIPSOLVER_EIG_MODE_VECTOR, &
       HIPSOLVER_FILL_MODE_UPPER, N, c_loc(dA(1,1)), lda, c_loc(dW(1)), dWork, lwork, dInfo(1), params))

  ! Copy results back to host
  call hipCheck(hipMemcpy(hA, dA, hipMemcpyDeviceToHost))
  call hipCheck(hipMemcpy(hW, dW, hipMemcpyDeviceToHost))

  ! Verify A0 * v_k = lambda_k * v_k for each eigenpair
  do k = 1,N
    lhs = matmul(hA0, hA(:,k))
    rhs = hW(k) * hA(:,k)
    do i = 1,N
        error = abs(lhs(i) - rhs(i))
        if(error .gt. error_max) then
            write(*,*) "FAILED! Error bigger than max! Error = ", error, " eigenpair ", k
            call exit(1)
        end if
    end do
  end do

  ! Clean up
  call hipCheck(hipFree(dWork))
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dW))
  call hipCheck(hipFree(dInfo))
  call hipsolverCheck(hipsolverDestroySyevjInfo(params))
  call hipsolverCheck(hipsolverDestroy(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dsyevj

Singular value decomposition#

gesvd computes the singular value decomposition A = U*S*V**T. The character(c_char) job codes choose which singular-vector matrices are computed. The program requests singular values only ('N') and checks the convention-independent invariant sum(sigma_i**2) == ||A||_F**2.

!!!!!!!!!!!!!!
! hipsolver dgesvd example (singular value decomposition)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Computes the singular values of A (jobu = jobv = 'N', values only) and checks
! the convention-independent invariant sum(sigma_i^2) == ||A||_F^2.
!
! Note: jobu/jobv are `signed char` job codes passed by value; the hipfort
! binding now types them as character(c_char) (they were previously type(c_ptr),
! which made this routine uncallable). devInfo lives on the DEVICE.
!!!!!!!!!!!!!!
!
program hipsolver_dgesvd
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver

  implicit none

  integer(c_int), parameter :: M = 2, N = 2, lda = 2, ldu = 2, ldv = 2
  integer(c_int), parameter :: mn = 2   ! min(M,N)

  ! A = [[1,2],[3,4]] (column-major); ||A||_F^2 = 1+9+4+16 = 30.
  real(c_double), target :: hA(M,N) = reshape((/1.0d0, 3.0d0, 2.0d0, 4.0d0/), (/M,N/))
  real(c_double), target :: hS(mn)

  type(c_ptr) :: dA, dS, dU, dV, dWork, dRwork, dInfo, handle = c_null_ptr
  integer(c_int) :: lwork
  integer(c_size_t) :: szA = M*N, szS = mn, szU = M*M, szV = N*N, szR = mn
  real(c_double) :: frob, ssum, error
  real(c_double), parameter :: rtol = 1.0d-9
  integer :: i

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

  call hipsolverCheck(hipsolverCreate(handle))
  call hipCheck(hipMalloc(dA,    szA * 8))
  call hipCheck(hipMalloc(dS,    szS * 8))
  call hipCheck(hipMalloc(dU,    szU * 8))
  call hipCheck(hipMalloc(dV,    szV * 8))
  call hipCheck(hipMalloc(dRwork, szR * 8))
  call hipCheck(hipMalloc(dInfo, 4_c_size_t))
  call hipCheck(hipMemcpy(dA, c_loc(hA(1,1)), szA * 8, hipMemcpyHostToDevice))

  call hipsolverCheck(hipsolverDgesvd_bufferSize(handle, 'N', 'N', M, N, lwork))
  call hipCheck(hipMalloc(dWork, max(int(lwork,c_size_t) * 8, 1_c_size_t)))

  ! Singular values only (jobu = jobv = 'N').
  call hipsolverCheck(hipsolverDgesvd(handle, 'N', 'N', M, N, dA, lda, dS, &
                                      dU, ldu, dV, ldv, dWork, lwork, dRwork, dInfo))

  call hipCheck(hipMemcpy(c_loc(hS(1)), dS, szS * 8, hipMemcpyDeviceToHost))

  frob = 1.0d0 + 9.0d0 + 4.0d0 + 16.0d0   ! ||A||_F^2
  ssum = 0.0d0
  do i = 1, mn
     ssum = ssum + hS(i)**2
  end do
  error = abs(ssum - frob) / frob
  if (error > rtol) then
     write(*,*) "FAILED! sum(sigma^2) = ", ssum, " expected ||A||_F^2 = ", frob
     call exit(1)
  end if

  call hipCheck(hipFree(dA)); call hipCheck(hipFree(dS)); call hipCheck(hipFree(dU))
  call hipCheck(hipFree(dV)); call hipCheck(hipFree(dRwork)); call hipCheck(hipFree(dWork))
  call hipCheck(hipFree(dInfo))
  call hipsolverCheck(hipsolverDestroy(handle))

  write(*,*) "PASSED!"

end program hipsolver_dgesvd

gesvdj computes the same decomposition with a Jacobi algorithm. The program requests all vectors and reconstructs A from the factors, which avoids the sign and order ambiguity of the singular vectors.

!!!!!!!!!!!!!/
! hipsolverDgesvdj example (double-precision Jacobi SVD)
! see: https:!rocm.docs.amd.com/projects/hipSOLVER/en/latest/
!
! Self-verifying: compute the singular value decomposition A = U*S*V**T with the
! Jacobi method and confirm the factors reconstruct the original matrix. Using
! the reconstruction avoids sign/order ambiguity in the singular vectors.
!
! gesvdj has no generated array-pointer overload, so the matrix/vector arguments
! are passed to the generic interface via c_loc of their device pointers.
!!!!!!!!!!!!!!/
!
program dgesvdj
  use iso_c_binding
  use hipfort
  use hipfort_check
  use hipfort_hipsolver
  use hipfort_hipsolver_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
  integer(c_int), parameter :: econ = 0

  ! Nonsingular 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) :: hA0(3,3)          ! original kept for verification
  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
  real(c_double) :: recon(3,3)

  type(c_ptr) :: handle = c_null_ptr
  type(c_ptr) :: params = c_null_ptr
  real(c_double), pointer :: dA(:,:), dU(:,:), dV(:,:)
  real(c_double), pointer :: dS(:)
  integer(c_int), pointer :: dInfo(:)
  type(c_ptr) :: dWork
  integer(c_int) :: lwork

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

  hA0 = hA ! keep original for the reconstruction check

  call hipsolverCheck(hipsolverCreate(handle))
  call hipsolverCheck(hipsolverCreateGesvdjInfo(params))

  ! 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(dInfo, 1))

  ! Query workspace size and allocate it
  call hipsolverCheck(hipsolverDgesvdj_bufferSize(handle, HIPSOLVER_EIG_MODE_VECTOR, econ, &
       M, N, c_loc(dA(1,1)), lda, c_loc(dS(1)), c_loc(dU(1,1)), ldu, c_loc(dV(1,1)), ldv, lwork, params))
  call hipCheck(hipMalloc(dWork, int(lwork,c_size_t) * 8))

  ! Compute the singular value decomposition
  call hipsolverCheck(hipsolverDgesvdj(handle, HIPSOLVER_EIG_MODE_VECTOR, econ, &
       M, N, c_loc(dA(1,1)), lda, c_loc(dS(1)), c_loc(dU(1,1)), ldu, c_loc(dV(1,1)), ldv, &
       dWork, lwork, dInfo(1), params))

  ! Copy results 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**T and compare with the original
  do j = 1,N
    do i = 1,M
      recon(i,j) = sum(hU(i,:) * hS(:) * hV(j,:))
    end do
  end do

  do j = 1,N
    do i = 1,M
      error = abs(recon(i,j) - hA0(i,j))
      if(error .gt. error_max) then
          write(*,*) "FAILED! Error bigger than max! Error = ", error, " (", i, ",", j, ")"
          call exit(1)
      end if
    end do
  end do

  ! Clean up
  call hipCheck(hipFree(dWork))
  call hipCheck(hipFree(dA))
  call hipCheck(hipFree(dS))
  call hipCheck(hipFree(dU))
  call hipCheck(hipFree(dV))
  call hipCheck(hipFree(dInfo))
  call hipsolverCheck(hipsolverDestroyGesvdjInfo(params))
  call hipsolverCheck(hipsolverDestroy(handle))
  call hipCheck(hipDeviceReset())

  write(*,*) "PASSED!"

end program dgesvdj
Contents