11module test_fftpack_dct
22
3- use fftpack, only: rk, dcosti, dcost, dct, idct
3+ use fftpack, only: rk, dcosti, dcost, dct, idct, dcosqi, dcosqf, dcosqb
44 use testdrive, only: new_unittest, unittest_type, error_type, check
55 implicit none
66 private
@@ -32,21 +32,54 @@ subroutine test_classic_dct(error)
3232 call dcost(4 , x, w)
3333 call check(error, all (x == [real (kind= rk) :: 15 , - 4 , 0 , - 1.0000000000000009_rk ]), " `dcosti` failed." )
3434 if (allocated (error)) return
35-
3635 call dcost(4 , x, w)
37- call check(error, all (x/ (2.0_rk * (4.0_rk - 1.0_rk )) == [real (kind= rk) :: 1 , 2 , 3 , 4 ]), " `dcost` failed." )
36+ call check(error, all (x/ (2 * 3 ) == [real (kind= rk) :: 1 , 2 , 3 , 4 ]), " `dcost` failed." )
37+
38+ x = [1 , 2 , 3 , 4 ]
39+ call dcosqi(4 , w)
40+ call dcosqf(4 , x, w)
41+ call check(error, sum (abs (x - [11.999626276085150_rk , - 9.1029432177492193_rk , &
42+ 2.6176618435106480_rk , - 1.5143449018465791_rk ])) < eps, &
43+ " `dcosqf` failed." )
44+ if (allocated (error)) return
45+ call dcosqb(4 , x, w)
46+ call check(error, sum (abs (x/ (4 * 4 ) - [real (kind= rk) :: 1 , 2 , 3 , 4 ])) < eps, &
47+ " `dcosqb` failed." )
3848
3949 end subroutine test_classic_dct
4050
4151 subroutine test_modernized_dct (error )
4252 type (error_type), allocatable , intent (out ) :: error
53+ real (kind= rk) :: eps = 1.0e-10_rk
4354 real (kind= rk) :: x(3 ) = [9 , - 9 , 3 ]
4455
45- call check(error, all (dct(x, 2 ) == [real (kind= rk) :: 0 , 18 ]), " `dct(x, 2)` failed." )
56+ ! DCT-1
57+ call check(error, sum (abs (dct(x,2 ,1 ) - [0.0_rk , 18.0_rk ])) < eps, " `dct(x,2,1)` failed." )
58+ if (allocated (error)) return
59+ call check(error, sum (abs (dct(x,3 ,1 ) - dct(x,t= 1 ))) < eps, " `dct(x,3,1)` failed." )
60+ if (allocated (error)) return
61+ call check(error, sum (abs (dct(x,4 ,1 ) - [real (kind= rk) :: - 3 , - 3.0000000000000036_rk , 15 , 33 ])) < eps,&
62+ " `dct(x,4,1)` failed." )
63+ ! DCT-2
64+ x = [9 , - 9 , 3 ]
65+ call check(error, sum (abs (dct(x,3 ,2 ) - [12.0_rk , 20.784609690826525_rk , 60.0_rk ])) < eps,&
66+ " `dct(x,3,2)` failed." )
67+ call check(error, sum (abs (dct(x,3 ) - [12.0_rk , 20.784609690826525_rk , 60.0_rk ])) < eps,&
68+ " `dct(x,3)` failed." )
69+ call check(error, sum (abs (dct(x,4 ,2 ) - [12.0_rk , 14.890858416882008_rk , 42.426406871192853_rk ,&
70+ 58.122821125684993_rk ])) < eps, " `dct(x,4,2)` failed." )
71+
72+ ! DCT-3
73+ x = [9 , - 9 , 3 ]
74+ call check(error, sum (abs (dct(x,2 ,3 ) - [- 3.7279220613578570_rk , 21.727922061357859_rk ])) < eps, &
75+ " `dct(x,2,3)` failed." )
4676 if (allocated (error)) return
47- call check(error, all (dct(x, 3 ) == dct(x)), " `dct(x, 3)` failed." )
77+ call check(error, sum (abs (dct(x,3 ,3 ) - dct(x,t= 3 ))) < eps, &
78+ " `dct(x,3,3)` failed." )
4879 if (allocated (error)) return
49- call check(error, all (dct(x, 4 ) == [real (kind= rk) :: - 3 , - 3.0000000000000036_rk , 15 , 33 ]), " `dct(x, 4)` failed." )
80+ call check(error, sum (abs (dct(x,4 ,3 ) - [- 3.3871908980838743_rk , - 2.1309424696909023_rk , &
81+ 11.645661095452331_rk , 29.872472272322447_rk ])) < eps, &
82+ " `dct(x,n=4,t=3)` failed." )
5083
5184 end subroutine test_modernized_dct
5285
@@ -55,15 +88,37 @@ subroutine test_modernized_idct(error)
5588 real (kind= rk) :: eps = 1.0e-10_rk
5689 real (kind= rk) :: x(4 ) = [1 , 2 , 3 , 4 ]
5790
58- call check(error, all (idct(dct(x))/ (2.0_rk * (4.0_rk - 1.0_rk )) == [real (kind= rk) :: 1 , 2 , 3 , 4 ]), &
59- " `idct(dct(x))/(2.0_rk*(4.0_rk-1.0_rk))` failed." )
91+ ! inverse DCT-1
92+ call check(error, sum (abs (idct(dct(x,t= 1 ),t= 1 )/ (2 * 3 ) - x)) < eps, " `idct(dct(x,t=1),t=1)/(2*3)` failed." )
93+ if (allocated (error)) return
94+ call check(error, sum (abs (idct(dct(x,t= 1 ), 2 , 1 )/ (2 * 1 ) - [5.5_rk , 9.5_rk ])) < eps,&
95+ " `idct(dct(x,t=1), 2, 1)/(2*1)` failed." )
96+ if (allocated (error)) return
97+ call check(error, sum (abs (idct(dct(x,2 ,1 ), 4 , 1 )/ (2 * 3 ) - [0.16666666666666666_rk , 0.33333333333333331_rk ,&
98+ 0.66666666666666663_rk , 0.83333333333333315_rk ])) < eps,&
99+ " `idct(dct(x,2,1), 4, 1)/(2*3)` failed." )
100+
101+ ! inverse DCT-2
102+ x = [1 , 2 , 3 , 4 ]
103+ call check(error, sum (abs (idct(dct(x,t= 2 ))/ (4 * 4 ) - x)) < eps, &
104+ " `idct(dct(x, t=2))/(4*4)` failed." )
105+ if (allocated (error)) return
106+ call check(error, sum (abs (idct(dct(x,t= 2 ),n= 2 ) - [22.156460020898692_rk , 57.843539979101308_rk ])) < eps,&
107+ " `idct(dct(x,t=2),n=2)` failed." )
108+ if (allocated (error)) return
109+ call check(error, sum (abs (idct(dct(x,n= 2 ,t= 2 ),n= 4 ) - [6.7737481404944937_rk , 9.8352155994152106_rk ,&
110+ 14.164784400584789_rk , 17.226251859505506_rk ])) < eps, " `idct(dct(x,n=2,t=2),n=4)` failed." )
111+
112+ ! inverse DCT-3
113+ x = [1 , 2 , 3 , 4 ]
114+ call check(error, sum (abs (idct(dct(x,t= 3 ),t= 3 )/ (4 * 4 ) - x)) < eps, &
115+ " `idct(dct(x, t=3), t=3)/(4*4)` failed." )
60116 if (allocated (error)) return
61- call check(error, all ( idct(dct(x), 2 )/ (2.0_rk * ( 2.0_rk - 1.0_rk )) == [ real (kind = rk) :: 5.5 , 9.5 ]) , &
62- " `idct(dct(x), 2 )/(2.0_rk*(2.0_rk-1.0_rk) )` failed." )
117+ call check(error, sum ( abs ( idct(dct(x,t = 3 ),n = 2 ,t = 3 )/ (4 * 2 ) - [ 1.4483415291679655_rk , 7.4608849947753271_rk ])) < eps , &
118+ " `idct(dct(x,t=3),n=2,t=3 )/(4*2 )` failed." )
63119 if (allocated (error)) return
64- call check(error, all (idct(dct(x, 2 ), 4 )/ (2.0_rk * (4.0_rk - 1.0_rk )) == &
65- [0.16666666666666666_rk , 0.33333333333333331_rk , 0.66666666666666663_rk , 0.83333333333333315_rk ]), &
66- " `idct(dct(x, 2), 4)/(2.0_rk*(4.0_rk-1.0_rk))` failed." )
120+ call check(error, sum (abs (idct(dct(x,n= 2 ,t= 3 ),n= 4 ,t= 3 )/ (4 * 4 ) - [0.5_rk , 0.70932417358418376_rk , 1.0_rk , &
121+ 0.78858050747473762_rk ])) < eps, " `idct(dct(x,n=2,t=3),n=4, t=3)/(4*4)` failed." )
67122
68123 end subroutine test_modernized_idct
69124
0 commit comments