> restart:
> # LDL^t ALGORITHM 6.5
> #
> # To factor the positive definite n by n matrix A into LDL**T,
> # where L is a lower triangular matrix with ones along the diagonal
> # and D is a diagonal matrix with positive entries on the
> # diagonal.
> #
> # INPUT:   the dimension n: 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(`This is the LDL^t Method for Positive Definite Matrices.\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:
> # Step 1
> for I1 from 1 to N do
> # Step 2
> for J from 1 to I1-1 do
> V[J-1] := A[I1-1,J-1]*D1[J-1]:
> od:
> # Step 3
> D1[I1-1] := A[I1-1,I1-1]:
> for J from 1 to I1-1 do
> D1[I1-1] := D1[I1-1]-A[I1-1,J-1]*V[J-1]:
> od:
> # Step 4
> for J from I1+1 to N do
> for K from 1 to I1-1 do
> A[J-1,I1-1] := A[J-1,I1-1]-A[J-1,K-1]*V[K-1]:
> od:
> A[J-1,I1-1] := A[J-1,I1-1]/D1[I1-1]:
> od:
> od:
> # Step 5
> 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, `LDL^t FACTORIZATION\n\n`):
> fprintf(OUP, `The matrix L output by rows:\n`):
> for I1 from 1 to N do
> for J from 1 to I1-1 do
> fprintf(OUP, `  %12.8f`, A[I1-1,J-1]):
> od:
> fprintf(OUP, `\n`):
> od:
> fprintf(OUP, `The diagonal of D:\n`):
> for I1 from 1 to N do
> fprintf(OUP, `  %12.8f`, D1[I1-1]):
> od:
> fprintf(OUP, `\n`):
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,` created successfully`):
> fi:
> fi: