linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / tests / unit_tests / schur.dia.ref
1
2 // =============================================================================
3
4 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
5
6 // Copyright (C) ????-2008 - INRIA Michael Baudin
7
8 //
9
10 //  This file is distributed under the same license as the Scilab package.
11
12 // =============================================================================
13
14 function r=Err(x)
15         r=norm(x,1)
16 endfunction
17 rand('normal')
18
19 ilib_verbose(0);
20
21 libName = [];
22
23 if getos() == "Windows" then
24     libName = SCI + "/bin/elem_func";
25 end
26 //define tools
27
28 function A=testmat1(a,n)
29         //eigen values are given by a dilation of nth roots of 1
30         A=diag(a*ones(1,n-1),1)+diag((1/a)*ones(1,n-1),-1)
31         A(1,n)=1/a;A(n,1)=a
32 endfunction
33 //==========================================================================
34
35 //==============================    schur     ============================== 
36
37 //==========================================================================
38
39 clear sel
40
41 function t=sel(R),t=real(R)<0 ,endfunction
42
43 //Empty matrix
44
45 A=[];
46
47 if schur(A)<>[] then bugmes();quit;end
48
49 if schur(A,'real')<>[] then bugmes();quit;end
50
51 if schur(A,'complex')<>[] then bugmes();quit;end
52
53 if schur(A,'c')<>[] then bugmes();quit;end
54
55 if schur(A,'d')<>[] then bugmes();quit;end
56
57 if schur(A,sel)<>[] then bugmes();quit;end
58
59 [U,S]=schur(A);
60
61 if U<>[]|S<>[] then bugmes();quit;end
62
63 [U,S]=schur(A,'real');
64
65 if U<>[]|S<>[] then bugmes();quit;end
66
67 [U,S]=schur(A,'complex');
68
69 if U<>[]|S<>[] then bugmes();quit;end
70
71 [U,N]=schur(A,'c');
72
73 if U<>[]|N<>0 then bugmes();quit;end
74
75 [U,N]=schur(A,'d');
76
77 if U<>[]|N<>0 then bugmes();quit;end
78
79 [U,N]=schur(A,sel);
80
81 if U<>[]|N<>0 then bugmes();quit;end
82
83 [U,N,S]=schur(A,'c');
84
85 if U<>[]|N<>0|S<>[] then bugmes();quit;end
86
87 [U,N,S]=schur(A,'d');
88
89 if U<>[]|N<>0|S<>[] then bugmes();quit;end
90
91 [U,N,S]=schur(A,sel);
92
93 if U<>[]|N<>0|S<>[] then bugmes();quit;end
94
95 //Rectangular matrix
96
97 if execstr('schur(rand(2,3))','errcatch')==0 then bugmes();quit;end
98
99 if execstr('[U,S]=schur(rand(2,3))','errcatch')==0 then bugmes();quit;end
100
101 if execstr('schur(rand(2,3)+%i*eye())','errcatch')==0 then bugmes();quit;end
102
103 if execstr('[U,S]=schur(rand(2,3)+%i*eye())','errcatch')==0 then bugmes();quit;end
104
105 //Small dimension
106
107 A=testmat1(3,5);Ac=testmat1(3+%i,5);
108
109 //Real
110
111 [U,S]=schur(A);
112
113 if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
114
115 if Err(U*S*U'-A)>200*%eps then bugmes();quit;end
116
117 if Err(schur(A)-S) >%eps then bugmes();quit;end
118
119 [U,S]=schur(A,'real');
120
121 if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
122
123 if Err(U*S*U'-A)>200*%eps then bugmes();quit;end
124
125 if Err(schur(A)-S) >%eps then bugmes();quit;end
126
127 [U,S]=schur(A,'complex');
128
129 if Err(triu(S)-S)>%eps then bugmes();quit;end
130
131 if Err(U*S*U'-A)>200*%eps then bugmes();quit;end
132
133 if Err(schur(A,'complex')-S) >%eps then bugmes();quit;end
134
135 [U,n]=schur(A,'c');S=U'*A*U;
136
137 if n<>2 then bugmes();quit;end
138
139 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
140
141 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
142
143 [U,n]=schur(A,'d');S=U'*A*U;
144
145 if n<>0 then bugmes();quit;end
146
147 if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
148
149 [U,n]=schur(A,sel);S=U'*A*U;
150
151 if n<>2 then bugmes();quit;end
152
153 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
154
155 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
156
157 //Complex
158
159 [U,S]=schur(Ac);
160
161 if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
162
163 if Err(U*S*U'-Ac)>200*%eps then bugmes();quit;end
164
165 if Err(schur(Ac)-S) >%eps then bugmes();quit;end
166
167 [U,S]=schur(Ac,'complex');
168
169 if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
170
171 if Err(U*S*U'-Ac)>200*%eps then bugmes();quit;end
172
173 if Err(schur(Ac)-S) >%eps then bugmes();quit;end
174
175 [U,n]=schur(Ac,'c');S=U'*Ac*U;
176
177 if n<>3 then bugmes();quit;end
178
179 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
180
181 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
182
183 [U,n]=schur(Ac,'d');S=U'*A*U;
184
185 if n<>0 then bugmes();quit;end
186
187 if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
188
189 [U,n]=schur(Ac,sel);S=U'*Ac*U;
190
191 if n<>3 then bugmes();quit;end
192
193 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
194
195 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
196
197 //Large dimension
198
199 A=testmat1(3,50);Ac=testmat1(3+%i,50);
200
201 //Real
202
203 [U,S]=schur(A);
204
205 if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
206
207 if Err(U*S*U'-A)>1000*%eps then bugmes();quit;end
208
209 if Err(schur(A)-S) >%eps then bugmes();quit;end
210
211 [U,S]=schur(A,'real');
212
213 if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
214
215 if Err(U*S*U'-A)>1000*%eps then bugmes();quit;end
216
217 if Err(schur(A)-S) >%eps then bugmes();quit;end
218
219 [U,S]=schur(A,'complex');
220
221 if Err(triu(S)-S)>%eps then bugmes();quit;end
222
223 if Err(U*S*U'-A)>1000*%eps then bugmes();quit;end
224
225 if Err(schur(A,'complex')-S) >%eps then bugmes();quit;end
226
227 [U,n]=schur(A,'c');S=U'*A*U;
228
229 if n<>25 then bugmes();quit;end
230
231 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
232
233 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
234
235 [U,n]=schur(A,'d');S=U'*A*U;
236
237 if n<>0 then bugmes();quit;end
238
239 if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
240
241 [U,n]=schur(A,sel);S=U'*A*U;
242
243 if n<>25 then bugmes();quit;end
244
245 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
246
247 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
248
249 //Complex
250
251 [U,S]=schur(Ac);
252
253 if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
254
255 if Err(U*S*U'-Ac)>1000*%eps then bugmes();quit;end
256
257 if Err(schur(Ac)-S) >%eps then bugmes();quit;end
258
259 [U,S]=schur(Ac,'complex');
260
261 if Err(triu(S,-1)-S)>%eps then bugmes();quit;end
262
263 if Err(U*S*U'-Ac)>1000*%eps then bugmes();quit;end
264
265 if Err(schur(Ac)-S) >%eps then bugmes();quit;end
266
267 [U,n]=schur(Ac,'c');S=U'*Ac*U;
268
269 if n<>25 then bugmes();quit;end
270
271 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
272
273 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
274
275 [U,n]=schur(Ac,'d');S=U'*Ac*U;
276
277 if n<>0 then bugmes();quit;end
278
279 if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
280
281 [U,n]=schur(Ac,sel);S=U'*Ac*U;
282
283 if n<>25 then bugmes();quit;end
284
285 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
286
287 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
288
289 // Lib part
290
291 cd TMPDIR;
292
293 // equal to schur(A, 'c');
294
295 C=[ 'int mytest1(double* _real, double* _img)'
296     '{' 
297     '    return *_real < 0;'
298     '}'];
299
300 mputl(C,TMPDIR+'/mytest.c');
301
302 ulink();
303
304 lp=ilib_for_link('mytest1','mytest.c',[],'c');
305
306 exec loader.sce;
307
308 [U,n]=schur(A,'mytest1');S=U'*A*U;
309
310 if n<>25 then bugmes();quit;end
311
312 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
313
314 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
315
316 // equal to schur(A, 'd');
317
318 C=[ 'extern double dpythags(double,double);' // YaSp function
319     ''
320     'int mytest2(double* _real, double* _img)'
321     '{' 
322     '    return dpythags(*_real, *_img) < 1;'
323     '}'];
324
325 mputl(C,TMPDIR+'/mytest.c');
326
327 ulink();
328
329 lp=ilib_for_link('mytest2','mytest.c', libName,'c');
330
331 exec loader.sce;
332
333 [U,n]=schur(A,'mytest2');S=U'*A*U;
334
335 if n<>0 then bugmes();quit;end
336
337 if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
338
339 // equal to schur(Ac, 'c');
340
341 C=[ '#include ""doublecomplex.h""'
342     ''
343     'int mytest3(doublecomplex* _complex)'
344     '{'
345     '    return _complex->r < 0 ? 1 : 0;'
346     '}'];
347
348 mputl(C,TMPDIR+'/mytest.c');
349
350 ulink();
351
352 lp=ilib_for_link('mytest3','mytest.c',[],'c');
353
354 exec loader.sce;
355
356 [U,n]=schur(Ac, 'mytest3');S=U'*Ac*U;
357
358 if n<>25 then bugmes();quit;end
359
360 if or(real(spec(S(1:n,1:n)))>=0) then bugmes();quit;end
361
362 if or(real(spec(S(n+1:$,n+1:$)))<0) then bugmes();quit;end
363
364 // equal to schur(Ac, 'd');
365
366 C=[ '#include ""doublecomplex.h""'
367     ''
368     'extern double dpythags(double,double);' // YaSp function
369     ''
370     'int mytest4(doublecomplex* _complex)'
371     '{'
372     '    if(dpythags(_complex->r, _complex->i) < 1)'
373     '    {'
374     '        return 1;'
375     '    }'
376     '    else'
377     '    {'
378     '        return 0;'
379     '    }'
380     '}'];
381
382 mputl(C,TMPDIR+'/mytest.c');
383
384 ulink();
385
386 lp=ilib_for_link('mytest4','mytest.c', libName,'c');
387
388 exec loader.sce;
389
390 [U,n]=schur(Ac,'mytest4');S=U'*Ac*U;
391
392 if n<>0 then bugmes();quit;end
393
394 if or(abs(spec(S(n+1:$,n+1:$)))<1) then bugmes();quit;end
395
396 //==========================================================================
397
398 //==============================    schur part II   ======================== 
399
400 //==========================================================================
401
402 //Empty matrix
403
404 [As,Es]=schur([],[]);
405
406 if As<>[]|Es<>[] then bugmes();quit;end
407
408 [As,dim]=schur([],[],'c');
409
410 if As<>[]|dim<>0 then bugmes();quit;end
411
412 [As,dim]=schur([],[],'d');
413
414 if As<>[]|dim<>0 then bugmes();quit;end
415
416 [As,dim]=schur([],[],sel);
417
418 if As<>[]|dim<>0 then bugmes();quit;end
419
420 [As,Es,Q,Z]=schur([],[]);
421
422 if As<>[]|Es<>[]|Q<>[]|Z<>[] then bugmes();quit;end
423
424 [As,Es,dim]=schur([],[],'c');
425
426 if As<>[]|Es<>[]|dim<>0 then bugmes();quit;end
427
428 [As,Es,dim]=schur([],[],'d');
429
430 if As<>[]|Es<>[]|dim<>0 then bugmes();quit;end
431
432 [As,Es,dim]=schur([],[],sel);
433
434 if As<>[]|Es<>[]|dim<>0 then bugmes();quit;end
435
436 [Z,dim]=schur([],[],'c');
437
438 if Z<>[]|dim<>0 then bugmes();quit;end
439
440 [Z,dim]=schur([],[],'d');
441
442 if Z<>[]|dim<>0 then bugmes();quit;end
443
444 [Z,dim]=schur([],[],sel);
445
446 if Z<>[]|dim<>0 then bugmes();quit;end
447
448 //Rectangular matrix
449
450 if execstr('[As,Es]=schur(rand(2,3),rand(2,3))','errcatch')==0 then  bugmes();quit;end
451
452 if execstr('[As,Es,Q,Z]=schur(rand(2,3),rand(2,3))','errcatch')==0 then  bugmes();quit;end
453
454 if execstr('[As,Es,dim]=schur(rand(2,3),rand(2,3),''c'')','errcatch')==0 then  bugmes();quit;end
455
456 if execstr('[Z,dim]=schur(rand(2,3),rand(2,3),sel)','errcatch')==0 then  bugmes();quit;end
457
458 //Small dimension
459
460 //----Real------------
461
462 A=testmat1(1,5);E=testmat1(-2,5) 
463 [As,Es,Q,Z]=schur(A,E);
464
465 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
466
467 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
468
469 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
470
471 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
472
473 [As1,Es1]=schur(A,E);
474
475 if Err(As1-As)>10*%eps then bugmes();quit;end
476
477 if Err(Es1-Es)>10*%eps then bugmes();quit;end
478
479 // Ordered 'c'
480
481 dim=schur(A,E,'c');
482
483 if dim<>5 then bugmes();quit;end
484
485 [Z,dim]=schur(A,E,'c');
486
487 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
488
489 [Q,Z1,dim]=schur(A,E,'c');
490
491 if Err(Z1-Z)>10*%eps then bugmes();quit;end
492
493 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
494
495 if dim<>5 then bugmes();quit;end
496
497 [As,Es,Z,dim]=schur(A,E,'c');
498
499 if dim<>5 then bugmes();quit;end
500
501 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
502
503 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
504
505 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
506
507 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
508
509 // Ordered 'd'
510
511 dim=schur(A,E,'d');
512
513 if dim<>5 then bugmes();quit;end
514
515 [Z,dim]=schur(A,E,'d');
516
517 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
518
519 [Q,Z1,dim]=schur(A,E,'d');
520
521 if Err(Z1-Z)>10*%eps then bugmes();quit;end
522
523 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
524
525 if dim<>5 then bugmes();quit;end
526
527 [As,Es,Z,dim]=schur(A,E,'d');
528
529 if dim<>5 then bugmes();quit;end
530
531 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
532
533 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
534
535 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
536
537 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
538
539 //ordered sel
540
541 clear sel
542
543 function t=sel(Alpha,Beta)
544     t = real(Alpha) > -0.2 * real(Beta);
545 endfunction
546 dim=schur(A,E,sel);
547
548 if dim<>2 then bugmes();quit;end
549
550 [Z,dim]=schur(A,E,sel);
551
552 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
553
554 [Q,Z1,dim]=schur(A,E,sel);
555
556 if Err(Z1-Z)>10*%eps then bugmes();quit;end
557
558 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
559
560 if dim<>2 then bugmes();quit;end
561
562 [As,Es,Z,dim]=schur(A,E,sel);
563
564 if dim<>2 then bugmes();quit;end
565
566 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
567
568 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
569
570 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
571
572 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
573
574 //----Complex------------
575
576 A=testmat1(1+%i,5);E=testmat1(-2-3*%i,5) 
577 [As,Es,Q,Z]=schur(A,E);
578
579 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
580
581 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
582
583 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
584
585 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
586
587 [As1,Es1]=schur(A,E);
588
589 if Err(As1-As)>10*%eps then bugmes();quit;end
590
591 if Err(Es1-Es)>10*%eps then bugmes();quit;end
592
593 // Ordered 'c'
594
595 dim=schur(A,E,'c');
596
597 if dim<>5 then bugmes();quit;end
598
599 [Z,dim]=schur(A,E,'c');
600
601 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
602
603 [Q,Z1,dim]=schur(A,E,'c');
604
605 if Err(Z1-Z)>10*%eps then bugmes();quit;end
606
607 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
608
609 if dim<>5 then bugmes();quit;end
610
611 [As,Es,Z,dim]=schur(A,E,'c');
612
613 if dim<>5 then bugmes();quit;end
614
615 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
616
617 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
618
619 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
620
621 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
622
623 // Ordered 'd'
624
625 dim=schur(A,E,'d');
626
627 if dim<>5 then bugmes();quit;end
628
629 [Z,dim]=schur(A,E,'d');
630
631 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
632
633 [Q,Z1,dim]=schur(A,E,'d');
634
635 if Err(Z1-Z)>10*%eps then bugmes();quit;end
636
637 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
638
639 if dim<>5 then bugmes();quit;end
640
641 [As,Es,Z,dim]=schur(A,E,'d');
642
643 if dim<>5 then bugmes();quit;end
644
645 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
646
647 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
648
649 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
650
651 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
652
653 //ordered sel
654
655 clear sel
656
657 function t=sel(Alpha,Beta),t=imag(Alpha)>0 ,endfunction
658
659 dim=schur(A,E,sel);
660
661 if dim<>3 then bugmes();quit;end
662
663 [Z,dim]=schur(A,E,sel);
664
665 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
666
667 [Q,Z1,dim]=schur(A,E,sel);
668
669 if Err(Z1-Z)>10*%eps then bugmes();quit;end
670
671 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
672
673 if dim<>3 then bugmes();quit;end
674
675 [As,Es,Z,dim]=schur(A,E,sel);
676
677 if dim<>3 then bugmes();quit;end
678
679 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
680
681 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
682
683 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
684
685 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
686
687 //Large dimension
688
689 //----Real------------
690
691 A=testmat1(1,50);E=testmat1(-2,50) 
692 [As,Es,Q,Z]=schur(A,E);
693
694 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
695
696 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
697
698 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
699
700 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
701
702 [As1,Es1]=schur(A,E);
703
704 if Err(As1-As)>10*%eps then bugmes();quit;end
705
706 if Err(Es1-Es)>10*%eps then bugmes();quit;end
707
708 // Ordered 'c'
709
710 dim=schur(A,E,'c');
711
712 if dim<>50 then bugmes();quit;end
713
714 [Z,dim]=schur(A,E,'c');
715
716 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
717
718 [Q,Z1,dim]=schur(A,E,'c');
719
720 if Err(Z1-Z)>10*%eps then bugmes();quit;end
721
722 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
723
724 if dim<>50 then bugmes();quit;end
725
726 [As,Es,Z,dim]=schur(A,E,'c');
727
728 if dim<>50 then bugmes();quit;end
729
730 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
731
732 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
733
734 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
735
736 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
737
738 // Ordered 'd'
739
740 dim=schur(A,E,'d');
741
742 if dim<>50 then bugmes();quit;end
743
744 [Z,dim]=schur(A,E,'d');
745
746 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
747
748 [Q,Z1,dim]=schur(A,E,'d');
749
750 if Err(Z1-Z)>10*%eps then bugmes();quit;end
751
752 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
753
754 if dim<>50 then bugmes();quit;end
755
756 [As,Es,Z,dim]=schur(A,E,'d');
757
758 if dim<>50 then bugmes();quit;end
759
760 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
761
762 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
763
764 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
765
766 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
767
768 //ordered sel
769
770 clear sel
771
772 function t=sel(Alpha,Beta)
773         t = real(Alpha) > -0.2 * real(Beta);
774 endfunction
775 dim=schur(A,E,sel);
776
777 if dim<>12 then bugmes();quit;end
778
779 [Z,dim]=schur(A,E,sel);
780
781 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
782
783 [Q,Z1,dim]=schur(A,E,sel);
784
785 if Err(Z1-Z)>10*%eps then bugmes();quit;end
786
787 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
788
789 if dim<>12 then bugmes();quit;end
790
791 [As,Es,Z,dim]=schur(A,E,sel);
792
793 if dim<>12 then bugmes();quit;end
794
795 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
796
797 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
798
799 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
800
801 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
802
803 // Lib part
804
805 cd TMPDIR;
806
807 // equal to schur(A, E, 'c');
808
809 C=[ 'extern double dlamch_(char *CMACH, unsigned long int);' // dlamch_ == C2F(dlamch)
810     ''
811     'int mytest5(double* _real, double* _img, double* _beta)'
812     '{'
813     '    double dblP = dlamch_((char*)""p"", 1L);'
814     '    int iTest1 = (*_real < 0 && *_beta > 0);'
815     '    int iTest2 = (*_real > 0 && *_beta < 0);'
816     '    int iTest3 = (fabs(*_beta) > fabs(*_real) * dblP);'
817     ''
818     '    return (iTest1 || iTest2 && iTest3);'
819     '}'];
820
821 mputl(C,TMPDIR+'/mytest.c');
822
823 ulink();
824
825 lp=ilib_for_link('mytest5','mytest.c',[],'c');
826
827 exec loader.sce;
828
829 dim=schur(A,E,'mytest5');
830
831 if dim<>50 then bugmes();quit;end
832
833 [Z,dim]=schur(A,E,'mytest5');
834
835 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
836
837 [Q,Z1,dim]=schur(A,E,'mytest5');
838
839 if Err(Z1-Z)>10*%eps then bugmes();quit;end
840
841 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
842
843 if dim<>50 then bugmes();quit;end
844
845 [As,Es,Z,dim]=schur(A,E,'mytest5');
846
847 if dim<>50 then bugmes();quit;end
848
849 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
850
851 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
852
853 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
854
855 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
856
857 // equal to schur(A, E, 'd');
858
859 C=[ 'extern double dpythags(double,double);' // YaSp function
860     ''
861     'int mytest6(double* _real, double* _img, double* _beta)'
862     '{'
863     '    double dblPythag =  dpythags(*_real, *_img);'
864     ''
865     '    return (dblPythag < fabs(*_beta));'
866     '}'];
867
868 mputl(C,TMPDIR+'/mytest.c');
869
870 ulink();
871
872 lp=ilib_for_link('mytest6','mytest.c', libName,'c');
873
874 exec loader.sce;
875
876 dim=schur(A,E,'mytest6');
877
878 if dim<>50 then bugmes();quit;end
879
880 [Z,dim]=schur(A,E,'mytest6');
881
882 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
883
884 [Q,Z1,dim]=schur(A,E,'mytest6');
885
886 if Err(Z1-Z)>10*%eps then bugmes();quit;end
887
888 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
889
890 if dim<>50 then bugmes();quit;end
891
892 [As,Es,Z,dim]=schur(A,E,'mytest6');
893
894 if dim<>50 then bugmes();quit;end
895
896 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
897
898 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
899
900 if Err(As-Q'*A*Z) >200*%eps then bugmes();quit;end
901
902 if Err(Es-Q'*E*Z) >200*%eps then bugmes();quit;end
903
904 //----Complex------------
905
906 A=testmat1(1+%i,50);E=testmat1(-2-3*%i,50) 
907 [As,Es,Q,Z]=schur(A,E);
908
909 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
910
911 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
912
913 if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
914
915 if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
916
917 [As1,Es1]=schur(A,E);
918
919 if Err(As1-As)>10*%eps then bugmes();quit;end
920
921 if Err(Es1-Es)>10*%eps then bugmes();quit;end
922
923 // Ordered 'c'
924
925 dim=schur(A,E,'c');
926
927 if dim<>50 then bugmes();quit;end
928
929 [Z,dim]=schur(A,E,'c');
930
931 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
932
933 [Q,Z1,dim]=schur(A,E,'c');
934
935 if Err(Z1-Z)>10*%eps then bugmes();quit;end
936
937 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
938
939 if dim<>50 then bugmes();quit;end
940
941 [As,Es,Z,dim]=schur(A,E,'c');
942
943 if dim<>50 then bugmes();quit;end
944
945 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
946
947 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
948
949 if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
950
951 if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
952
953 // Ordered 'd'
954
955 dim=schur(A,E,'d');
956
957 if dim<>50 then bugmes();quit;end
958
959 [Z,dim]=schur(A,E,'d');
960
961 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
962
963 [Q,Z1,dim]=schur(A,E,'d');
964
965 if Err(Z1-Z)>10*%eps then bugmes();quit;end
966
967 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
968
969 if dim<>50 then bugmes();quit;end
970
971 [As,Es,Z,dim]=schur(A,E,'d');
972
973 if dim<>50 then bugmes();quit;end
974
975 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
976
977 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
978
979 if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
980
981 if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
982
983 //ordered sel
984
985 clear sel
986
987 function t=sel(Alpha,Beta),t=imag(Alpha)>0 ,endfunction
988
989 dim=schur(A,E,sel);
990
991 if dim<>32 then bugmes();quit;end
992
993 [Z,dim]=schur(A,E,sel);
994
995 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
996
997 [Q,Z1,dim]=schur(A,E,sel);
998
999 if Err(Z1-Z)>10*%eps then bugmes();quit;end
1000
1001 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
1002
1003 if dim<>32 then bugmes();quit;end
1004
1005 [As,Es,Z,dim]=schur(A,E,sel);
1006
1007 if dim<>32 then bugmes();quit;end
1008
1009 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
1010
1011 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
1012
1013 if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
1014
1015 if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
1016
1017 //Lib part
1018
1019 // equal to schur(A, E, 'c');
1020
1021 C=[ '#include ""doublecomplex.h"";'
1022     ''
1023     'int mytest7(doublecomplex* _complex)'
1024     '{'
1025     '    return (_complex->r < 0);'
1026     '}'];
1027
1028 mputl(C,TMPDIR+'/mytest.c');
1029
1030 ulink();
1031
1032 lp=ilib_for_link('mytest7','mytest.c',[],'c');
1033
1034 exec loader.sce;
1035
1036 dim=schur(A,E,'mytest7');
1037
1038 if dim<>50 then bugmes();quit;end
1039
1040 [Z,dim]=schur(A,E,'mytest7');
1041
1042 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
1043
1044 [Q,Z1,dim]=schur(A,E,'mytest7');
1045
1046 if Err(Z1-Z)>10*%eps then bugmes();quit;end
1047
1048 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
1049
1050 if dim<>50 then bugmes();quit;end
1051
1052 [As,Es,Z,dim]=schur(A,E,'mytest7');
1053
1054 if dim<>50 then bugmes();quit;end
1055
1056 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
1057
1058 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
1059
1060 if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
1061
1062 if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
1063
1064 // equal to schur(A, E, 'd');
1065
1066 C=[ '#include ""doublecomplex.h"";'
1067     ''
1068     'extern double dpythags(double,double);' // YaSp function
1069     ''
1070     'int mytest8(doublecomplex* _alpha, doublecomplex* _beta)'
1071     '{'
1072     '    double dblP1 = dpythags(_alpha->r, _alpha->i);'
1073     '    double dblP2 = dpythags(_beta->r, _beta->i);'
1074     ''
1075     '    return (dblP1 <  dblP2);'
1076     '}'];
1077
1078 mputl(C,TMPDIR+'/mytest.c');
1079
1080 ulink();
1081
1082 lp=ilib_for_link('mytest8','mytest.c', libName,'c');
1083
1084 exec loader.sce;
1085
1086 dim=schur(A,E,'mytest8');
1087
1088 if dim<>50 then bugmes();quit;end
1089
1090 [Z,dim]=schur(A,E,'mytest8');
1091
1092 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
1093
1094 [Q,Z1,dim]=schur(A,E,'mytest8');
1095
1096 if Err(Z1-Z)>10*%eps then bugmes();quit;end
1097
1098 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
1099
1100 if dim<>50 then bugmes();quit;end
1101
1102 [As,Es,Z,dim]=schur(A,E,'mytest8');
1103
1104 if dim<>50 then bugmes();quit;end
1105
1106 if Err(Q*Q'-eye(Q)) >200*%eps then bugmes();quit;end
1107
1108 if Err(Z*Z'-eye(Z)) >200*%eps then bugmes();quit;end
1109
1110 if Err(As-Q'*A*Z) >1000*%eps then bugmes();quit;end
1111
1112 if Err(Es-Q'*E*Z) >1000*%eps then bugmes();quit;end
1113