* Bug 16636 fixed: det(sparse) now actually implemented
[scilab.git] / scilab / modules / overloading / macros / %sp_det.sci
1 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 // Copyright (C) 2013 - Scilab Enterprises - Charlotte HECQUET
3 // Copyright (C) 2012 - 2016 - Scilab Enterprises
4 // Copyright (C) 2021 - Samuel GOUGEON - Le Mans Universit√©
5 //
6 // This file is hereby licensed under the terms of the GNU GPL v2.0,
7 // pursuant to article 5.3.4 of the CeCILL v.2.1.
8 // This file was originally licensed under the terms of the CeCILL v2.1,
9 // and continues to be available under such terms.
10 // For more information, see the COPYING file which you should have received
11 // along with this program.
12
13 function [res1, res2] = %sp_det(A)
14     [lhs, rhs] = argn(0);
15     if length(A)==0 then
16         [res1, res2] = (1,1);
17         if lhs>1
18             res1 = 0
19         end
20         return
21     end
22     if nnz(A)==0 then
23         [res1, res2] = (0,0)
24         return
25     end
26
27     wstatus = warning("query"); warning("off")  // to avoid the "is singular" message
28     hand = umf_lufact(A);          // umfpack is used for complex sparse matrix
29     [L,U,P,Q,r] = umf_luget(hand);
30     umf_ludel(hand);
31     warning(wstatus)
32
33     dU = clean(diag(U), 0, %eps);
34     r = clean(r, 0, %eps);
35     if isreal(A) then
36         p = r .* dU;
37         tmp = sum(log10(full(abs(p))))
38         res1 = int(tmp)
39         if res1 < 0
40             res1 = res1 - 1
41         end
42         res2 = 10^(tmp-res1) * prod(sign(p))
43     else
44         tmp = sum(log([r full(dU)]))
45         phase = imag(tmp)
46         e10 = real(tmp)/log(10) + log10(abs(cos(phase)))
47         res1 = int(e10)
48         res2 = 10^(real(tmp)/log(10)-res1) * exp(%i*phase)
49         while abs(res2) >= 10
50             res2 = res2 / 10
51             res1 = res1 + 1
52         end
53         while abs(res2)<1 & abs(res2)<>0
54             res2 = res2 * 10
55             res1 = res1 - 1
56         end
57     end
58     if res1 == -%inf
59         [res1, res2] = (0,0)
60     end
61     if lhs == 1 then
62         res1 = res2 * 10^res1
63     end
64 endfunction