|
| 1 | +#:include "common.fypp" |
| 2 | +#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES |
| 3 | + |
| 4 | +module test_blas_lapack |
| 5 | + use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test |
| 6 | + use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64 |
| 7 | + use stdlib_linalg, only: eye |
| 8 | + use stdlib_linalg_blas |
| 9 | + use stdlib_linalg_lapack |
| 10 | + |
| 11 | + implicit none |
| 12 | + |
| 13 | + real(sp), parameter :: sptol = 1000 * epsilon(1._sp) |
| 14 | + real(dp), parameter :: dptol = 1000 * epsilon(1._dp) |
| 15 | +#:if WITH_QP |
| 16 | + real(qp), parameter :: qptol = 1000 * epsilon(1._qp) |
| 17 | +#:endif |
| 18 | + |
| 19 | + |
| 20 | + |
| 21 | +contains |
| 22 | + |
| 23 | + !> Collect all exported unit tests |
| 24 | + subroutine collect_blas_lapack(testsuite) |
| 25 | + !> Collection of tests |
| 26 | + type(unittest_type), allocatable, intent(out) :: testsuite(:) |
| 27 | + |
| 28 | + testsuite = [ & |
| 29 | + #:for k1, t1 in REAL_KINDS_TYPES |
| 30 | + new_unittest("test_gemv${t1[0]}$${k1}$", test_gemv${t1[0]}$${k1}$), & |
| 31 | + new_unittest("test_getri${t1[0]}$${k1}$", test_gemv${t1[0]}$${k1}$), & |
| 32 | + #:endfor |
| 33 | + new_unittest("test_idamax", test_idamax) & |
| 34 | + ] |
| 35 | + |
| 36 | + end subroutine collect_blas_lapack |
| 37 | + |
| 38 | + |
| 39 | + #:for k1, t1 in REAL_KINDS_TYPES |
| 40 | + subroutine test_gemv${t1[0]}$${k1}$(error) |
| 41 | + !> Error handling |
| 42 | + type(error_type), allocatable, intent(out) :: error |
| 43 | + |
| 44 | + ${t1}$ :: A(3,3),x(3),y(3),ylap(3),yintr(3),alpha,beta |
| 45 | + call random_number(alpha) |
| 46 | + call random_number(beta) |
| 47 | + call random_number(A) |
| 48 | + call random_number(x) |
| 49 | + call random_number(y) |
| 50 | + ylap = y |
| 51 | + call gemv('No transpose',size(A,1),size(A,2),alpha,A,size(A,1),x,1,beta,ylap,1) |
| 52 | + yintr = alpha*matmul(A,x)+beta*y |
| 53 | + |
| 54 | + call check(error, sum(abs(ylap - yintr)) < sptol, & |
| 55 | + "blas vs. intrinsics axpy: sum() < sptol failed") |
| 56 | + if (allocated(error)) return |
| 57 | + |
| 58 | + end subroutine test_gemv${t1[0]}$${k1}$ |
| 59 | + |
| 60 | + ! Find matrix inverse from LU decomposition |
| 61 | + subroutine test_getri${t1[0]}$${k1}$(error) |
| 62 | + !> Error handling |
| 63 | + type(error_type), allocatable, intent(out) :: error |
| 64 | + |
| 65 | + integer(ilp), parameter :: n = 3 |
| 66 | + ${t1}$ :: A(n,n) |
| 67 | + ${t1}$,allocatable :: work(:) |
| 68 | + integer(ilp) :: ipiv(n),info,lwork,nb |
| 69 | + |
| 70 | + |
| 71 | + A = eye(n) |
| 72 | + |
| 73 | + ! Factorize matrix (overwrite result) |
| 74 | + call getrf(size(A,1),size(A,2),A,size(A,1),ipiv,info) |
| 75 | + call check(error, info==0, "lapack getrf returned info/=0") |
| 76 | + if (allocated(error)) return |
| 77 | + |
| 78 | + ! Get optimal worksize (returned in work(1)) (apply 2% safety parameter) |
| 79 | + nb = stdlib_ilaenv(1,'${t1[0]}$getri',' ',n,-1,-1,-1) |
| 80 | + lwork = nint(1.02*n*nb,kind=ilp) |
| 81 | + allocate (work(lwork)) |
| 82 | + |
| 83 | + ! Invert matrix |
| 84 | + call getri(n,a,n,ipiv,work,lwork,info) |
| 85 | + |
| 86 | + call check(error, info==0, "lapack getri returned info/=0") |
| 87 | + if (allocated(error)) return |
| 88 | + |
| 89 | + call check(error, sum(abs(A - eye(3))) < sptol, & |
| 90 | + "lapack eye inversion: tolerance check failed") |
| 91 | + if (allocated(error)) return |
| 92 | + |
| 93 | + end subroutine test_getri${t1[0]}$${k1}$ |
| 94 | + #:endfor |
| 95 | + |
| 96 | + ! Return |
| 97 | + subroutine test_idamax(error) |
| 98 | + !> Error handling |
| 99 | + type(error_type), allocatable, intent(out) :: error |
| 100 | + |
| 101 | + integer(ilp), parameter :: n = 5 |
| 102 | + integer(ilp) :: imax |
| 103 | + real(dp) :: x(n) |
| 104 | + |
| 105 | + x = [1,2,3,4,5] |
| 106 | + |
| 107 | + imax = stdlib_idamax(n,x,1) |
| 108 | + |
| 109 | + call check(error, imax==5, "blas idamax returned wrong location") |
| 110 | + |
| 111 | + end subroutine test_idamax |
| 112 | + |
| 113 | +end module test_blas_lapack |
| 114 | + |
| 115 | + |
| 116 | +program tester |
| 117 | + use, intrinsic :: iso_fortran_env, only : error_unit |
| 118 | + use testdrive, only : run_testsuite, new_testsuite, testsuite_type |
| 119 | + use test_blas_lapack, only : collect_blas_lapack |
| 120 | + implicit none |
| 121 | + integer :: stat, is |
| 122 | + type(testsuite_type), allocatable :: testsuites(:) |
| 123 | + character(len=*), parameter :: fmt = '("#", *(1x, a))' |
| 124 | + |
| 125 | + stat = 0 |
| 126 | + testsuites = [ & |
| 127 | + new_testsuite("blas_lapack", collect_blas_lapack) & |
| 128 | + ] |
| 129 | + |
| 130 | + do is = 1, size(testsuites) |
| 131 | + write(error_unit, fmt) "Testing:", testsuites(is)%name |
| 132 | + call run_testsuite(testsuites(is)%collect, error_unit, stat) |
| 133 | + end do |
| 134 | + |
| 135 | + if (stat > 0) then |
| 136 | + write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!" |
| 137 | + error stop |
| 138 | + end if |
| 139 | +end program |
| 140 | + |
0 commit comments