Re: Cholesky factorisation [message #18977] |
Tue, 22 February 2000 00:00 |
Pavel Romashkin
Messages: 166 Registered: April 1999
|
Senior Member |
|
|
I can't help hearing some sarcasm regarding IDL in David McClain's postings.
Why is that? After all, accepting computer's numerical precision is a matter of
faith, too, especially in life-critical applications Raytheon is developing. I
find using upper-level language like IDL very convenient, and availability of
"lower-level" tools like pointers is nice. C is great, but so is IDL I think.
If we don't trust functions written by others, we might as well wright
everything ourselves in C...
Cheers,
Pavel
David McClain wrote:
> I use the method shown in Numerical Recipes, 2nd Ed. Coding in C and being
> clever allows it to run in under 3 microsec on a 6x6 array, which is a size
> commonly needed by one of our analyses. Doing it externally allows for both
> speed and error handling as you would want it. It also guarantees
> correctness, as verified by you, the coder; something that RSI would have
> you accept as a matter of faith...
>
> David McClain, Sr. Scientist
> Raytheon Systems Co.
> Tucson, AZ
|
|
|
Re: Cholesky factorisation [message #18989 is a reply to message #18977] |
Tue, 22 February 2000 00:00  |
David McClain
Messages: 17 Registered: January 1999
|
Junior Member |
|
|
I use the method shown in Numerical Recipes, 2nd Ed. Coding in C and being
clever allows it to run in under 3 microsec on a 6x6 array, which is a size
commonly needed by one of our analyses. Doing it externally allows for both
speed and error handling as you would want it. It also guarantees
correctness, as verified by you, the coder; something that RSI would have
you accept as a matter of faith...
David McClain, Sr. Scientist
Raytheon Systems Co.
Tucson, AZ
Peter Scarth <p.scarth@mailbox.uq.edu.au> wrote in message
news:01bf7cda$4593eba0$eb2f14ac@default...
> Howdy all,
> I'm trying to determine if a symmetric matrix is positive-definite. If
this
> sounds like gobbdlygook you might like to stop reading here. If the
> idl2matlab translate-o-matic existed, it would probably need to know about
> this to translate the handy \ (mldivide) operator. There a number of tests
> for positive-definiteness and the one I'm looking at is to attempt
cholesky
> factorisation of the matrix. IDL's CHOLDC procedure halts and returns
> CHOLDC: choldc failed (note that matlab can return a value to flag a
> failure in the decomposition).
> The metho that I an using to get the decomposition looks something like
> this:
>
> --------------
> pro choldc2,a,p
> sa=size(a)
> nd=sa[0] & m=sa[1] & n=sa[2]
>
> ; Some initial tests.....
> ; .....
>
> p=fltarr(n)
> for i=0,n-1 do begin
> for j=i,n-1 do begin
> sum=a(j,i)
> if i gt 0 then sum=sum-total(a(0:i-1,i)*a(0:i-1,j))
> if (i eq j) then begin
> if (sum le 0) then begin
> p=-1
> return
> endif else p(i)=sqrt(sum)
> endif else a(i,j)=sum/p(i)
> endfor
> endfor
> return
> end
> ----------------------
>
> it works ok, returns -1 in p if the method fails but only runs at 1/10 the
> speed of the built in version.
> Is it possible to vectorise this further, or has someone out there in
> cyberland found a more elegant solution to this problem already?
>
> Thanks,
>
>
> Peter Scarth
> Biophysical Remote Sensing Group
> The University of Queensland.
>
|
|
|