(*:Version: Mathematica 2.0 *)

(*:Name: LinearAlgebra`Cholesky` *)

(*:Title: Symmetric Matrix Routines *)

(*:Author:
	Jerry B. Keiper (Wolfram Research Inc.)
	Revised by Hon Wah Tam (Wolfram Research Inc.)
*)

(*:Keywords:
	symmetric matrices, Cholesky decomposition
*)

(*:Requirements:
	The input matrix is assumed to be symmetric positive definite.
	Cholesky.m does not test for symmetry.
*)

(*:Warnings: Does not check whether A is symmetric positive definite, or
             even complex. *)

(*:Sources:

*)

(*:Summary:
For a given symmetric positive definite matrix A, CholeskyDecomposition[A]
returns an upper-triangular matrix U such that A = Transpose[U] .
U.  This function uses the row-oriented Cholesky decomposition
method to compute U row by row.  That A is symmetric positive
definite stipulates that A is real.  The answers are totally garbage
if A is complex.  If A is not positive definite, U will contain
Indeterminate entries.
*)


BeginPackage["LinearAlgebra`Cholesky`"]

CholeskyDecomposition::usage =
"For a symmetric positive definite matrix A,
CholeskyDecomposition[A] returns an upper-triangular matrix U
such that A = Transpose[U] . U."


Begin["LinearAlgebra`Cholesky`private`"]


column[matrix_List, index_Integer]:= Map[(#[[index]])&, matrix]

diagonalprod[m_List, u_List, index_]:= 
	Module[{partial = Take[column[u, index], index-1]},
		Sqrt[m[[index, index]]-partial . partial]]

offdiagonalprod[m_List, u_List, partial1_List, index1_, index2_, diag_]:= 
	Module[{ partial2 = Take[column[u, index2], index1-1]},
	     (m[[index1, index2]]-partial1 . partial2)/diag]

makerow[m_List, u_List, rowindex_]:=
	Module[{part1 = Table[0, {rowindex-1}],
               partial1 = Take[column[u, rowindex], rowindex-1 ],
	       diag = diagonalprod[m, u, rowindex], part2, iter},
	part2 = Table[offdiagonalprod[m, u, partial1, rowindex, iter, diag], 
		{iter, rowindex+1, Length[m]}];
	Join[part1, {diag}, part2]]

CholeskyDecomposition[m_List] := Module[{len = Length[m], u, iter1},
	u = DiagonalMatrix[Table[0, {len}]];
	For[iter1=1, iter1<= len, ++iter1, 
		u[[iter1]] = makerow[m, u, iter1]];
		u]

End[]  (* LinearAlgebra`Cholesky`Private` *)

Protect[CholeskyDecomposition]

EndPackage[]  (* LinearAlgebra`Cholesky` *)

(*:Limitations: 
Does not decompose complex matrices.
*)


(*:Examples:
CholeskyDecomposition[ Table[ 1/(i+j),{i,1,5},{j,1,5}] ]

*)