> restart;
> # 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 the upper triangular
> # matrix 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.
> alg064 := proc() local AA, NAME, INP, OK, N, I, J, A, FLAG, ISW, XL, M, KK, S, K, JJ, SS, OUP;
> printf(`This is the general LU factorization method.\n`);
> printf(`The array will be input from a text file in the order:\n`);
> printf(`A(1,1), A(1,2), ..., A(1,N), A(2,1), A(2,2), ..., 
> A(2,N),\n`);
> printf(`..., A(N,1), A(N,2), ..., A(N,N)\n\n`);
> printf(`Place as many entries as desired on each line, but separate `);
> printf(`entries with\n`);
> printf(`at least one blank.\n\n\n`);
> printf(`Has the input file been created? - enter Y or N.\n`);
> AA := scanf(`%c`)[1];
> if AA = "Y" or AA = "y" then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:   A:\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> INP := fopen(NAME,READ,TEXT);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the dimension n - an integer.\n`);
> N := scanf(`%d`)[1];
> if N > 0 then
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := fscanf(INP, `%f`)[1];
> od;
> od;
> OK := TRUE;
> fclose(INP);
> else printf(`The number must be a positive integer.\n`);
> fi;
> od;
> printf(`Choice of diagonals:\n`);
> printf(`1. Diagonal of L consists of ones\n`); 
> printf(`2. Diagonal of U consists of ones\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 1 then
> ISW := 0;
> else
> ISW := 1;
> fi
> else 
> printf(`The program will end so the input file can be created.\n`);
> OK := FALSE;
> fi;
> if OK = TRUE then
> for I from 1 to N do
> XL[I-1] := 1;
> od;
> # Step 1
> if abs(A[0,0]) <= 1.0e-20 then
> OK := FALSE;
> else
> # The entries below the main diagonal will placed in the corresponding
> # entries in the matrix A
> A[0,0] := A[0,0]/XL[0];
> # Step 2
> for J from 2 to N do
> if ISW = 0 then
> # 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];
> else
> # 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];
> fi;
> od;
> # Step 3
> M := N-1;
> I := 2;
> while I <= M and OK = TRUE do
> # Step 4
> KK := I-1;
> S := 0;
> for K from 1 to KK do
> S := S-A[I-1,K-1]*A[K-1,I-1];
> od;
> A[I-1,I-1] := (A[I-1,I-1]+S)/XL[I-1];
> if abs(A[I-1,I-1]) <= 1.0e-20 then
> OK := FALSE;
> else
> # Step 5
> JJ := I+1;
> for J from JJ to N do
> SS := 0;
> S := 0;
> for K from 1 to KK do
> SS := SS-A[I-1,K-1]*A[K-1,J-1];
> S := S-A[J-1,K-1]*A[K-1,I-1];
> od;
> if ISW = 0 then
> # 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];
> else
> # 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];
> fi;
> od;
> fi;
> I := I+1;
> od;
> if OK = TRUE then
> # Step 6
> S := 0;
> for K from 1 to M do
> S := S-A[N-1,K-1]*A[K-1,N-1];
> od;
> 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
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`For example   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `GENERAL LU FACTORIZATION\n\n`);
> if ISW = 0 then 
> fprintf(OUP, `The diagonal of L consists of all entries = 1.0\n`);
> else
> fprintf(OUP, `The diagonal of U consists of all entries = 1.0\n`);
> fi;
> fprintf(OUP, `\nEntries of L below/on diagonal and entries of U above`);
> fprintf(OUP, `/on diagonal\n`);
> fprintf(OUP, `- output by rows in overwrite format:\n`);
> for I from 1 to N do
> for J from 1 to N do
> fprintf(OUP, ` %11.8f`, A[I-1,J-1]);
> od;
> fprintf(OUP, `\n`);
> od;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> fi;
> if OK = FALSE then
> printf(`System has no unique solution\n`);
> fi;
> fi;
> RETURN(0);
> end;
> alg064();