> 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.
> print(`This is the general LU factorization method.\n`):
> print(`Choice of input method`):
> print(`1. input from keyboard - not recommended for large matrices`):
> print(`2. input from a text file`):
> print(`Please enter 1 or 2.`):
> FLAG := scanf(`%d`)[1]: print(`Your input is`): print(FLAG):
> if FLAG = 2 then
>    print(`The array will be input from a text file in the order`):
>    print(`A(1,1), A(1,2), ..., A(1,N), A(2,1), A(2,2), ..., 
>    A(2,N)`):
>    print(`..., A(N,1), A(N,2), ..., A(N,N)\n`):
>    print(`Place as many entries as desired on each line, but separate `):
>    print(`entries with`):
>    print(`at least one blank.`):
>    print(`Has the input file been created? - enter 1 for yes or 2 for no.`):
>    AA := scanf(`%d`)[1]: print(`Your response is`): print(AA):
>    if AA = 1 then
>       print(`Input the file name in the form - drive:\\name.ext`):
>       print(`for example:   A:\\DATA.DTA`):
>       NAME := scanf(`%s`)[1]: print(`The file name is`): print(NAME):
>       INP := fopen(NAME,READ,TEXT):
>       OK := FALSE:
>       while OK = FALSE do
>          print(`Input the dimension n - an integer.`):
>          N := scanf(`%d`)[1]: print(`N is`): print(N):
>          if N > 0 then
>             for I1 from 1 to N do
>             for J from 1 to N do
>                A[I1-1,J-1] := fscanf(INP, `%f`)[1]:
>             od:
>             od:
>             OK := TRUE:
>             fclose(INP):
>          else print(`The number must be a positive integer.\n`):
>          fi:
>       od:
>    else 
>       print(`The program will end so the input file can be created.\n`):
>    fi:
> else
>    OK := FALSE:
>    while OK = FALSE do
>       print(`Input the dimension n - an integer.`):
>       N := scanf(`%d`)[1]: print(`N= `): print(N):
>       if N > 0 then
>          for I1 from 1 to N do
>          for J from 1 to N do
>             print(`input entry in position `,I1,J): 
>             A[I1-1,J-1] := scanf(`%f`)[1]:print(`Data is `):print(A[I1-1,J-1]):
>          od:
>          od:
>          OK := TRUE:
>       else print(`The number must be a positive integer.\n`):
>       fi:
>    od:
> fi:
> if OK = TRUE then
>    OUP := default:
>    fprintf(OUP, `The original matrix - output by rows:\n`):
>    for I1 from 1 to N do
>    for J from 1 to N do
>       fprintf(OUP, ` %11.8f`, A[I1-1,J-1]):
>    od:
>    fprintf(OUP, `\n`):
>    od:
>    print(`Choice of diagonals:\n`):
>    print(`1. Diagonal of L consists of ones\n`): 
>    print(`2. Diagonal of U consists of ones\n`):
>    print(`Please enter 1 or 2.\n`):
>    FLAG := scanf(`%d`)[1]:print(`Your choice is `):print(FLAG):
>    if FLAG = 1 then
>       ISW := 0:
>    else
>       ISW := 1:
>    fi:
>    for I1 from 1 to N do
>       XL[I1-1] := 1:
>    od:
> # Step 1
>    if abs(A[0,0]) <= 1.0e-20 then
>       OK := FALSE:
>    else
> # The entries below the main diagonal will be 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:
>       I1 := 2:
>       while I1 <= M and OK = TRUE do
> # Step 4
>          KK := I1-1:
>          S := 0:
>          for K from 1 to KK do
>             S := S-A[I1-1,K-1]*A[K-1,I1-1]:
>          od:
>          A[I1-1,I1-1] := (A[I1-1,I1-1]+S)/XL[I1-1]:
>          if abs(A[I1-1,I1-1]) <= 1.0e-20 then
>             OK := FALSE:
>          else
> # Step 5
>             JJ := I1+1:
>             for J from JJ to N do
>                SS := 0:
>                S := 0:
>                for K from 1 to KK do
>                   SS := SS-A[I1-1,K-1]*A[K-1,J-1]:
>                   S := S-A[J-1,K-1]*A[K-1,I1-1]:
>                od:
>                if ISW = 0 then
> # Ith row of U  
>                   A[I1-1,J-1] := (A[I1-1,J-1]+SS)/XL[I1-1]:
> # Ith column of L
>                   A[J-1,I1-1] := (A[J-1,I1-1]+S)/A[I1-1,I1-1]:
>                else
> # Ith row of U
>                   A[I1-1,J-1] := (A[I1-1,J-1]+SS)/A[I1-1,I1-1]:
> # Ith column of L
>                   A[J-1,I1-1] := (A[J-1,I1-1]+S)/XL[I1-1]:
>                fi:
>             od:
>          fi:
>          I1 := I1+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
>          print(`Choice of output method:\n`):
>          print(`1. Output to screen\n`):
>          print(`2. Output to text file\n`):
>          print(`Please enter 1 or 2\n`):
>          FLAG := scanf(`%d`)[1]:print(`Your input is `):print(FLAG):
>          if FLAG = 2 then
>             print(`Input the file name in the form - drive:\\name.ext\n`):
>             print(`For example   A:\\OUTPUT.DTA\n`):
>             NAME := scanf(`%s`)[1]:print(`The output file is `):print(NAME):
>             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 I1 from 1 to N do
>          for J from 1 to N do
>             fprintf(OUP, ` %11.8f`, A[I1-1,J-1]):
>          od:
>          fprintf(OUP, `\n`):
>          od:
>          if OUP <> default then
>             fclose(OUP):
>             print(`Output file `,NAME,` created successfully`):
>          fi:
>       fi:
>    fi:
>    if OK = FALSE then
>       print(`Factorization cannot be done\n`):
>    fi:
> fi: