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