(* CHOLESKI'S ALGORITHM 6.6
*
*  To factor the positive definite n by n matrix A into LL^T, 
*   where L is the lower triangular. 
*
*   Input:  dimension n; the entries A(i,J), 1<=i, J<=n, of A.
*
*  Output:  the entries L(i,J), 1<=J<=i, 1<=i<=n of L.
*)
Print["\n"];
Print["A(1,1), A(2,1), ..., A(1,n+1), A(2,1), A(2,2), ...\n"];
Print["A(2,n+1), ..., A(n,1), A(n,2), ..., A(n,n+1)\n"];
Print["\n"];
Print["Place as many entries as desired on each line, but \n"];
Print["separate entries with at least one blank\n"];
Print["\n"];
Print["\n"];
AA = InputString["This is the Choleski Factorization Method\n
     The array will be input from a text file in the order:(see screen)\n
     Has the input file been created?\n
     Enter 'yes' or 'no'\n"];
If[AA == "yes" || AA == "y" || AA == "Y",
   NAME = InputString["Input the file name in the form - \n
      drive:\ name.ext      for example\n
         A:\\DATA.DTA\n"];
   INP = OpenRead[NAME];
   OK = 0;
   While[OK == 0,
      n = Input["Input the dimension n - an integer\n"];
      If[n > 0,
	 For[i = 1,i <= n,i++,
	    For[J = 1,J <= n,J++,
	       A[i-1,J-1]=Read[INP,Number];
   	    ];
	 ];
	 OK = 1;
	 Close[INP],
	 Input["Number must be a positive integer\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ],
   Input["This program will end so the input file\n
   can be created.\n
   \n
   Press 1 [enter] to continue\n"];
   OK = 0;
];
If[OK == 1,
   (* Step 1 *)
   A[0,0] = Sqrt[A[0,0]];
   (* Step 2 *)
   For[J = 2,J <= n,J++,
      A[J-1,0]=A[J-1,0]/A[0,0];
   ];
   (* Step 3 *)
   NN = n-1;
   For[i = 2,i<= NN,i++,
      (* Step 4 *)
      KK = i-1;
      S = 0;
      For[K = 1,K <= KK,K++,
	 S = S-A[i-1,K-1]*A[i-1,K-1];
      ];
      A[i-1,i-1] = Sqrt[A[i-1,i-1]+S];
      (* Step 5 *)
      JJ = i+1;
      For[J = JJ,J <= n,J++,
	 S = 0;
	 KK = i-1;
	 For[K = 1,K <= KK,K++,
	    S = S-A[J-1,K-1]*A[i-1,K-1];
	 ];         
         A[J-1,i-1] = (A[J-1,i-1]+S)/A[i-1,i-1];
      ];
   ];
   (* Step 6 *)
   S = 0;
   For[K = 1,K <= NN,K++,
      S = S-A[n-1,K-1]*A[n-1,K-1];
   ];
   A[n-1,n-1] = Sqrt[A[n-1,n-1]+S]; 
   (* Step 7 *)
   FLAG = Input["Select output destination\n
        1. Screen\n
        2. Text file\n
        Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
         For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP,"CHOLESKI FACTORIZATION\n"];
   Write[OUP,"The matrix L output by rows:\n"];
   For[i = 1, i <= n, i++,
      For[J = 1, J <= i, J++,
         Write[OUP,SetPrecision[A[i-1,J-1],9]];
      ];
      Write[OUP,"\n"];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];	
Quit[];