(* LDL^t ALGORITHM 6.5
*
*  To factor the positive definite n by n matrix A into LDL**T, 
*   where L is the lower triangular matrix with ones along the
*   diagonal and D is a diagonal matrix with positive entries 
*   on the diagonal.
*
*   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 and D(i),
*	     1<=i<=n of D
*)
Print["\n"];
Print["A(1,1), A(2,1), ..., A(1,n), A(2,1), A(2,2), ...\n"];
Print["A(2,n), ..., A(n,1), A(n,2), ..., A(n,n)\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 LDL^t Method for Positive Definite Matricies\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 *)
   For[i = 1, i <= n, i++,
      (* Step 2 *)
      For[J = 1, J <= i-1, J++,
	 V[J-1]=A[i-1,J-1]*d[J-1];
      ];
      (* Step 3 *)
      d[i-1]=A[i-1,i-1];
      For[J = 1, J <= i-1, J++,
	 d[i-1]=d[i-1]-A[i-1,J-1]*V[J-1];
      ];
      (* Step 4 *)
      For[J = i+1, J <= n, J++,
	 For[K = 1, K <= i-1, K++,
	    A[J-1,i-1]=A[J-1,i-1]-A[J-1,K-1]*V[K-1];
	 ];
	 A[J-1,i-1]=A[J-1,i-1]/d[i-1];
      ];
   ];
   (* Step 5 *)
   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,"LDL^t FACTORIZATION\n"]; 
   Write[OUP,"\n"];
   For[i = 1, i <= n, i++,
      For[J = 1, J <= i-1, J++,
         Write[OUP,SetPrecision[A[i-1,J-1],9]];
      ];
      Write[OUP,"\n"];
   ];
   Write[OUP,"The diagonal of D:\n"];
   For[i = 1, i <= n, i++,
      Write[OUP,d[i-1]];
   ];
   Write[OUP,"\n"];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP];
   ];
];	
Quit[];