e323097ff3d033fc55b2102ca99ad855d72aa78a
[scilab.git] / scilab / modules / linear_algebra / tests / unit_tests / svd.dia.ref
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) ????-2008 - INRIA
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 //==========================================================================
12 //==============================     svd      ==============================
13 //==========================================================================
14 //Empty matrix
15 A=[];
16 [U,S,V]=svd(A);
17 if U<>[]|V<>[]|S<>[] then bugmes();quit;end
18 S=svd(A);
19 if S<>[] then bugmes();quit;end
20 [U,S,V]=svd(A,"e");
21 if U<>[]|V<>[]|S<>[] then bugmes();quit;end
22 S=svd(A,"e");
23 if S<>[] then bugmes();quit;end
24 //Matrix with inf or nan
25 if execstr('svd([%inf 1;2 3])','errcatch')==0 then bugmes();quit;end
26 if execstr('svd([1 %nan;2 3])','errcatch')==0 then bugmes();quit;end
27 if execstr('svd([%inf %i;2 3])','errcatch')==0 then bugmes();quit;end
28 if execstr('svd([%i %i;%nan 3])','errcatch')==0 then bugmes();quit;end
29 //Small dimension
30 //---------------
31 A=rand(3,5);
32 Ac=A+%i*rand(A);
33 //Real Case
34 [U,S,V]=svd(A);
35 if Err(U*S*V'-A)>200*%eps then bugmes();quit;end
36 if Err(svd(A)-diag(S))> 200*%eps then bugmes();quit;end
37 [U,S,V]=svd(A,"e");
38 if Err(U*S*V'-A)>200*%eps then bugmes();quit;end
39 A=A';
40 [U,S,V]=svd(A);
41 if Err(U*S*V'-A)>200*%eps then bugmes();quit;end
42 if Err(svd(A)-diag(S))> 200*%eps then bugmes();quit;end
43 [U,S,V]=svd(A,"e");
44 if Err(U*S*V'-A)>200*%eps then bugmes();quit;end
45 //Complex Case
46 [U,S,V]=svd(Ac);
47 if Err(U*S*V'-Ac)>200*%eps then bugmes();quit;end
48 if Err(svd(Ac)-diag(S))> 200*%eps then bugmes();quit;end
49 [U,S,V]=svd(Ac,"e");
50 if Err(U*S*V'-Ac)>200*%eps then bugmes();quit;end
51 Ac=Ac';
52 [U,S,V]=svd(Ac);U*S*V'-Ac;
53 if Err(U*S*V'-Ac)>200*%eps then bugmes();quit;end
54 if Err(svd(Ac)-diag(S))> 200*%eps then bugmes();quit;end
55 [U,S,V]=svd(Ac,"e");
56 if Err(U*S*V'-Ac)>200*%eps then bugmes();quit;end
57 //Large dimension
58 //---------------
59 A=rand(150,60);
60 Ac=A+rand(A)*%i;
61 //Real Case
62 [U,S,V]=svd(A);
63 if Err(U*S*V'-A)>10000*%eps then bugmes();quit;end
64 if Err(svd(A)-diag(S))> 10000*%eps then bugmes();quit;end
65 [U,S,V]=svd(A,"e");
66 if Err(U*S*V'-A)>10000*%eps then bugmes();quit;end
67 A=A';
68 [U,S,V]=svd(A);
69 if Err(U*S*V'-A)>10000*%eps then bugmes();quit;end
70 if Err(svd(A)-diag(S))> 10000*%eps then bugmes();quit;end
71 [U,S,V]=svd(A,"e");
72 if Err(U*S*V'-A)>10000*%eps then bugmes();quit;end
73 //Complex Case
74 [U,S,V]=svd(Ac);
75 if Err(U*S*V'-Ac)>10000*%eps then bugmes();quit;end
76 if Err(svd(Ac)-diag(S))> 10000*%eps then bugmes();quit;end
77 [U,S,V]=svd(Ac,"e");
78 if Err(U*S*V'-Ac)>10000*%eps then bugmes();quit;end
79 Ac=Ac';
80 [U,S,V]=svd(Ac);U*S*V'-Ac;
81 if Err(U*S*V'-Ac)>10000*%eps then bugmes();quit;end
82 if Err(svd(Ac)-diag(S))> 10000*%eps then bugmes();quit;end
83 [U,S,V]=svd(Ac,"e");
84 if Err(U*S*V'-Ac)>10000*%eps then bugmes();quit;end
85 //==========================================================================
86 //==============================     svd part II     =======================
87 //==========================================================================
88 //Empty matrix
89 if svd([])<>[] then bugmes();quit;end
90 if svd([],"e")<>[] then bugmes();quit;end
91 [U,S]=svd([]);
92 if U<>[]|S<>[]  then bugmes();quit;end
93 [U,S,V]=svd([]);
94 if U<>[]|S<>[]|V<>[]  then bugmes();quit;end
95 [U,S,V,rk]=svd([]);
96 if U<>[]|S<>[]|V<>[]|rk<>0  then bugmes();quit;end
97 [U,S,V,rk]=svd([],%eps);
98 if U<>[]|S<>[]|V<>[]|rk<>0  then bugmes();quit;end
99 if execstr('[U,S,V,rk]=svd([],'"e'")','errcatch') == 0 then bugmes();quit;end
100 //Small dimension
101 //Real Case Fat
102 A=rand(3,5);
103 S=svd(A);
104 if or(S<0) then bugmes();quit;end
105 if gsort(S)<>S  then bugmes();quit;end
106 [U,S1]=svd(A);
107 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
108 if Err(U'*U-eye())>10*%eps  then bugmes();quit;end
109 [U1,S1]=svd(A,"e");
110 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
111 if Err(U1'*U1-eye())>200*%eps  then bugmes();quit;end
112 [U1,S1,V]=svd(A);
113 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
114 if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
115 if Err(U1*S1*V'-A) >200*%eps  then bugmes();quit;end
116 [U1,S1,V1]=svd(A,"e");
117 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
118 if Err(U-U1)>10*%eps  then bugmes();quit;end
119 if Err(U1*S1*V1'-A) >200*%eps  then bugmes();quit;end
120 [U1,S1,V1,rk]=svd(A);
121 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
122 if Err(U-U1)>200*%eps  then bugmes();quit;end
123 if Err(V-V1) >200*%eps  then bugmes();quit;end
124 if rk<>3 then bugmes();quit;end
125 //Real Case Tall
126 A=rand(5,3);
127 S=svd(A);
128 if or(S<0) then bugmes();quit;end
129 if gsort(S)<>S  then bugmes();quit;end
130 [U,S1]=svd(A);
131 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
132 if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
133 [U1,S1]=svd(A,"e");
134 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
135 if Err(U1'*U1-eye())>200*%eps  then bugmes();quit;end
136 [U1,S1,V]=svd(A);
137 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
138 if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
139 if Err(U1*S1*V'-A) >200*%eps  then bugmes();quit;end
140 [U1,S1,V1]=svd(A,"e");
141 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
142 if size(U1,2)<>3 then bugmes();quit;end
143 if Err(U1*S1*V1'-A) >200*%eps  then bugmes();quit;end
144 [U1,S1,V1,rk]=svd(A);
145 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
146 if Err(U-U1)>200*%eps  then bugmes();quit;end
147 if Err(V-V1) >200*%eps  then bugmes();quit;end
148 if rk<>3 then bugmes();quit;end
149 //Complex Case Fat
150 A=rand(3,5)+%i*rand(3,5);
151 S=svd(A);
152 if or(S<0) then bugmes();quit;end
153 if gsort(S)<>S  then bugmes();quit;end
154 [U,S1]=svd(A);
155 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
156 if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
157 [U1,S1]=svd(A,"e");
158 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
159 if Err(U1'*U1-eye())>200*%eps  then bugmes();quit;end
160 [U1,S1,V]=svd(A);
161 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
162 if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
163 if Err(U1*S1*V'-A) >30*%eps  then bugmes();quit;end
164 [U1,S1,V1]=svd(A,"e");
165 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
166 if Err(U-U1)>200*%eps  then bugmes();quit;end
167 if Err(U1*S1*V1'-A) >200*%eps  then bugmes();quit;end
168 [U1,S1,V1,rk]=svd(A);
169 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
170 if Err(U-U1)>200*%eps  then bugmes();quit;end
171 if Err(V-V1) >200*%eps  then bugmes();quit;end
172 if rk<>3 then bugmes();quit;end
173 //Complex Case Tall
174 A=rand(5,3)+%i*rand(5,3);
175 S=svd(A);
176 if or(S<0) then bugmes();quit;end
177 if gsort(S)<>S  then bugmes();quit;end
178 [U,S1]=svd(A);
179 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
180 if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
181 [U1,S1]=svd(A,"e");
182 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
183 if Err(U1'*U1-eye())>200*%eps  then bugmes();quit;end
184 [U1,S1,V]=svd(A);
185 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
186 if Err(U'*U-eye())>200*%eps  then bugmes();quit;end
187 if Err(U1*S1*V'-A) >200*%eps  then bugmes();quit;end
188 [U1,S1,V1]=svd(A,"e");
189 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
190 if size(U1,2)<>3 then bugmes();quit;end
191 if Err(U1*S1*V1'-A) >200*%eps  then bugmes();quit;end
192 [U1,S1,V1,rk]=svd(A);
193 if Err(S-diag(S1))>200*%eps  then bugmes();quit;end
194 if Err(U-U1)>200*%eps  then bugmes();quit;end
195 if Err(V-V1) >200*%eps  then bugmes();quit;end
196 if rk<>3 then bugmes();quit;end
197 //Large dimension
198 //Real Case Fat
199 A=rand(30,50);
200 S=svd(A);
201 if or(S<0) then bugmes();quit;end
202 if gsort(S)<>S  then bugmes();quit;end
203 [U,S1]=svd(A);
204 if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
205 if Err(U'*U-eye())>1000*%eps  then bugmes();quit;end
206 [U1,S1]=svd(A,"e");
207 if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
208 if Err(U1'*U1-eye())>1000*%eps  then bugmes();quit;end
209 [U1,S1,V]=svd(A);
210 if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
211 if Err(U'*U-eye())>1000*%eps  then bugmes();quit;end
212 if Err(U1*S1*V'-A) >1000*%eps  then bugmes();quit;end
213 [U1,S1,V1]=svd(A,"e");
214 if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
215 if Err(U-U1)>10*%eps  then bugmes();quit;end
216 if Err(U1*S1*V1'-A) >1000*%eps  then bugmes();quit;end
217 [U1,S1,V1,rk]=svd(A);
218 if Err(S-diag(S1))>1000*%eps  then bugmes();quit;end
219 if Err(U-U1)>5000*%eps  then bugmes();quit;end
220 if Err(V-V1) >5000*%eps  then bugmes();quit;end
221 if rk<>30 then bugmes();quit;end
222 //Real Case Tall
223 A=rand(50,30);
224 S=svd(A);
225 if or(S<0) then bugmes();quit;end
226 if gsort(S)<>S  then bugmes();quit;end
227 [U,S1]=svd(A);
228 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
229 if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
230 [U1,S1]=svd(A,"e");
231 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
232 if Err(U1'*U1-eye())>5000*%eps  then bugmes();quit;end
233 [U1,S1,V]=svd(A);
234 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
235 if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
236 if Err(U1*S1*V'-A) >5000*%eps  then bugmes();quit;end
237 [U1,S1,V1]=svd(A,"e");
238 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
239 if size(U1,2)<>30 then bugmes();quit;end
240 if Err(U1*S1*V1'-A) >5000*%eps  then bugmes();quit;end
241 [U1,S1,V1,rk]=svd(A);
242 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
243 if Err(U-U1)>5000*%eps  then bugmes();quit;end
244 if Err(V-V1) >5000*%eps  then bugmes();quit;end
245 if rk<>30 then bugmes();quit;end
246 //Complex Case Fat
247 A=rand(30,50)+%i*rand(30,50);
248 S=svd(A);
249 if or(S<0) then bugmes();quit;end
250 if gsort(S)<>S  then bugmes();quit;end
251 [U,S1]=svd(A);
252 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
253 if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
254 [U1,S1]=svd(A,"e");
255 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
256 if Err(U1'*U1-eye())>5000*%eps  then bugmes();quit;end
257 [U1,S1,V]=svd(A);
258 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
259 if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
260 if Err(U1*S1*V'-A) >5000*%eps  then bugmes();quit;end
261 [U1,S1,V1]=svd(A,"e");
262 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
263 if Err(U-U1)>5000*%eps  then bugmes();quit;end
264 if Err(U1*S1*V1'-A) >5000*%eps  then bugmes();quit;end
265 [U1,S1,V1,rk]=svd(A);
266 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
267 if Err(U-U1)>5000*%eps  then bugmes();quit;end
268 if Err(V-V1) >5000*%eps  then bugmes();quit;end
269 if rk<>30 then bugmes();quit;end
270 //Complex Case Tall
271 A=rand(50,30)+%i*rand(50,30);
272 S=svd(A);
273 if or(S<0) then bugmes();quit;end
274 if gsort(S)<>S  then bugmes();quit;end
275 [U,S1]=svd(A);
276 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
277 if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
278 [U1,S1]=svd(A,"e");
279 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
280 if Err(U1'*U1-eye())>5000*%eps  then bugmes();quit;end
281 [U1,S1,V]=svd(A);
282 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
283 if Err(U'*U-eye())>5000*%eps  then bugmes();quit;end
284 if Err(U1*S1*V'-A) >5000*%eps  then bugmes();quit;end
285 [U1,S1,V1]=svd(A,"e");
286 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
287 if size(U1,2)<>30 then bugmes();quit;end
288 if Err(U1*S1*V1'-A) >5000*%eps  then bugmes();quit;end
289 [U1,S1,V1,rk]=svd(A);
290 if Err(S-diag(S1))>5000*%eps  then bugmes();quit;end
291 if Err(U-U1)>5000*%eps  then bugmes();quit;end
292 if Err(V-V1) >5000*%eps  then bugmes();quit;end
293 if rk<>30 then bugmes();quit;end