home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / LINEARAL.PAK / CHOLESKY.M next >
Encoding:
Text File  |  1992-07-29  |  2.4 KB  |  91 lines
  1.  
  2.  
  3. (*:Version: Mathematica 2.0 *)
  4.  
  5. (*:Name: LinearAlgebra`Cholesky` *)
  6.  
  7. (*:Title: Symmetric Matrix Routines *)
  8.  
  9. (*:Author:
  10.     Jerry B. Keiper (Wolfram Research Inc.)
  11.     Revised by Hon Wah Tam (Wolfram Research Inc.)
  12. *)
  13.  
  14. (*:Keywords:
  15.     symmetric matrices, Cholesky decomposition
  16. *)
  17.  
  18. (*:Requirements:
  19.     The input matrix is assumed to be symmetric positive definite.
  20.     Cholesky.m does not test for symmetry.
  21. *)
  22.  
  23. (*:Warnings: Does not check whether A is symmetric positive definite, or
  24.              even complex. *)
  25.  
  26. (*:Sources:
  27.  
  28. *)
  29.  
  30. (*:Summary:
  31. For a given symmetric positive definite matrix A, CholeskyDecomposition[A]
  32. returns an upper-triangular matrix U such that A = Transpose[U] .
  33. U.  This function uses the row-oriented Cholesky decomposition
  34. method to compute U row by row.  That A is symmetric positive
  35. definite stipulates that A is real.  The answers are totally garbage
  36. if A is complex.  If A is not positive definite, U will contain
  37. Indeterminate entries.
  38. *)
  39.  
  40.  
  41. BeginPackage["LinearAlgebra`Cholesky`"]
  42.  
  43. CholeskyDecomposition::usage =
  44. "For a symmetric positive definite matrix A,
  45. CholeskyDecomposition[A] returns an upper-triangular matrix U
  46. such that A = Transpose[U] . U."
  47.  
  48.  
  49. Begin["LinearAlgebra`Cholesky`private`"]
  50.  
  51.  
  52. column[matrix_List, index_Integer]:= Map[(#[[index]])&, matrix]
  53.  
  54. diagonalprod[m_List, u_List, index_]:= 
  55.     Module[{partial = Take[column[u, index], index-1]},
  56.         Sqrt[m[[index, index]]-partial . partial]]
  57.  
  58. offdiagonalprod[m_List, u_List, partial1_List, index1_, index2_, diag_]:= 
  59.     Module[{ partial2 = Take[column[u, index2], index1-1]},
  60.          (m[[index1, index2]]-partial1 . partial2)/diag]
  61.  
  62. makerow[m_List, u_List, rowindex_]:=
  63.     Module[{part1 = Table[0, {rowindex-1}],
  64.                partial1 = Take[column[u, rowindex], rowindex-1 ],
  65.            diag = diagonalprod[m, u, rowindex], part2, iter},
  66.     part2 = Table[offdiagonalprod[m, u, partial1, rowindex, iter, diag], 
  67.         {iter, rowindex+1, Length[m]}];
  68.     Join[part1, {diag}, part2]]
  69.  
  70. CholeskyDecomposition[m_List] := Module[{len = Length[m], u, iter1},
  71.     u = DiagonalMatrix[Table[0, {len}]];
  72.     For[iter1=1, iter1<= len, ++iter1, 
  73.         u[[iter1]] = makerow[m, u, iter1]];
  74.         u]
  75.  
  76. End[]  (* LinearAlgebra`Cholesky`Private` *)
  77.  
  78. Protect[CholeskyDecomposition]
  79.  
  80. EndPackage[]  (* LinearAlgebra`Cholesky` *)
  81.  
  82. (*:Limitations: 
  83. Does not decompose complex matrices.
  84. *)
  85.  
  86.  
  87. (*:Examples:
  88. CholeskyDecomposition[ Table[ 1/(i+j),{i,1,5},{j,1,5}] ]
  89.  
  90. *)
  91.