Re: Cholesky factorisation

["David McClain" <dmcclain(at)azstarnet.com>]

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
01bf7cda$4593eba0$eb2f14ac@default">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.
>