1//<-- CLI SHELL MODE -->
2// =============================================================================
3// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4// Copyright (C) ????-2008 - INRIA
5//
6//  This file is distributed under the same license as the Scilab package.
7// =============================================================================
8function r=Err(x),r=norm(x,1),endfunction
9rand('normal')
10
11//==========================================================================
12//==============================      lu      ==============================
13//==========================================================================
14//Empty matrix
15A=[];
16[L,U]=lu(A);
17if L<>[]|U<>[] then pause,end
18[L,U,E]=lu(A);
19if L<>[]|U<>[]|E<>[] then pause,end
20//Non full rank
21A=rand(5,2);A=A*A';Ac=rand(5,2)+%i*rand(5,2);Ac=Ac*Ac';
22[L,U,E]=lu(A);
23if Err(L*U-E*A) >200*%eps then pause,end
24[L,U,E]=lu(Ac);
25if Err(L*U-E*Ac) >200*%eps then pause,end
26
27//Small dimension
28//---------------
29//Square
30A=rand(5,5);Ac=A+%i*rand(A);
31//Real case
32[L,U]=lu(A);
33if Err(L*U-A) >200*%eps then pause,end
34[L,U,E]=lu(A);
35if Err(L*U-E*A) >200*%eps then pause,end
36//Complex case
37[L,U]=lu(Ac);
38if Err(L*U-Ac) >200*%eps then pause,end
39[L,U,E]=lu(Ac);
40if Err(L*U-E*Ac) >200*%eps then pause,end
41//Fat
42A=rand(3,5);Ac=A+%i*rand(A);
43//Real case
44[L,U]=lu(A);
45if Err(L*U-A) >200*%eps then pause,end
46[L,U,E]=lu(A);
47if Err(L*U-E*A) >200*%eps then pause,end
48//Complex case
49[L,U]=lu(Ac);
50if Err(L*U-Ac) >200*%eps then pause,end
51[L,U,E]=lu(Ac);
52if Err(L*U-E*Ac) >200*%eps then pause,end
53//Tall
54A=rand(5,3);Ac=A+%i*rand(A);
55//Real case
56[L,U]=lu(A);
57if Err(L*U-A) >200*%eps then pause,end
58[L,U,E]=lu(A);
59if Err(L*U-E*A) >200*%eps then pause,end
60//Complex case
61[L,U]=lu(Ac);
62if Err(L*U-Ac) >200*%eps then pause,end
63[L,U,E]=lu(Ac);
64if Err(L*U-E*Ac) >200*%eps then pause,end
65
66//large dimension
67//---------------
68//Square
69A=rand(50,50);Ac=A+%i*rand(A);
70//Real case
71[L,U]=lu(A);
72if Err(L*U-A) >1000*%eps then pause,end
73[L,U,E]=lu(A);
74if Err(L*U-E*A) >1000*%eps then pause,end
75//Complex case
76[L,U]=lu(Ac);
77if Err(L*U-Ac) >1000*%eps then pause,end
78[L,U,E]=lu(Ac);
79if Err(L*U-E*Ac) >1000*%eps then pause,end
80//Fat
81A=rand(30,50);Ac=A+%i*rand(A);
82//Real case
83[L,U]=lu(A);
84if Err(L*U-A) >1000*%eps then pause,end
85[L,U,E]=lu(A);
86if Err(L*U-E*A) >1000*%eps then pause,end
87//Complex case
88[L,U]=lu(Ac);
89if Err(L*U-Ac) >1000*%eps then pause,end
90[L,U,E]=lu(Ac);
91if Err(L*U-E*Ac) >1000*%eps then pause,end
92//Tall
93A=rand(50,30);Ac=A+%i*rand(A);
94//Real case
95[L,U]=lu(A);
96if Err(L*U-A) >1000*%eps then pause,end
97[L,U,E]=lu(A);
98if Err(L*U-E*A) >1000*%eps then pause,end
99//Complex case
100[L,U]=lu(Ac);
101if Err(L*U-Ac) >1000*%eps then pause,end
102[L,U,E]=lu(Ac);
103if Err(L*U-E*Ac) >1000*%eps then pause,end
104
105
106
107