alg066.txt

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

> restart;
> # CHOLESKI'S ALGORITHM 6.6
> #
> # To factor the positive definite n by n matrix A into LL**T,
> # where L is lower triangular.
> #
> # 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.
> #
> alg066 := proc() local AA, NAME, INP, OK, N, I, J, A, NN, KK, S, K, JJ, FLAG, OUP;
> printf(`This is Choleski 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;
> else 
> printf(`The program will end so the input file can be created.\n`);
> OK := FALSE;
> fi;
> if OK = TRUE then
> # Step 1
> A[0,0] := sqrt(A[0,0]);
> # Step 2
> for J from 2 to N do
> A[J-1,0] := A[J-1,0]/A[0,0];
> od;
> # Step 3
> NN := N-1;
> for I from 2 to NN do
> # Step 4
> KK := I-1;
> S := 0;
> for K from 1 to KK do
> S := S-A[I-1,K-1]*A[I-1,K-1];
> od;
> A[I-1,I-1] := sqrt(A[I-1,I-1]+S);
> # Step 5
> JJ := I+1;
> for J from JJ to N do
> S := 0;
> KK := I-1;
> for K from 1 to KK do
> S := S - A[J-1,K-1]*A[I-1,K-1];
> od;
> A[J-1,I-1] := (A[J-1,I-1]+S)/A[I-1,I-1];
> od;
> od;
> # Step 6
> S := 0;
> for K from 1 to NN do
> S := S-A[N-1,K-1]*A[N-1,K-1];
> od;
> A[N-1,N-1] := sqrt(A[N-1,N-1]+S);
> # 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, `CHOLESKI 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 do
> fprintf(OUP, `  %12.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;
> RETURN(0);
> end;
> alg066();