* Bug 16636 fixed: det(sparse) now actually implemented
[scilab.git] / scilab / modules / linear_algebra / tests / unit_tests / det.tst
1 // =============================================================================
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) ????-2008 - INRIA Michael Baudin
4 // Copyright (C) 2015 - Scilab Enterprises - John Gliksberg
5 // Copyright (C) 2021 - Le Mans Universit√© - Samuel GOUGEON
6 //
7 //  This file is distributed under the same license as the Scilab package.
8 // =============================================================================
9 // <-- CLI SHELL MODE -->
10 // <-- NO CHECK REF -->
11
12 //==========================================================================
13 //==============================   det        ==============================
14 //==========================================================================
15
16 //Small dimension
17 //Real
18 A=[1 1; 1 2];
19 assert_checkalmostequal(det(A), 1);
20 [e,m]=det(A);
21 assert_checkalmostequal(e, 0);
22 assert_checkalmostequal(m, 1);
23 //Complex
24 A=A+%i;
25 assert_checkalmostequal(real(det(A)), 1);
26 assert_checkalmostequal(imag(det(A)), 1);
27 [e,m]=det(A);
28 assert_checkalmostequal(e, 0);
29 assert_checkalmostequal(real(m), 1);
30 assert_checkalmostequal(imag(m), 1);
31 //Sparse
32 A=[1 1; 1 2];
33 A=sparse(A);
34 assert_checkalmostequal(det(A), 1);
35 [e,m]=det(A);
36 assert_checkalmostequal(e, 0);
37 assert_checkalmostequal(m, 1);
38 //Polynomials
39 A=[1+%s 1; 1 2+%s];
40 assert_checkequal(det(A), 1+3*%s+%s*%s);
41 //Rationals
42 A=[1+%s 1/%s; 1 2+%s];
43 assert_checkequal(det(A).num, -1+2*%s+3*%s^2+%s^3);
44 assert_checkequal(det(A).den, %s);
45 //Sparse complex
46 A=[1 1; 1 2];
47 A=A+%i;
48 A=sparse(A);
49 assert_checkalmostequal(real(det(A)), 1);
50 assert_checkalmostequal(imag(det(A)), 1);
51 [e,m]=det(A);
52 assert_checkalmostequal(e, 0);
53 assert_checkalmostequal(real(m), 1);
54 assert_checkalmostequal(imag(m), 1);
55
56 //Large dimension
57 //Real
58 v=rand(1,21);
59 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
60 assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
61 [e,m]=det(A);
62 assert_checktrue(abs(m*(10^e) - prod(v)) < 1D-7);
63 //Complex
64 v=(v+rand(v)*%i)/2;
65 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
66 assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
67 [e,m]=det(A);
68 assert_checktrue(abs(m*(10^e) - prod(v)) < 1D-7);
69 //Sparse
70 v=rand(1,21);
71 v=sparse(v);
72 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
73 assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
74 [e,m]=det(A);
75 assert_checktrue(abs(m*(10^e) - prod(v)) < 1d-7);
76 //Polynomials
77 v=rand(1,21);
78 v=v+%s;
79 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
80 assert_checktrue(abs(coeff(det(A)-prod(v))) < 1D-7);
81 //Rationals
82 v=rand(1,21);
83 v=v/%s;
84 A=rand(21,21); A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
85 assert_checktrue(abs(coeff(det(A).num - prod(v).num)) < 1D-7);
86 //Sparse complex
87 v=rand(1,21);
88 v=(v+rand(v)*%i)/2;
89 v=sparse(v);
90 A=rand(21,21);
91 A=(triu(A,1)+diag(v))*(tril(A,-1)+diag(ones(1,21)));
92 A=sparse(A);
93 assert_checktrue(abs(det(A) - prod(v)) < 1D-7);
94 [e,m]=det(A);
95 assert_checktrue(abs(m*(10^e) - prod(v)) < 1d-7);
96
97 //Error messages
98 A=[1 1; 1 2];
99 errmsg1 = msprintf(_("%s: Wrong type for input argument #%d: Square matrix expected.\n"), "det", 1);
100 assert_checkerror("det([1,2;3,4;5,6])", errmsg1, 20);
101 errmsg2 = msprintf(_("%s: Wrong number of input argument(s): %d expected.\n"), "det", 1);
102 assert_checkerror("det(A,1)", errmsg2, 77);
103
104 // Check det == 0 for simple cases
105 A = 0;
106 assert_checkalmostequal(det(A), 0);
107 A = [1 2
108 1 2];
109 assert_checkalmostequal(det(A), 0);
110 A = [1 2 3
111 1 2 3
112 1 2 3];
113 assert_checkalmostequal(det(A), 0);
114 b = rand(1, 5);
115 A = [b
116 b
117 b
118 b
119 b];
120 assert_checkalmostequal(det(A), 0, %eps, norm(A)*%eps);
121
122 a = sprand(1000,1000,0.1); a(:,1:2) = 1;
123 assert_checkequal(det(a), 0);
124
125 // Underflow and overflow management for sparse matrices. [e, m] = det(s)  syntax
126 // -------------------------------------------------------------------------------
127 // http://bugzilla.scilab.org/16636 :
128 // With a real matrix
129 // ..................
130 // Underflow:
131 n = 75;
132 while n < 3000
133     s = sparse(triu(rand(n,n)));
134     ref = prod(diag(s));
135     if n < 1000
136         assert_checkalmostequal(det(s), ref, 10*n*%eps);
137     else
138         assert_checkequal(det(s), 0);
139     end
140     [e, m] = det(s);
141     ref = sum(log10(full(diag(s))));
142     assert_checkequal(e, int(ref)-1);
143     assert_checkalmostequal(m, 10^(ref-int(ref)+1), 10*n*%eps);
144     n = n*2;
145 end
146 // Overflow:
147 n = 75;
148 while n < 3000
149     s = sparse(triu(rand(n,n)))*6;
150     ref = prod(diag(s));
151     if n < 1000
152         assert_checkalmostequal(det(s), ref, 10*n*%eps);
153     else
154         assert_checkequal(det(s), %inf);
155     end
156     [e, m] = det(s);
157     ref = sum(log10(full(diag(s))));
158     assert_checkequal(e, int(ref));
159     assert_checkalmostequal(m, 10^(ref-int(ref)), 10*n*%eps);
160     n = n*2;
161 end
162 // With a complex matrix
163 // .....................
164 // Underflow:
165 n = 75;
166 while n < 3000
167     s = sparse(triu(complex(rand(n,n),rand(n,n))));
168     ref = prod(diag(s));
169     if n < 2000
170         assert_checkalmostequal(det(s), ref, 100*n*%eps);
171     else
172         assert_checkequal(det(s), 0*%i);
173     end
174     [e, m] = det(s);
175     [eref, mref] = det(full(s));
176     assert_checkequal(e, eref);
177     assert_checkalmostequal(m, mref, 100*n*%eps);
178     n = n*2;
179 end
180 // Overflow:
181 n = 75;
182 while n < 3000
183     s = sparse(triu(complex(rand(n,n),rand(n,n))))*6;
184     ref = prod(diag(s));
185     if n < 600
186         assert_checkalmostequal(det(s), ref, 100*n*%eps);
187     else
188         assert_checkequal(abs(real(det(s))), %inf);
189     end
190     [e, m] = det(s);
191     [eref, mref] = det(full(s));
192     assert_checkequal(e, eref);
193     assert_checkalmostequal(m, mref, 200*n*%eps);
194     n = n*2;
195 end