From c25e3ff248526edffb632e0d4f317da8e13ed3d4 Mon Sep 17 00:00:00 2001
From: Federico Perini <federico.perini@gmail.com>
Date: Mon, 11 Mar 2024 17:24:07 +0100
Subject: [PATCH] add tests

---
 test/linalg/CMakeLists.txt        |   2 +
 test/linalg/test_blas_lapack.fypp | 140 ++++++++++++++++++++++++++++++
 2 files changed, 142 insertions(+)
 create mode 100644 test/linalg/test_blas_lapack.fypp

diff --git a/test/linalg/CMakeLists.txt b/test/linalg/CMakeLists.txt
index 4a315f545..3d590a9d2 100644
--- a/test/linalg/CMakeLists.txt
+++ b/test/linalg/CMakeLists.txt
@@ -1,9 +1,11 @@
 set(
   fppFiles
   "test_linalg.fypp"
+  "test_blas_lapack.fypp"
   "test_linalg_matrix_property_checks.fypp"
 )
 fypp_f90("${fyppFlags}" "${fppFiles}" outFiles)
 
 ADDTEST(linalg)
 ADDTEST(linalg_matrix_property_checks)
+ADDTEST(blas_lapack)
diff --git a/test/linalg/test_blas_lapack.fypp b/test/linalg/test_blas_lapack.fypp
new file mode 100644
index 000000000..0e9d338a8
--- /dev/null
+++ b/test/linalg/test_blas_lapack.fypp
@@ -0,0 +1,140 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+
+module test_blas_lapack
+    use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test
+    use stdlib_kinds, only: sp, dp, xdp, qp, int8, int16, int32, int64
+    use stdlib_linalg, only: eye
+    use stdlib_linalg_blas
+    use stdlib_linalg_lapack
+
+    implicit none
+
+    real(sp), parameter :: sptol = 1000 * epsilon(1._sp)
+    real(dp), parameter :: dptol = 1000 * epsilon(1._dp)
+#:if WITH_QP
+    real(qp), parameter :: qptol = 1000 * epsilon(1._qp)
+#:endif
+
+
+
+contains
+
+    !> Collect all exported unit tests
+    subroutine collect_blas_lapack(testsuite)
+        !> Collection of tests
+        type(unittest_type), allocatable, intent(out) :: testsuite(:)
+
+        testsuite = [ &
+  #:for k1, t1 in REAL_KINDS_TYPES
+            new_unittest("test_gemv${t1[0]}$${k1}$", test_gemv${t1[0]}$${k1}$), &
+            new_unittest("test_getri${t1[0]}$${k1}$", test_gemv${t1[0]}$${k1}$), &
+  #:endfor
+            new_unittest("test_idamax", test_idamax) &
+            ]
+
+    end subroutine collect_blas_lapack
+
+
+  #:for k1, t1 in REAL_KINDS_TYPES
+    subroutine test_gemv${t1[0]}$${k1}$(error)
+        !> Error handling
+        type(error_type), allocatable, intent(out) :: error
+
+        ${t1}$ :: A(3,3),x(3),y(3),ylap(3),yintr(3),alpha,beta
+        call random_number(alpha)
+        call random_number(beta)
+        call random_number(A)
+        call random_number(x)
+        call random_number(y)
+        ylap = y
+        call gemv('No transpose',size(A,1),size(A,2),alpha,A,size(A,1),x,1,beta,ylap,1)
+        yintr = alpha*matmul(A,x)+beta*y
+
+        call check(error, sum(abs(ylap - yintr)) < sptol, &
+            "blas vs. intrinsics axpy: sum() < sptol failed")
+        if (allocated(error)) return
+
+    end subroutine test_gemv${t1[0]}$${k1}$
+
+    ! Find matrix inverse from LU decomposition
+    subroutine test_getri${t1[0]}$${k1}$(error)
+        !> Error handling
+        type(error_type), allocatable, intent(out) :: error
+
+        integer(ilp), parameter :: n = 3
+        ${t1}$ :: A(n,n)
+        ${t1}$,allocatable :: work(:)
+        integer(ilp) :: ipiv(n),info,lwork,nb
+
+
+        A = eye(n)
+
+        ! Factorize matrix (overwrite result)
+        call getrf(size(A,1),size(A,2),A,size(A,1),ipiv,info)
+        call check(error, info==0, "lapack getrf returned info/=0")
+        if (allocated(error)) return
+
+        ! Get optimal worksize (returned in work(1)) (apply 2% safety parameter)
+        nb = stdlib_ilaenv(1,'${t1[0]}$getri',' ',n,-1,-1,-1)
+        lwork = nint(1.02*n*nb,kind=ilp)
+        allocate (work(lwork))
+
+        ! Invert matrix
+        call getri(n,a,n,ipiv,work,lwork,info)
+
+        call check(error, info==0, "lapack getri returned info/=0")
+        if (allocated(error)) return
+
+        call check(error, sum(abs(A - eye(3))) < sptol, &
+            "lapack eye inversion: tolerance check failed")
+        if (allocated(error)) return
+
+    end subroutine test_getri${t1[0]}$${k1}$
+  #:endfor
+
+    ! Return
+    subroutine test_idamax(error)
+        !> Error handling
+        type(error_type), allocatable, intent(out) :: error
+
+        integer(ilp), parameter :: n = 5
+        integer(ilp) :: imax
+        real(dp) :: x(n)
+
+        x = [1,2,3,4,5]
+
+        imax = stdlib_idamax(n,x,1)
+
+        call check(error, imax==5, "blas idamax returned wrong location")
+
+    end subroutine test_idamax
+
+end module test_blas_lapack
+
+
+program tester
+    use, intrinsic :: iso_fortran_env, only : error_unit
+    use testdrive, only : run_testsuite, new_testsuite, testsuite_type
+    use test_blas_lapack, only : collect_blas_lapack
+    implicit none
+    integer :: stat, is
+    type(testsuite_type), allocatable :: testsuites(:)
+    character(len=*), parameter :: fmt = '("#", *(1x, a))'
+
+    stat = 0
+    testsuites = [ &
+        new_testsuite("blas_lapack", collect_blas_lapack) &
+        ]
+
+    do is = 1, size(testsuites)
+        write(error_unit, fmt) "Testing:", testsuites(is)%name
+        call run_testsuite(testsuites(is)%collect, error_unit, stat)
+    end do
+
+    if (stat > 0) then
+        write(error_unit, '(i0, 1x, a)') stat, "test(s) failed!"
+        error stop
+    end if
+end program
+