1! ----------------------------------------------------------------- 2! Programmer(s): Cody J. Balos @ LLNL 3! ----------------------------------------------------------------- 4! SUNDIALS Copyright Start 5! Copyright (c) 2002-2021, Lawrence Livermore National Security 6! and Southern Methodist University. 7! All rights reserved. 8! 9! See the top-level LICENSE and NOTICE files for details. 10! 11! SPDX-License-Identifier: BSD-3-Clause 12! SUNDIALS Copyright End 13! ----------------------------------------------------------------- 14! This file tests the Fortran 2003 interface to the SUNDIALS 15! band SUNLinearSolver implementation. 16! ----------------------------------------------------------------- 17 18module test_fsunlinsol_band 19 use, intrinsic :: iso_c_binding 20 implicit none 21 22 integer(C_LONG), parameter :: N = 10 23 integer(C_LONG), parameter :: mu = 2 24 integer(C_LONG), parameter :: ml = 3 25 26contains 27 28 integer(C_INT) function unit_tests() result(fails) 29 use, intrinsic :: iso_c_binding 30 use fsundials_nvector_mod 31 use fsundials_matrix_mod 32 use fsundials_linearsolver_mod 33 use fnvector_serial_mod 34 use fsunmatrix_band_mod 35 use fsunlinsol_band_mod 36 use test_sunlinsol 37 38 implicit none 39 40 type(SUNLinearSolver), pointer :: LS ! test linear solver 41 type(SUNMatrix), pointer :: A ! test matrices 42 type(N_Vector), pointer :: x, y, b ! test vectors 43 real(C_DOUBLE), pointer :: colj(:) ! matrix column data 44 real(C_DOUBLE), pointer :: xdata(:), Adata(:) ! data arrays 45 real(C_DOUBLE) :: tmpr ! temporary real value 46 integer(C_LONG) :: j, k 47 integer(C_LONG) :: smu, kstart, kend, offset 48 integer(C_INT) :: tmp 49 50 fails = 0 51 smu = 0 52 53 A => FSUNBandMatrix(N, mu, ml) 54 x => FN_VNew_Serial(N) 55 y => FN_VNew_Serial(N) 56 b => FN_VNew_Serial(N) 57 58 ! fill A matrix with uniform random data in [0, 1/N) 59 Adata => FSUNBandMatrix_Data(A) 60 do j=1, N 61 offset = smu-mu+1 + j*(smu+ml+1) 62 colj(-mu:ml) => Adata(offset:mu+ml) 63 kstart = merge(-j, -mu, j < mu) 64 kend = merge(N-1-j, ml, j > N-1-ml) 65 do k=kstart, kend 66 call random_number(tmpr) 67 colj(k) = tmpr / N 68 end do 69 end do 70 71 ! fill x vector with uniform random data in [1, 2) 72 xdata => FN_VGetArrayPointer(x) 73 do j=1, N 74 call random_number(tmpr) 75 xdata(j) = ONE + tmpr 76 end do 77 78 ! scale/shift matrix to ensure diagonal dominance 79 fails = FSUNMatScaleAddI(ONE/(mu+ml+1), A) 80 if (fails /= 0) then 81 print *, 'FAIL: FSUNMatScaleAddI failure' 82 call FSUNMatDestroy(A) 83 call FN_VDestroy(x) 84 call FN_VDestroy(y) 85 call FN_VDestroy(b) 86 return 87 end if 88 89 ! create RHS vector for linear solve 90 fails = FSUNMatMatvec(A, x, b) 91 if (fails /= 0) then 92 print *, 'FAIL: FSUNMatMatvec failure' 93 call FSUNMatDestroy(A) 94 call FN_VDestroy(x) 95 call FN_VDestroy(y) 96 call FN_VDestroy(b) 97 return 98 end if 99 100 ! create band linear solver 101 LS => FSUNLinSol_Band(x, A) 102 103 ! run tests 104 fails = fails + Test_FSUNLinSolInitialize(LS, 0) 105 fails = fails + Test_FSUNLinSolSetup(LS, A, 0) 106 fails = fails + Test_FSUNLinSolSolve(LS, A, x, b, 100*UNIT_ROUNDOFF, 0) 107 108 fails = fails + Test_FSUNLinSolGetType(LS, SUNLINEARSOLVER_DIRECT, 0) 109 fails = fails + Test_FSUNLinSolLastFlag(LS, 0) 110 fails = fails + Test_FSUNLinSolSpace(LS, 0) 111 112 ! cleanup 113 tmp = FSUNLinSolFree(LS) 114 call FSUNMatDestroy(A) 115 call FN_VDestroy(x) 116 call FN_VDestroy(y) 117 call FN_VDestroy(b) 118 119 end function unit_tests 120 121end module 122 123integer(C_INT) function check_vector(X, Y, tol) result(failure) 124 use, intrinsic :: iso_c_binding 125 use fsundials_nvector_mod 126 use test_utilities 127 128 implicit none 129 type(N_Vector) :: x, y 130 real(C_DOUBLE) :: tol, maxerr 131 integer(C_LONG) :: i, xlen, ylen 132 real(C_DOUBLE), pointer :: xdata(:), ydata(:) 133 134 failure = 0 135 136 xdata => FN_VGetArrayPointer(x) 137 ydata => FN_VGetArrayPointer(y) 138 139 xlen = FN_VGetLength(x) 140 ylen = FN_VGetLength(y) 141 142 if (xlen /= ylen) then 143 print *, 'FAIL: check_vector: different data array lengths' 144 failure = 1 145 return 146 end if 147 148 do i = 1, xlen 149 failure = failure + FNEQTOL(xdata(i), ydata(i), FIVE*tol*abs(xdata(i))) 150 end do 151 152 if (failure > 0) then 153 maxerr = ZERO 154 do i = 1, xlen 155 maxerr = max(abs(xdata(i)-ydata(i))/abs(xdata(i)), maxerr) 156 end do 157 write(*,'(A,E14.7,A,E14.7,A)') & 158 "FAIL: check_vector failure: maxerr = ", maxerr, " (tol = ", FIVE*tol, ")" 159 end if 160 161end function check_vector 162 163program main 164 !======== Inclusions ========== 165 use, intrinsic :: iso_c_binding 166 use test_fsunlinsol_band 167 168 !======== Declarations ======== 169 implicit none 170 integer(C_INT) :: fails = 0 171 172 !============== Introduction ============= 173 print *, 'Band SUNLinearSolver Fortran 2003 interface test' 174 175 fails = unit_tests() 176 if (fails /= 0) then 177 print *, 'FAILURE: n unit tests failed' 178 stop 1 179 else 180 print *,'SUCCESS: all unit tests passed' 181 end if 182end program main 183