6e4cde9f8d8e8649c8b23a97c02291312de7599c
[scilab.git] / scilab / modules / linear_algebra / tests / unit_tests / bdiag.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),r=norm(x,1),endfunction
15
16 rand('normal')
17
18 //==========================================================================
19
20 //==============================   bdiag      ============================== 
21
22 //==========================================================================
23
24 if bdiag([])<>[] then bugmes();quit;end
25
26 [ab,x]=bdiag([]);
27
28 if ab<>[]|x<>[] then bugmes();quit;end
29
30 [ab,x,bs]=bdiag([]);
31
32 if ab<>[]|x<>[]|bs<>[] then bugmes();quit;end
33
34 if execstr('bdiag([1 2;3 4;5 6])','errcatch')==0 then bugmes();quit;end
35
36 //Small dimension
37
38 //---------------
39
40 //Real case
41
42 e=1.d-1;
43
44 A=[1 1  2 3 4 5
45    0 1  6 7 8 9
46    0 0  1 e 3 1
47    0 0 -e 1 5 9
48    0 0  0 0 2 e
49    0 0  0 0 0 3];
50
51 X1=[0.5,0.3,0,0.3,0.3,0.2;
52    1,0.6,0.5,0.1,0.7,0.4;
53    0.7,0.1,0.4,0.6,0.1,1;
54    0,0.6,0.2,0.3,0.4,0.5;
55    0.6,0.7,0.5,0.7,0.7,0.5;
56    0.3,0.3,0.4,0.5,0.9,0.6]X1  = 
57
58   0.5    0.3    0      0.3    0.3    0.2  
59   1      0.6    0.5    0.1    0.7    0.4  
60   0.7    0.1    0.4    0.6    0.1    1    
61   0      0.6    0.2    0.3    0.4    0.5  
62   0.6    0.7    0.5    0.7    0.7    0.5  
63   0.3    0.3    0.4    0.5    0.9    0.6  
64
65 A=inv(X1)*A*X1;
66
67 Ab1=bdiag(A);
68
69 if or(triu(Ab1,-1)<>Ab1) then bugmes();quit;end
70
71 [Ab2,X]=bdiag(A);
72
73 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
74
75 if Err(Ab2-inv(X)*A*X )>1d6*%eps then bugmes();quit;end 
76
77 [Ab2,X,bs]=bdiag(A);
78
79 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
80
81 if Err(Ab2-inv(X)*A*X )>2.d-10 then bugmes();quit;end 
82
83 if or(size(bs)<>[3,1]) then bugmes();quit;end
84
85 if sum(bs)<>size(A,1) then bugmes();quit;end
86
87 if or(bs<=0) then bugmes();quit;end
88
89 [Ab2,X,bs]=bdiag(A,1);
90
91 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
92
93 if Err(Ab2-inv(X)*A*X )>2d-7 then bugmes();quit;end 
94
95 if or(size(bs)<>[1,1]) then bugmes();quit;end
96
97 if sum(bs)<>size(A,1) then bugmes();quit;end
98
99 if or(bs<=0) then bugmes();quit;end
100
101 //Complex case
102
103 e=1.d-1;
104
105 A=[1 1  2 3 4 5
106    0 1  6 7 8 9
107    0 0  1 e 3 1
108    0 0 -e 1 5 9
109    0 0  0 0 2 e
110    0 0  0 0 0 3];
111
112 X1=[0.5,0.3,0,0.3,0.3,0.2;
113    1,0.6,0.5,0.1,0.7,0.4;
114    0.7,0.1,0.4,0.6,0.1,1;
115    0,0.6,0.2,0.3,0.4,0.5;
116    0.6,0.7,0.5,0.7,0.7,0.5;
117    0.3,0.3,0.4,0.5,0.9,0.6]+%i*eye(A);
118
119 A=inv(X1)*A*X1;
120
121 Ab1=bdiag(A);
122
123 if or(triu(Ab1)<>Ab1) then bugmes();quit;end
124
125 [Ab2,X]=bdiag(A);
126
127 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
128
129 if Err(Ab2-inv(X)*A*X )>1.d-8 then bugmes();quit;end 
130
131 [Ab2,X,bs]=bdiag(A);
132
133 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
134
135 if Err(Ab2-inv(X)*A*X )>1.d-8 then bugmes();quit;end 
136
137 if size(bs,2)<>1 then bugmes();quit;end
138
139 if sum(bs)<>size(A,1) then bugmes();quit;end
140
141 if or(bs<=0) then bugmes();quit;end
142
143 //Large dimension
144
145 //---------------
146
147 //Real case
148
149 A=rand(25,25);
150
151 Ab1=bdiag(A);
152
153 if or(triu(Ab1,-1)<>Ab1) then bugmes();quit;end
154
155 [Ab2,X]=bdiag(A);
156
157 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
158
159 if Err(Ab2-inv(X)*A*X )>10000*%eps then bugmes();quit;end 
160
161 [Ab2,X,bs]=bdiag(A);
162
163 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
164
165 if Err(Ab2-inv(X)*A*X )>10000*%eps then bugmes();quit;end 
166
167 if size(bs,2)<>1 then bugmes();quit;end
168
169 if sum(bs)<>size(A,1) then bugmes();quit;end
170
171 if or(bs<=0) then bugmes();quit;end
172
173 //Complex case
174
175 A=rand(25,25)+%i*rand(25,25);
176
177 Ab1=bdiag(A);
178
179 if or(triu(Ab1)<>Ab1) then bugmes();quit;end
180
181 [Ab2,X]=bdiag(A);
182
183 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
184
185 if Err(Ab2-inv(X)*A*X )>10000*%eps then bugmes();quit;end 
186
187 [Ab2,X,bs]=bdiag(A);
188
189 //if Err(Ab2-Ab1)>>10*%eps then bugmes();quit;end 
190
191 if Err(Ab2-inv(X)*A*X )>10000*%eps then bugmes();quit;end 
192
193 if size(bs,2)<>1 then bugmes();quit;end
194
195 if sum(bs)<>size(A,1) then bugmes();quit;end
196
197 if or(bs<=0) then bugmes();quit;end
198