(*DIRECT FACTORIZATION ALGORITHM 6.4
*
*   To factor the n by n matrix A= (A(i,J)) into the product of 
*   the lower triangular matrix L = (L(i,J)) and U = (U(i,J)),
*   that is A=LU, where the main diagonal of either L or U 
*   consists of all ones:
*
*   Input:  dimension n; the entries A(i,J), 1<=i, J<=n, of A;
*	    The diagonal L(1,1), ..., L(n,n) of L or the diagonal
*	    U(1,1), ..., U(n,n) of U
*
*   Output:  the entries L(i,J), 1<=J<=i, 1<=i<=n of L and the
*	     entries U(i,J), i<=J<=n, 1<=i<=n of U
*)
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 General LU 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"];
      ];
      FLAG = Input["Choice of diagonals:\n
             1. Diagonal of L consists of ones\n
	     2. Diagonal of U consists of ones\n
             Please enter 1 or 2.\n"];
      If[FLAG == 1,
	 ISW = 0,
	 ISW = 1;
      ]; 	
   ],
   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,
   For[i = 1,i <= n,i++,
      XL[i-1] = 1;
   ];
   (* Step 1 *)
   If[Abs[A[0,0]] <=10^-20,
      OK = 0,
      (* The entries below the main diagonal will be placed 
         in the corresponding entries of A; the entries of U
         above the main diagonal will be placed in the corres-
         ponding entries of A; the main diagonal which was not
         inputed will become the main diagonal of A; the input
         main diagonal of L or U is placed in XL *)      
      A[0,0] = A[0,0]/XL[0];
      (* Step 2 *)
      For[J = 2,J <= n,J++,
	 If[ISW == 0,
            (* First row of U *)
	    A[0,J-1] = A[0,J-1]/XL[0];
            (* First column of L *)
	    A[J-1,0] = A[J-1,0]/A[0,0],
            (* first row of U *)
	    A[0,J-1] = A[0,J-1]/A[0,0];
            (* First column of L *)
	    A[J-1,0] = A[J-1,0]/XL[0];
	 ];
      ];
      (* Step 3 *)
      M = n-1;
      i = 2;
      While[i <= M && OK == 1,
         (* Step 4 *)
	 KK = i-1;
	 S = 0;
	 For[K = 1,K <= KK,K++,
	    S = S-A[i-1,K-1]*A[K-1,i-1];
	 ];
	 A[i-1,i-1] = (A[i-1,i-1]+S)/XL[i-1];
	 If[Abs[A[i-1,i-1]] <= 10^-20,
	    OK = 0,
            (* Step 5 *)
	    JJ = i+1;
	    For[J = JJ,J <= n,J++,
	       SS = 0;
	       S = 0;
	       For[K = 1,K <= KK,K++,
		  SS = SS-A[i-1,K-1]*A[K-1,J-1];
		  S = S-A[J-1,K-1]*A[K-1,i-1];
	       ];
	       If[ISW == 0,
                  (* Ith row of U *)
		  A[i-1,J-1] = (A[i-1,J-1]+SS)/XL[i-1];
                  (* Ith column of L *)
	  	  A[J-1,i-1] = (A[J-1,i-1]+S)/A[i-1,i-1],
                  (* Ith row of U *)
		  A[i-1,J-1] = (A[i-1,J-1]+SS)/A[i-1,i-1];
                  (* Ith column of L *)
	  	  A[J-1,i-1] = (A[J-1,i-1]+S)/XL[i-1];
	       ];
	    ];
         ];
	 i = i+1;
      ];
      If[OK == 1,
         (* Step 6 *)
	 S = 0;
	 For[K = 1,K <= M,K++,
	    S = S-A[n-1,K-1]*A[K-1,n-1]
	 ];
	 A[n-1,n-1] = (A[n-1,n-1]+S)/XL[n-1];
         (* If A(N-1,N-1)=0 then A=LU but the matrix is 
            singular. Process is complete, all entries 
            of A have been determined. *)
         (* 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,"GENERAL LU FACTORIZATION\n"]; 
         Write[OUP,"\n"];
	 If[ISW == 0,
            Write[OUP,"The diagonal of L consists of all entries = 1.0\n"],
	    Write[OUP,"The diagonal of U consists of all entries = 1.0\n"];
	 ];
         Write[OUP,"\n"];
         Write[OUP,"Entries of L below/on diagonal and entries\n
            of U above/on diagonal \n
            output by rows in overwrite format:\n"];
         For[i = 1,i <= n,i++,
	    For[J = 1,J <= n,J++,
               Write[OUP,SetPrecision[A[i-1,J-1],9]];
            ];
            Write[OUP,"\n"];
         ];
	 Write[OUP,"\n"];
         If[OUP == "OutputStream[",NAME," 3]",
            Print["Output file: ",NAME," created successfully\n"];
            Close[OUP];
         ];
      ];
   ];
   If[OK == 0,
      Print["System has no unique solution\n"];
   ];
];
Quit[];