1! ------------------------------------------------------------------
2! Programmer(s): Daniel R. Reynolds @ SMU
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! Program to test custom fnvector_fortran_mod implementation
15! ------------------------------------------------------------------
16
17! ------------------------------------------------------------------
18! Utility module for error-checking
19! ------------------------------------------------------------------
20module fnvector_test_mod
21  use, intrinsic :: iso_c_binding
22  use fnvector_fortran_mod
23  implicit none
24
25contains
26  integer(c_int) function check_ans(val, tol, Nvar, N, sunvec_x) result(failure)
27
28    implicit none
29    real(c_double), value :: val, tol
30    integer(c_long), value :: Nvar, N
31    Type(N_Vector) :: sunvec_x
32    Type(FVec), pointer :: x
33    integer(c_long) :: i, j
34
35    x => FN_VGetFVec(sunvec_x)
36    failure = 0
37    do j = 1,N
38       do i = 1,Nvar
39          if (dabs(x%data(i,j) - val) > tol)  failure = 1
40       end do
41    end do
42
43  end function check_ans
44end module fnvector_test_mod
45
46! ------------------------------------------------------------------
47program main
48
49  !======= Inclusions ===========
50  use, intrinsic :: iso_c_binding
51  use fnvector_fortran_mod
52  use fnvector_test_mod
53
54  !======= Declarations =========
55  implicit none
56
57  ! local variables
58  integer(c_int)             :: fails
59  integer(c_long)            :: i, j, loc
60  integer(c_long), parameter :: N = 1000
61  integer(c_long), parameter :: Nvar = 10
62  type(N_Vector), pointer    :: sU, sV, sW, sX, sY, sZ
63  type(FVec), pointer        :: U, V, W, X, Y, Z
64  real(c_double)             :: Udata(Nvar,N)
65  real(c_double)             :: fac
66  logical                    :: failure
67
68
69  !======= Internals ============
70
71  ! initialize failure total
72  fails = 0
73
74  ! create new vectors, using New, Make and Clone routines
75  sU => FN_VMake_Fortran(Nvar, N, Udata)
76  if (.not. associated(sU)) then
77     print *, 'ERROR: sunvec = NULL'
78     stop 1
79  end if
80  U => FN_VGetFVec(sU)
81
82  sV => FN_VNew_Fortran(Nvar, N)
83  if (.not. associated(sV)) then
84     print *, 'ERROR: sunvec = NULL'
85     stop 1
86  end if
87  V => FN_VGetFVec(sV)
88
89  sW => FN_VNew_Fortran(Nvar, N)
90  if (.not. associated(sW)) then
91     print *, 'ERROR: sunvec = NULL'
92     stop 1
93  end if
94  W => FN_VGetFVec(sW)
95
96  sX => FN_VNew_Fortran(Nvar, N)
97  if (.not. associated(sX)) then
98     print *, 'ERROR: sunvec = NULL'
99     stop 1
100  end if
101  X => FN_VGetFVec(sX)
102
103  sY => FN_VNew_Fortran(Nvar, N)
104  if (.not. associated(sY)) then
105     print *, 'ERROR: sunvec = NULL'
106     stop 1
107  end if
108  Y => FN_VGetFVec(sY)
109
110  call c_f_pointer(FN_VClone_Fortran(sU), sZ)
111  if (.not. associated(sZ)) then
112     print *, 'ERROR: sunvec = NULL'
113     stop 1
114  end if
115  Z => FN_VGetFVec(sZ)
116
117
118  ! check vector ID
119  if (FN_VGetVectorID(sU) /= SUNDIALS_NVEC_CUSTOM) then
120     fails = fails + 1
121     print *, '>>> FAILED test -- FN_VGetVectorID'
122     print *, '    Unrecognized vector type', FN_VGetVectorID(sU)
123  else
124     print *, 'PASSED test -- FN_VGetVectorID'
125  end if
126
127
128  ! check vector length
129  if (FN_VGetLength(sV) /= (N*Nvar)) then
130     fails = fails + 1
131     print *, '>>> FAILED test -- FN_VGetLength'
132     print *, '    ', FN_VGetLength(sV), ' /= ', N*Nvar
133  else
134     print *, 'PASSED test -- FN_VGetLength'
135  end if
136
137  ! test FN_VConst
138  Udata = 0.d0
139  call FN_VConst(1.d0, sU)
140  if (check_ans(1.d0, 1.d-14, Nvar, N, sU) /= 0) then
141     fails = fails + 1
142     print *, '>>> FAILED test -- FN_VConst'
143  else
144     print *, 'PASSED test -- FN_VConst'
145  end if
146
147  ! test FN_VLinearSum
148  call FN_VConst(1.d0, sX)
149  call FN_VConst(-2.d0, sY)
150  call FN_VLinearSum(1.d0, sX, 1.d0, sY, sY)
151  if (check_ans(-1.d0, 1.d-14, Nvar, N, sY) /= 0) then
152     fails = fails + 1
153     print *, '>>> FAILED test -- FN_VLinearSum Case 1a'
154  else
155     print *, 'PASSED test -- FN_VLinearSum Case 1a'
156  end if
157
158  call FN_VConst(1.d0, sX)
159  call FN_VConst(2.d0, sY)
160  call FN_VLinearSum(-1.d0, sX, 1.d0, sY, sY)
161  if (check_ans(1.d0, 1.d-14, Nvar, N, sY) /= 0) then
162     fails = fails + 1
163     print *, '>>> FAILED test -- FN_VLinearSum Case 1b'
164  else
165     print *, 'PASSED test -- FN_VLinearSum Case 1b'
166  end if
167
168  call FN_VConst(2.d0, sX)
169  call FN_VConst(-2.d0, sY)
170  call FN_VLinearSum(0.5d0, sX, 1.d0, sY, sY)
171  if (check_ans(-1.d0, 1.d-14, Nvar, N, sY) /= 0) then
172     fails = fails + 1
173     print *, '>>> FAILED test -- FN_VLinearSum Case 1c'
174  else
175     print *, 'PASSED test -- FN_VLinearSum Case 1c'
176  end if
177
178  call FN_VConst(2.d0, sX)
179  call FN_VConst(-1.d0, sY)
180  call FN_VLinearSum(1.d0, sX, 1.d0, sY, sX)
181  if (check_ans(1.d0, 1.d-14, Nvar, N, sX) /= 0) then
182     fails = fails + 1
183     print *, '>>> FAILED test -- FN_VLinearSum Case 2a'
184  else
185     print *, 'PASSED test -- FN_VLinearSum Case 2a'
186  end if
187
188  call FN_VConst(1.d0, sX)
189  call FN_VConst(2.d0, sY)
190  call FN_VLinearSum(1.d0, sX, -1.d0, sY, sX)
191  if (check_ans(-1.d0, 1.d-14, Nvar, N, sX) /= 0) then
192     fails = fails + 1
193     print *, '>>> FAILED test -- FN_VLinearSum Case 2b'
194  else
195     print *, 'PASSED test -- FN_VLinearSum Case 2b'
196  end if
197
198  call FN_VConst(2.d0, sX)
199  call FN_VConst(-0.5d0, sY)
200  call FN_VLinearSum(1.d0, sX, 2.d0, sY, sX)
201  if (check_ans(1.d0, 1.d-14, Nvar, N, sX) /= 0) then
202     fails = fails + 1
203     print *, '>>> FAILED test -- FN_VLinearSum Case 2c'
204  else
205     print *, 'PASSED test -- FN_VLinearSum Case 2c'
206  end if
207
208  call FN_VConst(-2.d0, sX)
209  call FN_VConst(1.d0, sY)
210  call FN_VConst(0.d0, sZ)
211  call FN_VLinearSum(1.d0, sX, 1.d0, sY, sZ)
212  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
213     fails = fails + 1
214     print *, '>>> FAILED test -- FN_VLinearSum Case 3'
215  else
216     print *, 'PASSED test -- FN_VLinearSum Case 3'
217  end if
218
219  call FN_VConst(2.d0, sX)
220  call FN_VConst(1.d0, sY)
221  call FN_VConst(0.d0, sZ)
222  call FN_VLinearSum(1.d0, sX, -1.d0, sY, sZ)
223  if (check_ans(1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
224     fails = fails + 1
225     print *, '>>> FAILED test -- FN_VLinearSum Case 4a'
226  else
227     print *, 'PASSED test -- FN_VLinearSum Case 4a'
228  end if
229
230  call FN_VConst(2.d0, sX)
231  call FN_VConst(1.d0, sY)
232  call FN_VConst(0.d0, sZ)
233  call FN_VLinearSum(-1.d0, sX, 1.d0, sY, sZ)
234  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
235     fails = fails + 1
236     print *, '>>> FAILED test -- FN_VLinearSum Case 4b'
237  else
238     print *, 'PASSED test -- FN_VLinearSum Case 4b'
239  end if
240
241  call FN_VConst(2.d0, sX)
242  call FN_VConst(-0.5d0, sY)
243  call FN_VConst(0.d0, sZ)
244  call FN_VLinearSum(1.d0, sX, 2.d0, sY, sZ)
245  if (check_ans(1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
246     fails = fails + 1
247     print *, '>>> FAILED test -- FN_VLinearSum Case 5a'
248  else
249     print *, 'PASSED test -- FN_VLinearSum Case 5a'
250  end if
251
252  call FN_VConst(0.5d0, sX)
253  call FN_VConst(-2.d0, sY)
254  call FN_VConst(0.d0, sZ)
255  call FN_VLinearSum(2.d0, sX, 1.d0, sY, sZ)
256  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
257     fails = fails + 1
258     print *, '>>> FAILED test -- FN_VLinearSum Case 5b'
259  else
260     print *, 'PASSED test -- FN_VLinearSum Case 5b'
261  end if
262
263  call FN_VConst(-2.d0, sX)
264  call FN_VConst(-0.5d0, sY)
265  call FN_VConst(0.d0, sZ)
266  call FN_VLinearSum(-1.d0, sX, 2.d0, sY, sZ)
267  if (check_ans(1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
268     fails = fails + 1
269     print *, '>>> FAILED test -- FN_VLinearSum Case 6a'
270  else
271     print *, 'PASSED test -- FN_VLinearSum Case 6a'
272  end if
273
274  call FN_VConst(0.5d0, sX)
275  call FN_VConst(2.d0, sY)
276  call FN_VConst(0.d0, sZ)
277  call FN_VLinearSum(2.d0, sX, -1.d0, sY, sZ)
278  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
279     fails = fails + 1
280     print *, '>>> FAILED test -- FN_VLinearSum Case 6b'
281  else
282     print *, 'PASSED test -- FN_VLinearSum Case 6b'
283  end if
284
285  call FN_VConst(1.d0, sX)
286  call FN_VConst(-0.5d0, sY)
287  call FN_VConst(0.d0, sZ)
288  call FN_VLinearSum(2.d0, sX, 2.d0, sY, sZ)
289  if (check_ans(1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
290     fails = fails + 1
291     print *, '>>> FAILED test -- FN_VLinearSum Case 7'
292  else
293     print *, 'PASSED test -- FN_VLinearSum Case 7'
294  end if
295
296  call FN_VConst(0.5d0, sX)
297  call FN_VConst(1.d0, sY)
298  call FN_VConst(0.d0, sZ)
299  call FN_VLinearSum(2.d0, sX, -2.d0, sY, sZ)
300  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
301     fails = fails + 1
302     print *, '>>> FAILED test -- FN_VLinearSum Case 8'
303  else
304     print *, 'PASSED test -- FN_VLinearSum Case 8'
305  end if
306
307  call FN_VConst(1.d0, sX)
308  call FN_VConst(-2.d0, sY)
309  call FN_VConst(0.d0, sZ)
310  call FN_VLinearSum(2.d0, sX, 0.5d0, sY, sZ)
311  if (check_ans(1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
312     fails = fails + 1
313     print *, '>>> FAILED test -- FN_VLinearSum Case 9'
314  else
315     print *, 'PASSED test -- FN_VLinearSum Case 9'
316  end if
317
318  ! test FN_VProd
319  call FN_VConst(2.d0, sX)
320  call FN_VConst(-0.5d0, sY)
321  call FN_VConst(0.d0, sZ)
322  call FN_VProd(sX, sY, sZ)
323  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
324     fails = fails + 1
325     print *, '>>> FAILED test -- FN_VProd'
326  else
327     print *, 'PASSED test -- FN_VProd'
328  end if
329
330  ! test FN_VDiv
331  call FN_VConst(1.d0, sX)
332  call FN_VConst(2.d0, sY)
333  call FN_VConst(0.d0, sZ)
334  call FN_VDiv(sX, sY, sZ)
335  if (check_ans(0.5d0, 1.d-14, Nvar, N, sZ) /= 0) then
336     fails = fails + 1
337     print *, '>>> FAILED test -- FN_VDiv'
338  else
339     print *, 'PASSED test -- FN_VDiv'
340  end if
341
342  ! test FN_VScale
343  call FN_VConst(0.5d0, sX)
344  call FN_VScale(2.d0, sX, sX)
345  if (check_ans(1.d0, 1.d-14, Nvar, N, sX) /= 0) then
346     fails = fails + 1
347     print *, '>>> FAILED test -- FN_VScale Case 1'
348  else
349     print *, 'PASSED test -- FN_VScale Case 1'
350  end if
351
352  call FN_VConst(-1.d0, sX)
353  call FN_VConst(0.d0, sZ)
354  call FN_VScale(1.d0, sX, sZ)
355  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
356     fails = fails + 1
357     print *, '>>> FAILED test -- FN_VScale Case 2'
358  else
359     print *, 'PASSED test -- FN_VScale Case 2'
360  end if
361
362  call FN_VConst(-1.d0, sX)
363  call FN_VConst(0.d0, sZ)
364  call FN_VScale(-1.d0, sX, sZ)
365  if (check_ans(1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
366     fails = fails + 1
367     print *, '>>> FAILED test -- FN_VScale Case 3'
368  else
369     print *, 'PASSED test -- FN_VScale Case 3'
370  end if
371
372  call FN_VConst(-0.5d0, sX)
373  call FN_VConst(0.d0, sZ)
374  call FN_VScale(2.d0, sX, sZ)
375  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
376     fails = fails + 1
377     print *, '>>> FAILED test -- FN_VScale Case 4'
378  else
379     print *, 'PASSED test -- FN_VScale Case 4'
380  end if
381
382  ! test FN_VAbs
383  call FN_VConst(-1.d0, sX)
384  call FN_VConst(0.d0, sZ)
385  call FN_VAbs(sX, sZ)
386  if (check_ans(1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
387     fails = fails + 1
388     print *, '>>> FAILED test -- FN_VAbs'
389  else
390     print *, 'PASSED test -- FN_VAbs'
391  end if
392
393  ! test FN_VInv
394  call FN_VConst(2.d0, sX)
395  call FN_VConst(0.d0, sZ)
396  call FN_VInv(sX, sZ)
397  if (check_ans(0.5d0, 1.d-14, Nvar, N, sZ) /= 0) then
398     fails = fails + 1
399     print *, '>>> FAILED test -- FN_VInv'
400  else
401     print *, 'PASSED test -- FN_VInv'
402  end if
403
404  ! test FN_VAddConst
405  call FN_VConst(1.d0, sX)
406  call FN_VConst(0.d0, sZ)
407  call FN_VAddConst(sX, -2.d0, sZ)
408  if (check_ans(-1.d0, 1.d-14, Nvar, N, sZ) /= 0) then
409     fails = fails + 1
410     print *, '>>> FAILED test -- FN_VAddConst'
411  else
412     print *, 'PASSED test -- FN_VAddConst'
413  end if
414
415  ! test FN_VDotProd
416  call FN_VConst(2.d0, sX)
417  call FN_VConst(0.5d0, sY)
418  if (dabs(FN_VDotProd(sX,sY) - (N*Nvar)) > 1.d-14) then
419     fails = fails + 1
420     print *, '>>> FAILED test -- FN_VDotProd (',FN_VDotProd(sX,sY),' /= ',N*Nvar,')'
421  else
422     print *, 'PASSED test -- FN_VDotProd'
423  end if
424
425  ! test FN_VMaxNorm
426  call FN_VConst(-0.5d0, sX)
427  X%data(Nvar,N) = -2.d0
428  if (dabs(FN_VMaxNorm(sX) - 2.d0) > 1.d-14) then
429     fails = fails + 1
430     print *, '>>> FAILED test -- FN_VMaxNorm (',FN_VMaxNorm(sX),' /= 2.d0)'
431  else
432     print *, 'PASSED test -- FN_VMaxNorm'
433  end if
434
435  ! test FN_VWrmsNorm
436  call FN_VConst(-0.5d0, sX)
437  call FN_VConst(0.5d0, sY)
438  if (dabs(FN_VWrmsNorm(sX,sY) - 0.25d0) > 1.d-14) then
439     fails = fails + 1
440     print *, '>>> FAILED test -- FN_VWrmsNorm (',FN_VWrmsNorm(sX,sY),' /= 0.25d0)'
441  else
442     print *, 'PASSED test -- FN_VWrmsNorm'
443  end if
444
445  ! test FN_VWrmsNormMask
446  call FN_VConst(-0.5d0, sX)
447  call FN_VConst(0.5d0, sY)
448  call FN_VConst(1.d0, sZ)
449  Z%data(Nvar,N) = 0.d0
450  fac = dsqrt(1.d0*(N*Nvar - 1)/(N*Nvar))*0.25d0
451  if (dabs(FN_VWrmsNormMask(sX,sY,sZ) - fac) > 1.d-14) then
452     fails = fails + 1
453     print *, '>>> FAILED test -- FN_VWrmsNormMask (',FN_VWrmsNormMask(sX,sY,sZ),' /= ',fac,')'
454  else
455     print *, 'PASSED test -- FN_VWrmsNormMask'
456  end if
457
458  ! test FN_VMin
459  call FN_VConst(2.d0, sX)
460  X%data(Nvar,N) = -2.d0
461  if (dabs(FN_VMin(sX) + 2.d0) > 1.d-14) then
462     fails = fails + 1
463     print *, '>>> FAILED test -- FN_VMin (',FN_VMin(sX),' /= -2.d0)'
464  else
465     print *, 'PASSED test -- FN_VMin'
466  end if
467
468  ! test FN_VWL2Norm
469  call FN_VConst(-0.5d0, sX)
470  call FN_VConst(0.5d0, sY)
471  fac = dsqrt(1.d0*N*Nvar)*0.25d0
472  if (dabs(FN_VWL2Norm(sX,sY) - fac) > 1.d-14) then
473     fails = fails + 1
474     print *, '>>> FAILED test -- FN_VWL2Norm (',FN_VWL2Norm(sX,sY),' /= ',fac,')'
475  else
476     print *, 'PASSED test -- FN_VWL2Norm'
477  end if
478
479  ! test FN_VL1Norm
480  call FN_VConst(-1.d0, sX)
481  fac = 1.d0*N*Nvar
482  if (dabs(FN_VL1Norm(sX) - fac) > 1.d-14) then
483     fails = fails + 1
484     print *, '>>> FAILED test -- FN_VL1Norm (',FN_VL1Norm(sX),' /= ',fac,')'
485  else
486     print *, 'PASSED test -- FN_VL1Norm'
487  end if
488
489  ! test FN_VCompare
490  call FN_VConst(-1.d0, sZ)
491  do j = 1,N
492     do i = 1,Nvar
493        loc = mod((j-1)*Nvar + i - 1, 3)
494        if (loc == 0)  X%data(i,j) = 0.d0
495        if (loc == 1)  X%data(i,j) = -1.d0
496        if (loc == 2)  X%data(i,j) = -2.d0
497     end do
498  end do
499  call FN_VCompare(1.d0, sX, sZ)
500  failure = .false.
501  do j = 1,N
502     do i = 1,Nvar
503        loc = mod((j-1)*Nvar + i - 1, 3)
504        if ((loc == 0) .and. (Z%data(i,j) /= 0.d0))  failure = .true.
505        if ((loc == 1) .and. (Z%data(i,j) /= 1.d0))  failure = .true.
506        if ((loc == 2) .and. (Z%data(i,j) /= 1.d0))  failure = .true.
507     end do
508  end do
509  if (failure) then
510     fails = fails + 1
511     print *, '>>> FAILED test -- FN_VCompare'
512  else
513     print *, 'PASSED test -- FN_VCompare'
514  end if
515
516  ! test FN_VInvTest
517  call FN_VConst(0.5d0, sX)
518  call FN_VConst(0.d0, sZ)
519  failure = (FN_VInvTest(sX, sZ) == 0)
520  if ((check_ans(2.d0, 1.d-14, Nvar, N, sZ) /= 0) .or. failure) then
521     fails = fails + 1
522     print *, '>>> FAILED test -- FN_VInvTest Case 1'
523  else
524     print *, 'PASSED test -- FN_VInvTest Case 1'
525  end if
526
527  failure = .false.
528  call FN_VConst(0.d0, sZ)
529  do j = 1,N
530     do i = 1,Nvar
531        loc = mod((j-1)*Nvar + i - 1, 2)
532        if (loc == 0)  X%data(i,j) = 0.d0
533        if (loc == 1)  X%data(i,j) = 0.5d0
534     end do
535  end do
536  if (FN_VInvTest(sX, sZ) == 1)  failure = .true.
537  do j = 1,N
538     do i = 1,Nvar
539        loc = mod((j-1)*Nvar + i - 1, 2)
540        if ((loc == 0) .and. (Z%data(i,j) /= 0.d0))  failure = .true.
541        if ((loc == 1) .and. (Z%data(i,j) /= 2.d0))  failure = .true.
542     end do
543  end do
544  if (failure) then
545     fails = fails + 1
546     print *, '>>> FAILED test -- FN_VInvTest Case 2'
547  else
548     print *, 'PASSED test -- FN_VInvTest Case 2'
549  end if
550
551  ! test FN_VConstrMask
552  call FN_VConst(-1.d0, sZ)
553  do j = 1,N
554     do i = 1,Nvar
555        loc = mod((j-1)*Nvar + i - 1, 7)
556        if (loc == 0) then  ! y = -2, test for < 0
557           Y%data(i,j) = -2.d0
558           X%data(i,j) = -2.d0
559        end if
560        if (loc == 1) then ! y = -1, test for <= 0
561           Y%data(i,j) = -1.d0
562           X%data(i,j) = -1.d0
563        end if
564        if (loc == 2) then ! y = -1, test for == 0
565           Y%data(i,j) = -1.d0
566           X%data(i,j) = 0.d0
567        end if
568        if (loc == 3) then ! y = 0, no test
569           Y%data(i,j) = 0.d0
570           X%data(i,j) = 0.5d0
571        end if
572        if (loc == 4) then ! y = 1, test for == 0
573           Y%data(i,j) = 1.d0
574           X%data(i,j) = 0.d0
575        end if
576        if (loc == 5) then ! y = 1, test for >= 0
577           Y%data(i,j) = 1.d0
578           X%data(i,j) = 1.d0
579        end if
580        if (loc == 6) then ! y = 2, test for > 0
581           Y%data(i,j) = 2.d0
582           X%data(i,j) = 2.d0
583        end if
584     end do
585  end do
586  failure = .false.
587  if (FN_VConstrMask(sY, sX, sZ) /= 1) then
588     failure = .true.
589  end if
590  if ((check_ans(0.d0, 1.d-14, Nvar, N, sZ) /= 0) .or. failure) then
591     fails = fails + 1
592     print *, '>>> FAILED test -- FN_VConstrMask Case 1'
593  else
594     print *, 'PASSED test -- FN_VConstrMask Case 1'
595  end if
596
597  call FN_VConst(-1.d0, sZ)
598  do j = 1,N
599     do i = 1,Nvar
600        loc = mod((j-1)*Nvar + i - 1, 5)
601        if (loc == 0) then  ! y = -2, test for < 0
602           Y%data(i,j) = -2.d0
603           X%data(i,j) = 2.d0
604        end if
605        if (loc == 1) then ! y = -1, test for <= 0
606           Y%data(i,j) = -1.d0
607           X%data(i,j) = 1.d0
608        end if
609        if (loc == 2) then ! y = 0, no test
610           Y%data(i,j) = 0.d0
611           X%data(i,j) = 0.5d0
612        end if
613        if (loc == 3) then ! y = 1, test for >= 0
614           Y%data(i,j) = 1.d0
615           X%data(i,j) = -1.d0
616        end if
617        if (loc == 4) then ! y = 2, test for > 0
618           Y%data(i,j) = 2.d0
619           X%data(i,j) = -2.d0
620        end if
621     end do
622  end do
623  failure = .false.
624  if (FN_VConstrMask(sY, sX, sZ) /= 0) then
625     failure = .true.
626  end if
627  do j = 1,N
628     do i = 1,Nvar
629        loc = mod((j-1)*Nvar + i - 1, 5)
630        if (loc == 2) then
631           if (Z%data(i,j) /= 0.d0) failure = .true.
632        else
633           if (Z%data(i,j) /= 1.d0) failure = .true.
634        end if
635     end do
636  end do
637  if (failure) then
638     fails = fails + 1
639     print *, '>>> FAILED test -- FN_VConstrMask Case 2'
640  else
641     print *, 'PASSED test -- FN_VConstrMask Case 2'
642  end if
643
644  ! test FN_VMinQuotient
645  call FN_VConst(2.d0, sX)
646  call FN_VConst(2.d0, sY)
647  X%data(Nvar,N) = 0.5d0
648  fac = 0.25d0
649  if (dabs(FN_VMinQuotient(sX,sY) - fac) > 1.d-14) then
650     fails = fails + 1
651     print *, '>>> FAILED test -- FN_VMinQuotient Case 1'
652  else
653     print *, 'PASSED test -- FN_VMinQuotient Case 1'
654  end if
655
656  call FN_VConst(2.d0, sX)
657  call FN_VConst(0.d0, sY)
658  fac = 1.d307
659  if (dabs(FN_VMinQuotient(sX,sY) - fac) > 1.d-14) then
660     fails = fails + 1
661     print *, '>>> FAILED test -- FN_VMinQuotient Case 2'
662  else
663     print *, 'PASSED test -- FN_VMinQuotient Case 2'
664  end if
665
666  ! test FN_VWSqrSumLocal
667  call FN_VConst(-1.d0, sX)
668  call FN_VConst(0.5d0, sY)
669  fac = 0.25d0*N*Nvar
670  if (dabs(FN_VWSqrSumLocal(sX,sY) - fac) > 1.d-14) then
671     fails = fails + 1
672     print *, '>>> FAILED test -- FN_VWSqrSumLocal (',FN_VWSqrSumLocal(sX,sY),' /= ',fac,')'
673  else
674     print *, 'PASSED test -- FN_VWSqrSumLocal'
675  end if
676
677
678  ! test FN_VWSqrSumMaskLocal
679  call FN_VConst(-1.d0, sX)
680  call FN_VConst(0.5d0, sY)
681  call FN_VConst(1.d0, sZ)
682  Z%data(Nvar,N) = 0.d0
683  fac = 0.25d0*(N*Nvar-1)
684  if (dabs(FN_VWSqrSumMaskLocal(sX,sY,sZ) - fac) > 1.d-14) then
685     fails = fails + 1
686     print *, '>>> FAILED test -- FN_VWSqrSumMaskLocal (',FN_VWSqrSumMaskLocal(sX,sY,sZ),' /= ',fac,')'
687  else
688     print *, 'PASSED test -- FN_VWSqrSumMaskLocal'
689  end if
690
691  ! free vectors
692  call FN_VDestroy(sU)
693  call FN_VDestroy(sV)
694  call FN_VDestroy(sW)
695  call FN_VDestroy(sX)
696  call FN_VDestroy(sY)
697  call FN_VDestroy(sZ)
698
699  ! print results
700  if (fails > 0) then
701     print '(a,i3,a)', 'FAIL: FNVector module failed ',fails,' tests'
702     stop 1
703  else
704     print *, 'SUCCESS: FNVector module passed all tests'
705  end if
706  print *, '  '
707
708end program main
709