linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / tests / unit_tests / schur.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) ????-2008 - INRIA Michael Baudin
4 //
5 //  This file is distributed under the same license as the Scilab package.
6 // =============================================================================
7 function r=Err(x)
8         r=norm(x,1)
9 endfunction
10 rand('normal')
11 ilib_verbose(0);
12
13 libName = [];
14 if getos() == "Windows" then
15     libName = SCI + "/bin/elem_func";
16 end
17
18 //define tools
19 function A=testmat1(a,n)
20         //eigen values are given by a dilation of nth roots of 1
21         A=diag(a*ones(1,n-1),1)+diag((1/a)*ones(1,n-1),-1)
22         A(1,n)=1/a;A(n,1)=a
23 endfunction
24
25 //==========================================================================
26 //==============================    schur     ============================== 
27 //==========================================================================
28 clear sel
29 function t=sel(R),t=real(R)<0 ,endfunction
30 //Empty matrix
31 A=[];
32 if schur(A)<>[] then pause,end
33 if schur(A,'real')<>[] then pause,end
34 if schur(A,'complex')<>[] then pause,end
35
36 if schur(A,'c')<>[] then pause,end
37 if schur(A,'d')<>[] then pause,end
38 if schur(A,sel)<>[] then pause,end
39
40 [U,S]=schur(A);
41 if U<>[]|S<>[] then pause,end
42 [U,S]=schur(A,'real');
43 if U<>[]|S<>[] then pause,end
44 [U,S]=schur(A,'complex');
45 if U<>[]|S<>[] then pause,end
46
47
48
49 [U,N]=schur(A,'c');
50 if U<>[]|N<>0 then pause,end
51 [U,N]=schur(A,'d');
52 if U<>[]|N<>0 then pause,end
53 [U,N]=schur(A,sel);
54 if U<>[]|N<>0 then pause,end
55
56 [U,N,S]=schur(A,'c');
57 if U<>[]|N<>0|S<>[] then pause,end
58 [U,N,S]=schur(A,'d');
59 if U<>[]|N<>0|S<>[] then pause,end
60 [U,N,S]=schur(A,sel);
61 if U<>[]|N<>0|S<>[] then pause,end
62
63 //Rectangular matrix
64 if execstr('schur(rand(2,3))','errcatch')==0 then pause,end
65 if execstr('[U,S]=schur(rand(2,3))','errcatch')==0 then pause,end
66
67 if execstr('schur(rand(2,3)+%i*eye())','errcatch')==0 then pause,end
68 if execstr('[U,S]=schur(rand(2,3)+%i*eye())','errcatch')==0 then pause,end
69
70 //Small dimension
71 A=testmat1(3,5);Ac=testmat1(3+%i,5);
72 //Real
73 [U,S]=schur(A);
74 if Err(triu(S,-1)-S)>%eps then pause,end
75 if Err(U*S*U'-A)>200*%eps then pause,end
76 if Err(schur(A)-S) >%eps then pause,end
77
78 [U,S]=schur(A,'real');
79 if Err(triu(S,-1)-S)>%eps then pause,end
80 if Err(U*S*U'-A)>200*%eps then pause,end
81 if Err(schur(A)-S) >%eps then pause,end
82
83 [U,S]=schur(A,'complex');
84 if Err(triu(S)-S)>%eps then pause,end
85 if Err(U*S*U'-A)>200*%eps then pause,end
86 if Err(schur(A,'complex')-S) >%eps then pause,end
87
88 [U,n]=schur(A,'c');S=U'*A*U;
89 if n<>2 then pause,end
90 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
91 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
92
93 [U,n]=schur(A,'d');S=U'*A*U;
94 if n<>0 then pause,end
95 if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
96
97 [U,n]=schur(A,sel);S=U'*A*U;
98 if n<>2 then pause,end
99 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
100 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
101
102
103 //Complex
104 [U,S]=schur(Ac);
105 if Err(triu(S,-1)-S)>%eps then pause,end
106 if Err(U*S*U'-Ac)>200*%eps then pause,end
107 if Err(schur(Ac)-S) >%eps then pause,end
108
109 [U,S]=schur(Ac,'complex');
110 if Err(triu(S,-1)-S)>%eps then pause,end
111 if Err(U*S*U'-Ac)>200*%eps then pause,end
112 if Err(schur(Ac)-S) >%eps then pause,end
113
114
115 [U,n]=schur(Ac,'c');S=U'*Ac*U;
116 if n<>3 then pause,end
117 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
118 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
119
120 [U,n]=schur(Ac,'d');S=U'*A*U;
121 if n<>0 then pause,end
122 if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
123
124 [U,n]=schur(Ac,sel);S=U'*Ac*U;
125 if n<>3 then pause,end
126 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
127 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
128
129
130 //Large dimension
131 A=testmat1(3,50);Ac=testmat1(3+%i,50);
132 //Real
133 [U,S]=schur(A);
134 if Err(triu(S,-1)-S)>%eps then pause,end
135 if Err(U*S*U'-A)>1000*%eps then pause,end
136 if Err(schur(A)-S) >%eps then pause,end
137
138 [U,S]=schur(A,'real');
139 if Err(triu(S,-1)-S)>%eps then pause,end
140 if Err(U*S*U'-A)>1000*%eps then pause,end
141 if Err(schur(A)-S) >%eps then pause,end
142
143 [U,S]=schur(A,'complex');
144 if Err(triu(S)-S)>%eps then pause,end
145 if Err(U*S*U'-A)>1000*%eps then pause,end
146 if Err(schur(A,'complex')-S) >%eps then pause,end
147
148
149 [U,n]=schur(A,'c');S=U'*A*U;
150 if n<>25 then pause,end
151 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
152 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
153
154 [U,n]=schur(A,'d');S=U'*A*U;
155 if n<>0 then pause,end
156 if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
157
158 [U,n]=schur(A,sel);S=U'*A*U;
159 if n<>25 then pause,end
160 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
161 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
162
163 //Complex
164 [U,S]=schur(Ac);
165 if Err(triu(S,-1)-S)>%eps then pause,end
166 if Err(U*S*U'-Ac)>1000*%eps then pause,end
167 if Err(schur(Ac)-S) >%eps then pause,end
168
169 [U,S]=schur(Ac,'complex');
170 if Err(triu(S,-1)-S)>%eps then pause,end
171 if Err(U*S*U'-Ac)>1000*%eps then pause,end
172 if Err(schur(Ac)-S) >%eps then pause,end
173
174 [U,n]=schur(Ac,'c');S=U'*Ac*U;
175 if n<>25 then pause,end
176 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
177 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
178
179 [U,n]=schur(Ac,'d');S=U'*Ac*U;
180 if n<>0 then pause,end
181 if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
182
183 [U,n]=schur(Ac,sel);S=U'*Ac*U;
184 if n<>25 then pause,end
185 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
186 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
187
188 // Lib part
189 cd TMPDIR;
190 // equal to schur(A, 'c');
191 C=[ 'int mytest1(double* _real, double* _img)'
192     '{' 
193     '    return *_real < 0;'
194     '}'];
195
196 mputl(C,TMPDIR+'/mytest.c');
197 ulink();
198 lp=ilib_for_link('mytest1','mytest.c',[],'c');
199 exec loader.sce;
200 [U,n]=schur(A,'mytest1');S=U'*A*U;
201 if n<>25 then pause,end
202 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
203 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
204
205 // equal to schur(A, 'd');
206 C=[ 'extern double dpythags(double,double);' // YaSp function
207     ''
208     'int mytest2(double* _real, double* _img)'
209     '{' 
210     '    return dpythags(*_real, *_img) < 1;'
211     '}'];
212
213 mputl(C,TMPDIR+'/mytest.c');
214 ulink();
215 lp=ilib_for_link('mytest2','mytest.c', libName,'c');
216
217 exec loader.sce;
218 [U,n]=schur(A,'mytest2');S=U'*A*U;
219 if n<>0 then pause,end
220 if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
221
222 // equal to schur(Ac, 'c');
223 C=[ '#include ""doublecomplex.h""'
224     ''
225     'int mytest3(doublecomplex* _complex)'
226     '{'
227     '    return _complex->r < 0 ? 1 : 0;'
228     '}'];
229 mputl(C,TMPDIR+'/mytest.c');
230 ulink();
231 lp=ilib_for_link('mytest3','mytest.c',[],'c');
232 exec loader.sce;
233 [U,n]=schur(Ac, 'mytest3');S=U'*Ac*U;
234 if n<>25 then pause,end
235 if or(real(spec(S(1:n,1:n)))>=0) then pause,end
236 if or(real(spec(S(n+1:$,n+1:$)))<0) then pause,end
237
238 // equal to schur(Ac, 'd');
239 C=[ '#include ""doublecomplex.h""'
240     ''
241     'extern double dpythags(double,double);' // YaSp function
242     ''
243     'int mytest4(doublecomplex* _complex)'
244     '{'
245     '    if(dpythags(_complex->r, _complex->i) < 1)'
246     '    {'
247     '        return 1;'
248     '    }'
249     '    else'
250     '    {'
251     '        return 0;'
252     '    }'
253     '}'];
254
255 mputl(C,TMPDIR+'/mytest.c');
256 ulink();
257 lp=ilib_for_link('mytest4','mytest.c', libName,'c');
258 exec loader.sce;
259 [U,n]=schur(Ac,'mytest4');S=U'*Ac*U;
260 if n<>0 then pause,end
261 if or(abs(spec(S(n+1:$,n+1:$)))<1) then pause,end
262
263
264 //==========================================================================
265 //==============================    schur part II   ======================== 
266 //==========================================================================
267 //Empty matrix
268 [As,Es]=schur([],[]);
269 if As<>[]|Es<>[] then pause,end
270
271 [As,dim]=schur([],[],'c');
272 if As<>[]|dim<>0 then pause,end
273 [As,dim]=schur([],[],'d');
274 if As<>[]|dim<>0 then pause,end
275 [As,dim]=schur([],[],sel);
276 if As<>[]|dim<>0 then pause,end
277
278 [As,Es,Q,Z]=schur([],[]);
279 if As<>[]|Es<>[]|Q<>[]|Z<>[] then pause,end
280
281 [As,Es,dim]=schur([],[],'c');
282 if As<>[]|Es<>[]|dim<>0 then pause,end
283 [As,Es,dim]=schur([],[],'d');
284 if As<>[]|Es<>[]|dim<>0 then pause,end
285 [As,Es,dim]=schur([],[],sel);
286 if As<>[]|Es<>[]|dim<>0 then pause,end
287
288 [Z,dim]=schur([],[],'c');
289 if Z<>[]|dim<>0 then pause,end
290 [Z,dim]=schur([],[],'d');
291 if Z<>[]|dim<>0 then pause,end
292 [Z,dim]=schur([],[],sel);
293 if Z<>[]|dim<>0 then pause,end
294
295
296 //Rectangular matrix
297 if execstr('[As,Es]=schur(rand(2,3),rand(2,3))','errcatch')==0 then  pause,end
298 if execstr('[As,Es,Q,Z]=schur(rand(2,3),rand(2,3))','errcatch')==0 then  pause,end
299 if execstr('[As,Es,dim]=schur(rand(2,3),rand(2,3),''c'')','errcatch')==0 then  pause,end
300 if execstr('[Z,dim]=schur(rand(2,3),rand(2,3),sel)','errcatch')==0 then  pause,end
301
302 //Small dimension
303 //----Real------------
304 A=testmat1(1,5);E=testmat1(-2,5) ;
305 [As,Es,Q,Z]=schur(A,E);
306 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
307 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
308 if Err(As-Q'*A*Z) >200*%eps then pause,end
309 if Err(Es-Q'*E*Z) >200*%eps then pause,end
310
311 [As1,Es1]=schur(A,E);
312 if Err(As1-As)>10*%eps then pause,end
313 if Err(Es1-Es)>10*%eps then pause,end
314
315 // Ordered 'c'
316 dim=schur(A,E,'c');
317 if dim<>5 then pause,end
318 [Z,dim]=schur(A,E,'c');
319 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
320
321 [Q,Z1,dim]=schur(A,E,'c');
322 if Err(Z1-Z)>10*%eps then pause,end
323 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
324 if dim<>5 then pause,end
325
326 [As,Es,Z,dim]=schur(A,E,'c');
327 if dim<>5 then pause,end
328 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
329 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
330 if Err(As-Q'*A*Z) >200*%eps then pause,end
331 if Err(Es-Q'*E*Z) >200*%eps then pause,end
332 // Ordered 'd'
333 dim=schur(A,E,'d');
334 if dim<>5 then pause,end
335 [Z,dim]=schur(A,E,'d');
336 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
337
338 [Q,Z1,dim]=schur(A,E,'d');
339 if Err(Z1-Z)>10*%eps then pause,end
340 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
341 if dim<>5 then pause,end
342
343 [As,Es,Z,dim]=schur(A,E,'d');
344 if dim<>5 then pause,end
345 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
346 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
347 if Err(As-Q'*A*Z) >200*%eps then pause,end
348 if Err(Es-Q'*E*Z) >200*%eps then pause,end
349
350 //ordered sel
351 clear sel
352 function t=sel(Alpha,Beta)
353     t = real(Alpha) > -0.2 * real(Beta);
354 endfunction
355
356 dim=schur(A,E,sel);
357 if dim<>2 then pause,end
358 [Z,dim]=schur(A,E,sel);
359 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
360
361 [Q,Z1,dim]=schur(A,E,sel);
362 if Err(Z1-Z)>10*%eps then pause,end
363 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
364 if dim<>2 then pause,end
365
366 [As,Es,Z,dim]=schur(A,E,sel);
367 if dim<>2 then pause,end
368 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
369 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
370 if Err(As-Q'*A*Z) >200*%eps then pause,end
371 if Err(Es-Q'*E*Z) >200*%eps then pause,end
372
373
374 //----Complex------------
375 A=testmat1(1+%i,5);E=testmat1(-2-3*%i,5) ;
376 [As,Es,Q,Z]=schur(A,E);
377 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
378 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
379 if Err(As-Q'*A*Z) >200*%eps then pause,end
380 if Err(Es-Q'*E*Z) >200*%eps then pause,end
381
382 [As1,Es1]=schur(A,E);
383 if Err(As1-As)>10*%eps then pause,end
384 if Err(Es1-Es)>10*%eps then pause,end
385
386 // Ordered 'c'
387 dim=schur(A,E,'c');
388 if dim<>5 then pause,end
389 [Z,dim]=schur(A,E,'c');
390 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
391
392 [Q,Z1,dim]=schur(A,E,'c');
393 if Err(Z1-Z)>10*%eps then pause,end
394 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
395 if dim<>5 then pause,end
396
397 [As,Es,Z,dim]=schur(A,E,'c');
398 if dim<>5 then pause,end
399 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
400 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
401 if Err(As-Q'*A*Z) >200*%eps then pause,end
402 if Err(Es-Q'*E*Z) >200*%eps then pause,end
403
404
405 // Ordered 'd'
406 dim=schur(A,E,'d');
407 if dim<>5 then pause,end
408 [Z,dim]=schur(A,E,'d');
409 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
410
411 [Q,Z1,dim]=schur(A,E,'d');
412 if Err(Z1-Z)>10*%eps then pause,end
413 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
414 if dim<>5 then pause,end
415
416 [As,Es,Z,dim]=schur(A,E,'d');
417 if dim<>5 then pause,end
418 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
419 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
420 if Err(As-Q'*A*Z) >200*%eps then pause,end
421 if Err(Es-Q'*E*Z) >200*%eps then pause,end
422
423 //ordered sel
424 clear sel
425 function t=sel(Alpha,Beta),t=imag(Alpha)>0 ,endfunction
426
427 dim=schur(A,E,sel);
428 if dim<>3 then pause,end
429 [Z,dim]=schur(A,E,sel);
430 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
431
432 [Q,Z1,dim]=schur(A,E,sel);
433 if Err(Z1-Z)>10*%eps then pause,end
434 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
435 if dim<>3 then pause,end
436
437 [As,Es,Z,dim]=schur(A,E,sel);
438 if dim<>3 then pause,end
439 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
440 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
441 if Err(As-Q'*A*Z) >200*%eps then pause,end
442 if Err(Es-Q'*E*Z) >200*%eps then pause,end
443
444 //Large dimension
445
446 //----Real------------
447 A=testmat1(1,50);E=testmat1(-2,50) ;
448 [As,Es,Q,Z]=schur(A,E);
449 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
450 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
451 if Err(As-Q'*A*Z) >200*%eps then pause,end
452 if Err(Es-Q'*E*Z) >200*%eps then pause,end
453
454 [As1,Es1]=schur(A,E);
455 if Err(As1-As)>10*%eps then pause,end
456 if Err(Es1-Es)>10*%eps then pause,end
457
458 // Ordered 'c'
459 dim=schur(A,E,'c');
460 if dim<>50 then pause,end
461 [Z,dim]=schur(A,E,'c');
462 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
463
464 [Q,Z1,dim]=schur(A,E,'c');
465 if Err(Z1-Z)>10*%eps then pause,end
466 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
467 if dim<>50 then pause,end
468
469 [As,Es,Z,dim]=schur(A,E,'c');
470 if dim<>50 then pause,end
471 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
472 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
473 if Err(As-Q'*A*Z) >200*%eps then pause,end
474 if Err(Es-Q'*E*Z) >200*%eps then pause,end
475 // Ordered 'd'
476 dim=schur(A,E,'d');
477 if dim<>50 then pause,end
478 [Z,dim]=schur(A,E,'d');
479 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
480
481 [Q,Z1,dim]=schur(A,E,'d');
482 if Err(Z1-Z)>10*%eps then pause,end
483 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
484 if dim<>50 then pause,end
485
486 [As,Es,Z,dim]=schur(A,E,'d');
487 if dim<>50 then pause,end
488 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
489 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
490 if Err(As-Q'*A*Z) >200*%eps then pause,end
491 if Err(Es-Q'*E*Z) >200*%eps then pause,end
492
493 //ordered sel
494 clear sel
495 function t=sel(Alpha,Beta)
496         t = real(Alpha) > -0.2 * real(Beta);
497 endfunction
498
499 dim=schur(A,E,sel);
500 if dim<>12 then pause,end
501 [Z,dim]=schur(A,E,sel);
502 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
503
504 [Q,Z1,dim]=schur(A,E,sel);
505 if Err(Z1-Z)>10*%eps then pause,end
506 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
507 if dim<>12 then pause,end
508
509 [As,Es,Z,dim]=schur(A,E,sel);
510 if dim<>12 then pause,end
511 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
512 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
513 if Err(As-Q'*A*Z) >200*%eps then pause,end
514 if Err(Es-Q'*E*Z) >200*%eps then pause,end
515
516 // Lib part
517 cd TMPDIR;
518
519 // equal to schur(A, E, 'c');
520 C=[ 'extern double dlamch_(char *CMACH, unsigned long int);' // dlamch_ == C2F(dlamch)
521     ''
522     'int mytest5(double* _real, double* _img, double* _beta)'
523     '{'
524     '    double dblP = dlamch_((char*)""p"", 1L);'
525     '    int iTest1 = (*_real < 0 && *_beta > 0);'
526     '    int iTest2 = (*_real > 0 && *_beta < 0);'
527     '    int iTest3 = (fabs(*_beta) > fabs(*_real) * dblP);'
528     ''
529     '    return (iTest1 || iTest2 && iTest3);'
530     '}'];
531 mputl(C,TMPDIR+'/mytest.c');
532 ulink();
533 lp=ilib_for_link('mytest5','mytest.c',[],'c');
534 exec loader.sce;
535
536 dim=schur(A,E,'mytest5');
537 if dim<>50 then pause,end
538 [Z,dim]=schur(A,E,'mytest5');
539 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
540
541 [Q,Z1,dim]=schur(A,E,'mytest5');
542 if Err(Z1-Z)>10*%eps then pause,end
543 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
544 if dim<>50 then pause,end
545
546 [As,Es,Z,dim]=schur(A,E,'mytest5');
547 if dim<>50 then pause,end
548 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
549 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
550 if Err(As-Q'*A*Z) >200*%eps then pause,end
551 if Err(Es-Q'*E*Z) >200*%eps then pause,end
552
553 // equal to schur(A, E, 'd');
554 C=[ 'extern double dpythags(double,double);' // YaSp function
555     ''
556     'int mytest6(double* _real, double* _img, double* _beta)'
557     '{'
558     '    double dblPythag =  dpythags(*_real, *_img);'
559     ''
560     '    return (dblPythag < fabs(*_beta));'
561     '}'];
562 mputl(C,TMPDIR+'/mytest.c');
563 ulink();
564 lp=ilib_for_link('mytest6','mytest.c', libName,'c');
565 exec loader.sce;
566
567 dim=schur(A,E,'mytest6');
568 if dim<>50 then pause,end
569
570 [Z,dim]=schur(A,E,'mytest6');
571 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
572
573 [Q,Z1,dim]=schur(A,E,'mytest6');
574 if Err(Z1-Z)>10*%eps then pause,end
575 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
576 if dim<>50 then pause,end
577
578 [As,Es,Z,dim]=schur(A,E,'mytest6');
579 if dim<>50 then pause,end
580 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
581 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
582 if Err(As-Q'*A*Z) >200*%eps then pause,end
583 if Err(Es-Q'*E*Z) >200*%eps then pause,end
584
585 //----Complex------------
586 A=testmat1(1+%i,50);E=testmat1(-2-3*%i,50) ;
587 [As,Es,Q,Z]=schur(A,E);
588 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
589 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
590 if Err(As-Q'*A*Z) >1000*%eps then pause,end
591 if Err(Es-Q'*E*Z) >1000*%eps then pause,end
592
593 [As1,Es1]=schur(A,E);
594 if Err(As1-As)>10*%eps then pause,end
595 if Err(Es1-Es)>10*%eps then pause,end
596
597 // Ordered 'c'
598 dim=schur(A,E,'c');
599 if dim<>50 then pause,end
600 [Z,dim]=schur(A,E,'c');
601 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
602
603 [Q,Z1,dim]=schur(A,E,'c');
604 if Err(Z1-Z)>10*%eps then pause,end
605 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
606 if dim<>50 then pause,end
607
608 [As,Es,Z,dim]=schur(A,E,'c');
609 if dim<>50 then pause,end
610 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
611 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
612 if Err(As-Q'*A*Z) >1000*%eps then pause,end
613 if Err(Es-Q'*E*Z) >1000*%eps then pause,end
614 // Ordered 'd'
615 dim=schur(A,E,'d');
616 if dim<>50 then pause,end
617 [Z,dim]=schur(A,E,'d');
618 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
619
620 [Q,Z1,dim]=schur(A,E,'d');
621 if Err(Z1-Z)>10*%eps then pause,end
622 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
623 if dim<>50 then pause,end
624
625 [As,Es,Z,dim]=schur(A,E,'d');
626 if dim<>50 then pause,end
627 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
628 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
629 if Err(As-Q'*A*Z) >1000*%eps then pause,end
630 if Err(Es-Q'*E*Z) >1000*%eps then pause,end
631
632 //ordered sel
633 clear sel
634 function t=sel(Alpha,Beta),t=imag(Alpha)>0 ,endfunction
635
636 dim=schur(A,E,sel);
637 if dim<>32 then pause,end
638 [Z,dim]=schur(A,E,sel);
639 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
640
641 [Q,Z1,dim]=schur(A,E,sel);
642 if Err(Z1-Z)>10*%eps then pause,end
643 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
644 if dim<>32 then pause,end
645
646 [As,Es,Z,dim]=schur(A,E,sel);
647 if dim<>32 then pause,end
648 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
649 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
650 if Err(As-Q'*A*Z) >1000*%eps then pause,end
651 if Err(Es-Q'*E*Z) >1000*%eps then pause,end
652
653 //Lib part
654 // equal to schur(A, E, 'c');
655 C=[ '#include ""doublecomplex.h"";'
656     ''
657     'int mytest7(doublecomplex* _complex)'
658     '{'
659     '    return (_complex->r < 0);'
660     '}'];
661 mputl(C,TMPDIR+'/mytest.c');
662 ulink();
663 lp=ilib_for_link('mytest7','mytest.c',[],'c');
664 exec loader.sce;
665
666 dim=schur(A,E,'mytest7');
667 if dim<>50 then pause,end
668
669 [Z,dim]=schur(A,E,'mytest7');
670 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
671
672 [Q,Z1,dim]=schur(A,E,'mytest7');
673 if Err(Z1-Z)>10*%eps then pause,end
674 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
675 if dim<>50 then pause,end
676
677 [As,Es,Z,dim]=schur(A,E,'mytest7');
678 if dim<>50 then pause,end
679 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
680 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
681 if Err(As-Q'*A*Z) >1000*%eps then pause,end
682 if Err(Es-Q'*E*Z) >1000*%eps then pause,end
683
684 // equal to schur(A, E, 'd');
685 C=[ '#include ""doublecomplex.h"";'
686     ''
687     'extern double dpythags(double,double);' // YaSp function
688     ''
689     'int mytest8(doublecomplex* _alpha, doublecomplex* _beta)'
690     '{'
691     '    double dblP1 = dpythags(_alpha->r, _alpha->i);'
692     '    double dblP2 = dpythags(_beta->r, _beta->i);'
693     ''
694     '    return (dblP1 <  dblP2);'
695     '}'];
696 mputl(C,TMPDIR+'/mytest.c');
697 ulink();
698 lp=ilib_for_link('mytest8','mytest.c', libName,'c');
699 exec loader.sce;
700
701 dim=schur(A,E,'mytest8');
702 if dim<>50 then pause,end
703
704 [Z,dim]=schur(A,E,'mytest8');
705 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
706
707 [Q,Z1,dim]=schur(A,E,'mytest8');
708 if Err(Z1-Z)>10*%eps then pause,end
709 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
710 if dim<>50 then pause,end
711
712 [As,Es,Z,dim]=schur(A,E,'mytest8');
713 if dim<>50 then pause,end
714 if Err(Q*Q'-eye(Q)) >200*%eps then pause,end
715 if Err(Z*Z'-eye(Z)) >200*%eps then pause,end
716 if Err(As-Q'*A*Z) >1000*%eps then pause,end
717 if Err(Es-Q'*E*Z) >1000*%eps then pause,end
718