From: spellucci@mathematik.tu-darmstadt.de (Peter Spellucci)
Subject: Re: Matrix inversion
Date: 1 Oct 1999 13:24:50 GMT
Newsgroups: sci.math.num-analysis
Keywords: basic code for Cholesky decomposition

In article <omaeq49oah.fsf@prancer.cs.utk.edu>,
 Victor Eijkhout <eijkhout@prancer.cs.utk.edu> writes:
|> lg@kt.dtu.dk (Lars Gregersen) writes:
|> 
|> > Cholesky decomposition do not pivoting which also make the bookkeeping
|> > code much simpler.
|> 
|> Yeah, but try implementing Cholesy with a matrix where only half is stored.
|> It's a headache.
|> 
|> -- 
|> Victor Eijkhout
why? 
write down the ordinary cholesky-decomposition, e.g. 
      ifail=0
      do    i=1,n
        h=a(i,i)
        do   j=1,i-1
          h=h-l(i,j)**2
        enddo
        if ( h .le. 0.d0 ) then
          ifail=1
          return
        endif
        h=dsqrt(h)
        l(i,i)=h
        h=1.d0/h
        do    k=i+1,n
          s=0.d0
          do    j=1,i-1
            s=s+l(i,j)*l(k,j)
          enddo
          s=(a(k,i)-s)*h
          l(k,i)=s
        enddo
      enddo
      return
you may even identify a and l here. now translate any index pair 

    (k,j)  into (k*(k-1)/2)+j  

done

hope this helps
peter