alg065.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> 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.
> alg065 := proc() local AA, NAME, INP, OK, N, I, J, A, V, D, K, FLAG, OUP;
> printf(`This is the LDL^t Method for Positive Definite Matrices.\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;
> else 
> printf(`The program will end so the input file can be created.\n`);
> OK := FALSE;
> fi;
> if OK = TRUE then
> # Step 1
> for I from 1 to N do
> # Step 2
> for J from 1 to I-1 do
> V[J-1] := A[I-1,J-1]*D[J-1];
> od;
> # Step 3
> D[I-1] := A[I-1,I-1];
> for J from 1 to I-1 do
> D[I-1] := D[I-1]-A[I-1,J-1]*V[J-1];
> od;
> # Step 4
> for J from I+1 to N do
> for K from 1 to I-1 do
> A[J-1,I-1] := A[J-1,I-1]-A[J-1,K-1]*V[K-1];
> od;
> A[J-1,I-1] := A[J-1,I-1]/D[I-1];
> od;
> od;
> # Step 5
> 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, `LDL^t FACTORIZATION\n\n`);
> fprintf(OUP, `The matrix L output by rows:\n`);
> for I from 1 to N do
> for J from 1 to I-1 do
> fprintf(OUP, `  %12.8f`, A[I-1,J-1]);
> od;
> fprintf(OUP, `\n`);
> od;
> fprintf(OUP, `The diagonal of D:\n`);
> for I from 1 to N do
> fprintf(OUP, `  %12.8f`, D[I-1]);
> od;
> fprintf(OUP, `\n`);
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg065();