/* CHOLESKY.G: procedure to compute Cholesky factorization of an
	       NxN symmetric postive semidefinite matrix A
	       Procedure returns lower triangular Cholesky factor of A

	       John Rust, Yale University, February, 1997 */

   proc (1)=cholesky(a);

   local i,j,k,n,t;

   n=rows(a);

   /* check to make sure input matrix is square */

   if (n /= cols(a)); "Error: input matrix must be square"; retp(-1); endif;

   /* check to make sure input matrix is symmetric */

   k=1; do until k > n;
       i=k+1; do until i > n;
	  if (a[k,i] /= a[i,k]); 
	     "Error: input matrix not symmetric";
	      retp(-1); 
	  endif;
       i=i+1; endo;	     
   k=k+1; endo;

   /* OK, now compute Cholesky factor, storing it back into a matrix */

   k=1; do until k > n;

      i=1; do until i > k-1;

	  if (i==1);
	     a[k,i]=a[k,i]/a[1,1];
          else;
	     a[k,i]=(a[k,i]-sumc(a[i,1:i-1].*a[k,1:i-1]))/a[i,i];
          endif;
      i=i+1; endo;

     if (k==1);
	 t=a[1,1];
     else;
         t=a[k,k]-sumc(a[k,1:k-1]^2);
     endif; 

     if (t == 0); 
	 a[k,k]=0;
     elseif (t < 0); /* can't take square root of a negative number! */
	 "Error: input matrix not positive semidefinite"; retp(-1); 
     else;
         a[k,k]=sqrt(t);
     endif;
 
   k=k+1; endo;

   /* finally zero out the upper triangle of the input matrix */

   k=1; do until k > n;
     i=k+1; do until i > n;
	a[k,i]=0;
     i=i+1; endo;
   k=k+1; endo;

   retp(a);

   endp;