[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/
- [Octave-bug-tracker] [bug #39014] Wrong determinant for some (large) matrices,
Marco Caliari <=