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.\134n`):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)\134n`): 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:\134\134name.ext`): print(`for example: A:\134\134DATA.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.\134n`): fi: od: else print(`The program will end so the input file can be created.\134n`): 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.\134n`): fi: od:fi:if OK = TRUE then OUP := default: fprintf(OUP, `The original matrix - output by rows:\134n`): 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, `\134n`): od: print(`Choice of diagonals:\134n`): print(`1. Diagonal of L consists of ones\134n`): print(`2. Diagonal of U consists of ones\134n`): print(`Please enter 1 or 2.\134n`): 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:\134n`): print(`1. Output to screen\134n`): print(`2. Output to text file\134n`): print(`Please enter 1 or 2\134n`): FLAG := scanf(`%d`)[1]:print(`Your input is `):print(FLAG): if FLAG = 2 then print(`Input the file name in the form - drive:\134\134name.ext\134n`): print(`For example A:\134\134OUTPUT.DTA\134n`): 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\134n\134n`): if ISW = 0 then fprintf(OUP, `The diagonal of L consists of all entries = 1.0\134n`): else fprintf(OUP, `The diagonal of U consists of all entries = 1.0\134n`): fi: fprintf(OUP, `\134nEntries of L below/on diagonal and entries of U above`): fprintf(OUP, `/on diagonal\134n`): fprintf(OUP, `- output by rows in overwrite format:\134n`): 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, `\134n`): 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\134n`): fi:fi:JSFH