1 /*
2 * $Id: testm.i,v 1.1 2005-09-18 22:06:12 dhmunro Exp $
3 * Certify functions in fft.i and matrix.i
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 func testm
12 {
13 elapsed= old_elapsed= [0., 0., 0.];
14 write, "testing fft routines...";
15 timer, old_elapsed;
16 force2d = 0;
17 fft_test, 95;
18 fft_test, 32;
19 fft_test, 12;
20 fft_test, 1024;
21 for (i=0 ; i<32 ; i++) fft_test, 4096;
22 force2d = 1; fft_test, 128; force2d = 0;
23 timer, elapsed;
24 timer_print, "FFT elapsed time", elapsed-old_elapsed;
25 write, "testing LU decomposition routine...";
26 timer, old_elapsed;
27 testLU, 75;
28 testLU, 16;
29 testLU, 216;
30 timer, elapsed;
31 timer_print, "LU elapsed time", elapsed-old_elapsed;
32 write, "testing QR decomposition routine...";
33 timer, old_elapsed;
34 testQR, 75, 75;
35 testQR, 16, 16;
36 testQR, 25, 21;
37 testQR, 21, 25;
38 testQR, 256, 240;
39 timer, elapsed;
40 timer_print, "QR elapsed time", elapsed-old_elapsed;
41 write, "testing SVD routine...";
42 timer, old_elapsed;
43 testSVD, 75, 75;
44 testSVD, 16, 16;
45 testSVD, 25, 21;
46 testSVD, 21, 25;
47 testSV, 21;
48 testSV, 64;
49 testSV, 128;
50 timer, elapsed;
51 timer_print, "SVD elapsed time", elapsed-old_elapsed;
52 }
53
fft_test(n)54 func fft_test(n)
55 {
56 index= 2*pi*(indgen(n)-1.0)/n;
57 z= sin(index*3);
58 zf= fft(z, 1);
59 z3= z2= array(0i, n);
60 z3(4)= -0.5i*n;
61 z3(-2)= 0.5i*n;
62 zb= fft(z, -1);
63 if (max(abs(zf-z3))>1.e-12*n || max(abs(zb-conj(z3)))>1.e-12*n)
64 write, "***WARNING*** failed 1D fft test";
65 if (n<=96 || force2d) {
66 z*= cos(index*2)(-,);
67 zf= fft(z, [0, 1]);
68 z2(3)= z2(-1)= 0.5*n;
69 zb= fft(z, [], [-1, 0]);
70 if (max(abs(zf-sin(index*3)*z2(-,)))>1.e-12*n ||
71 max(abs(zb-conj(z3)*cos(index*2)(-,)))>1.e-12*n)
72 write, "***WARNING*** failed first 2D fft test";
73 zf= fft(z, 1);
74 zb= fft(z, -1);
75 if (max(abs(zf-z3*z2(-,)))>1.e-12*n ||
76 max(abs(zb-conj(z3)*z2(-,)))>1.e-12*n)
77 write, "***WARNING*** failed second 2D fft test";
78 }
79 }
80
TDcheck(c,d,e,b,x,s)81 func TDcheck(c, d, e, b, x, s)
82 {
83 check= _( d(1)*x(1) + e(1)*x(2),
84 c(1:-1)*x(1:-2) + d(2:-1)*x(2:-1) + e(2:0)*x(3:0),
85 c(0)*x(-1) + d(0)*x(0) );
86 if (max(abs(check-b))>1.e-9*max(abs(b))) {
87 write, "***WARNING*** "+s+" tridiagonal solution doesn't check";
88 write, " max relative error is "+pr1((max(abs(check-b)))/max(abs(b)));
89 }
90 }
91
testTD(n)92 func testTD(n)
93 {
94 c= random(n-1);
95 d= random(n);
96 e= random(n-1);
97 b= random(n);
98 TDcheck,c,d,e,b,TDsolve(c,d,e,b), "1D";
99 b2= random(n);
100 x= TDsolve(c,d,e,[b,b2])
101 TDcheck,c,d,e,b, x(,1), "2D(1)";
102 TDcheck,c,d,e,b2, x(,2), "2D(2)";
103 x= TDsolve(c,d,e,transpose([b,b2]), which=2)
104 TDcheck,c,d,e,b, x(1,), "2D(1)/which";
105 TDcheck,c,d,e,b2, x(2,), "2D(2)/which";
106 }
107
LUcheck(a,b,x,s)108 func LUcheck(a, b, x, s)
109 {
110 check= a(,+)*x(+);
111 if (max(abs(check-b))>1.e-9*max(abs(b))) {
112 write, "***WARNING*** "+s+" LUsolve solution doesn't check";
113 write, " max relative error is "+pr1((max(abs(check-b)))/max(abs(b)));
114 }
115 }
116
testLU(n)117 func testLU(n)
118 {
119 a= random(n,n);
120 b= random(n);
121 LUcheck,a,b,LUsolve(a,b), "1D";
122 b2= random(n);
123 x= LUsolve(a,[b,b2])
124 LUcheck,a,b, x(,1), "2D(1)";
125 LUcheck,a,b2, x(,2), "2D(2)";
126 x= LUsolve(transpose(a),[b,b2])
127 LUcheck,transpose(a),b, x(,1), "t2D(1)";
128 LUcheck,transpose(a),b2, x(,2), "t2D(2)";
129 x= LUsolve(a,transpose([b,b2]), which=2)
130 LUcheck,a,b, x(1,), "2D(1)/which";
131 LUcheck,a,b2, x(2,), "2D(2)/which";
132 x= LUsolve(transpose(a),transpose([b,b2]), which=2)
133 LUcheck,transpose(a),b, x(1,), "t2D(1)/which";
134 LUcheck,transpose(a),b2, x(2,), "t2D(2)/which";
135 ai= LUsolve(a);
136 err= max(abs(ai(,+)*a(+,)-unit(n)))/max(abs(ai));
137 if (err>1.e-9) {
138 write, "***WARNING*** LUsolve inverse is fishy";
139 write, " max relative error is "+pr1(err);
140 }
141 }
142
QRcheck(a,b,x,s)143 func QRcheck(a, b, x, s)
144 {
145 check= a(,+)*x(+);
146 err= max(abs(check-b))/max(abs(b));
147 if (err>1.e-9 && dimsof(a)(2)<=dimsof(a)(3)) {
148 write, "***WARNING*** "+s+" QRsolve solution doesn't check";
149 write, " max relative error is "+pr1(err);
150 }
151 }
152
testQR(m,n)153 func testQR(m,n)
154 {
155 a= random(m,n);
156 b= random(m);
157 QRcheck,a,b,QRsolve(a,b), "1D";
158 b2= random(m);
159 x= QRsolve(a,[b,b2])
160 QRcheck,a,b, x(,1), "2D(1)";
161 QRcheck,a,b2, x(,2), "2D(2)";
162 x= QRsolve(a,transpose([b,b2]), which=2)
163 QRcheck,a,b, x(1,), "2D(1)/which";
164 QRcheck,a,b2, x(2,), "2D(2)/which";
165 }
166
SVcheck(a,b,x,s)167 func SVcheck(a, b, x, s)
168 {
169 check= a(,+)*x(+);
170 err= max(abs(check-b))/max(abs(b));
171 if (err>1.e-9) {
172 write, "***WARNING*** "+s+" SVsolve solution doesn't check";
173 write, " max relative error is "+pr1(err);
174 }
175 }
176
testSVD(m,n)177 func testSVD(m, n)
178 {
179 a= random(m,n);
180 s= SVdec(a, u, v);
181 if (anyof(s(dif)>0.0))
182 error, "***WARNING*** SVdec returned increasing singular values";
183 achk= u(,+) * (s*v)(+,);
184 sabs= max(abs(s));
185 err= max(abs(a-achk))/sabs;
186 if (err>1.e-9) {
187 write, "***WARNING*** SVdec decomposition doesn't check";
188 write, " max relative error is "+pr1(err);
189 }
190
191 s= SVdec(a, u, v, full=1);
192 if (anyof(s(dif)>0.0))
193 error, "***WARNING*** SVdec returned increasing singular values";
194 uu= u(,1:min(m,n));
195 vv= v(1:min(m,n),);
196 achk= uu(,+) * (s*vv)(+,);
197 sabs= max(abs(s));
198 err= max(abs(a-achk))/sabs;
199 if (err>1.e-9) {
200 write, "***WARNING*** SVdec decomposition doesn't check";
201 write, " max relative error is "+pr1(err);
202 }
203 err= max(abs(u(,+)*u(,+)-unit(m)))+max(abs(v(,+)*v(,+)-unit(n)));
204 if (err>1.e-9) {
205 write, "***WARNING*** SVdec decomposition not orthogonal";
206 write, " max relative error is "+pr1(err);
207 }
208 }
209
testSV(n)210 func testSV(n)
211 {
212 a= random(n,n);
213 b= random(n);
214 SVcheck,a,b,SVsolve(a,b), "1D";
215 b2= random(n);
216 x= SVsolve(a,[b,b2])
217 SVcheck,a,b, x(,1), "2D(1)";
218 SVcheck,a,b2, x(,2), "2D(2)";
219 x= SVsolve(a,transpose([b,b2]), which=2)
220 SVcheck,a,b, x(1,), "2D(1)/which";
221 SVcheck,a,b2, x(2,), "2D(2)/which";
222 }
223