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