linear_algebra plugged.
[scilab.git] / scilab / modules / linear_algebra / tests / unit_tests / bdiag.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),r=norm(x,1),endfunction
8 rand('normal')
9
10
11 //==========================================================================
12 //==============================   bdiag      ============================== 
13 //==========================================================================
14 if bdiag([])<>[] then pause,end
15 [ab,x]=bdiag([]);
16 if ab<>[]|x<>[] then pause,end
17 [ab,x,bs]=bdiag([]);
18 if ab<>[]|x<>[]|bs<>[] then pause,end
19 if execstr('bdiag([1 2;3 4;5 6])','errcatch')==0 then pause,end
20
21 //Small dimension
22 //---------------
23 //Real case
24 e=1.d-1;
25 A=[1 1  2 3 4 5
26    0 1  6 7 8 9
27    0 0  1 e 3 1
28    0 0 -e 1 5 9
29    0 0  0 0 2 e
30    0 0  0 0 0 3];
31 X1=[0.5,0.3,0,0.3,0.3,0.2;
32    1,0.6,0.5,0.1,0.7,0.4;
33    0.7,0.1,0.4,0.6,0.1,1;
34    0,0.6,0.2,0.3,0.4,0.5;
35    0.6,0.7,0.5,0.7,0.7,0.5;
36    0.3,0.3,0.4,0.5,0.9,0.6]
37 A=inv(X1)*A*X1;
38
39 Ab1=bdiag(A);
40 if or(triu(Ab1,-1)<>Ab1) then pause,end
41 [Ab2,X]=bdiag(A);
42 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
43 if Err(Ab2-inv(X)*A*X )>1d6*%eps then pause,end 
44
45 [Ab2,X,bs]=bdiag(A);
46 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
47 if Err(Ab2-inv(X)*A*X )>2.d-10 then pause,end 
48 if or(size(bs)<>[3,1]) then pause,end
49 if sum(bs)<>size(A,1) then pause,end
50 if or(bs<=0) then pause,end
51
52 [Ab2,X,bs]=bdiag(A,1);
53 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
54 if Err(Ab2-inv(X)*A*X )>2d-7 then pause,end 
55 if or(size(bs)<>[1,1]) then pause,end
56 if sum(bs)<>size(A,1) then pause,end
57 if or(bs<=0) then pause,end
58
59
60
61 //Complex case
62 e=1.d-1;
63 A=[1 1  2 3 4 5
64    0 1  6 7 8 9
65    0 0  1 e 3 1
66    0 0 -e 1 5 9
67    0 0  0 0 2 e
68    0 0  0 0 0 3];
69 X1=[0.5,0.3,0,0.3,0.3,0.2;
70    1,0.6,0.5,0.1,0.7,0.4;
71    0.7,0.1,0.4,0.6,0.1,1;
72    0,0.6,0.2,0.3,0.4,0.5;
73    0.6,0.7,0.5,0.7,0.7,0.5;
74    0.3,0.3,0.4,0.5,0.9,0.6]+%i*eye(A);
75 A=inv(X1)*A*X1;
76
77 Ab1=bdiag(A);
78 if or(triu(Ab1)<>Ab1) then pause,end
79 [Ab2,X]=bdiag(A);
80 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
81 if Err(Ab2-inv(X)*A*X )>1.d-8 then pause,end 
82
83 [Ab2,X,bs]=bdiag(A);
84 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
85 if Err(Ab2-inv(X)*A*X )>1.d-8 then pause,end 
86 if size(bs,2)<>1 then pause,end
87 if sum(bs)<>size(A,1) then pause,end
88 if or(bs<=0) then pause,end
89 //Large dimension
90 //---------------
91 //Real case
92 A=rand(25,25);
93 Ab1=bdiag(A);
94 if or(triu(Ab1,-1)<>Ab1) then pause,end
95 [Ab2,X]=bdiag(A);
96 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
97 if Err(Ab2-inv(X)*A*X )>10000*%eps then pause,end 
98
99 [Ab2,X,bs]=bdiag(A);
100 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
101 if Err(Ab2-inv(X)*A*X )>10000*%eps then pause,end 
102 if size(bs,2)<>1 then pause,end
103 if sum(bs)<>size(A,1) then pause,end
104 if or(bs<=0) then pause,end
105
106 //Complex case
107 A=rand(25,25)+%i*rand(25,25);
108 Ab1=bdiag(A);
109 if or(triu(Ab1)<>Ab1) then pause,end
110 [Ab2,X]=bdiag(A);
111 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
112 if Err(Ab2-inv(X)*A*X )>10000*%eps then pause,end 
113
114 [Ab2,X,bs]=bdiag(A);
115 //if Err(Ab2-Ab1)>>10*%eps then pause,end 
116 if Err(Ab2-inv(X)*A*X )>10000*%eps then pause,end 
117 if size(bs,2)<>1 then pause,end
118 if sum(bs)<>size(A,1) then pause,end
119 if or(bs<=0) then pause,end
120
121