octave-bug-tracker
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Octave-bug-tracker] [bug #39014] Wrong determinant for some (large) mat


From: Marco Caliari
Subject: [Octave-bug-tracker] [bug #39014] Wrong determinant for some (large) matrices
Date: Fri, 17 May 2013 13:37:42 +0000
User-agent: Mozilla/5.0 (X11; Linux i686; rv:13.0) Gecko/20100101 Firefox/13.0

URL:
  <http://savannah.gnu.org/bugs/?39014>

                 Summary: Wrong determinant for some (large) matrices
                 Project: GNU Octave
            Submitted by: caliari
            Submitted on: Fri 17 May 2013 01:37:41 PM GMT
                Category: Octave Function
                Severity: 3 - Normal
                Priority: 5 - Normal
              Item Group: None
                  Status: None
             Assigned to: None
         Originator Name: 
        Originator Email: 
             Open/Closed: Open
         Discussion Lock: Any
                 Release: 3.6.2
        Operating System: GNU/Linux

    _______________________________________________________

Details:

Dear all,

as reported in the thread
https://mailman.cae.wisc.edu/pipermail/help-octave/2013-May/058753.html, the
following code


A = rand(2100);
det(A*inv(A))


quite often gives Inf instead of 1. The problem seems to be in the way the
product of the diagonal elements of the R (or U) factor of the matrix is done.
I don't know if prod is sufficiently robust with respect to the product of
numbers very different in magnitude. I see


prod([1e-200*ones(1,20),1e200*ones(1,20)])
ans =  1.00000


but


prod([1e-200*ones(1,30),1e200*ones(1,30)])
ans = 0


Suppose prod is fixed if needed. Isn't then it possible to rewrite det
function with a simple m file basically doing


function ret = det(A)
if ishermitian(A)
   ret = prod(diag(chol(A)));
else
   ret = prod(diag([~,U]=lu(A)));
end


By the way, the actual implementation of prod allows the script above to work
fine on A*inv(A).

Cheers,

Marco




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?39014>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/




reply via email to

[Prev in Thread] Current Thread [Next in Thread]