alg095.txt

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

> restart;
> # HOUSEHOLDER'S ALGORITHM 9.5
> #
> # To obtain a symmetric tridiagonal matrix A(n-1) similar
> # to the symmetric matrix A = A(1), construct the following
> # matrices A(2),A(3),...,A(n-1) where A(K) = A(I,J)**K, for
> # each K = 1,2,...,n-1:
> #
> # INPUT:   Dimension n; matrix A.
> #
> # OUTPUT:  A(n-1) (At each step, A can be overwritten.)
> alg095 := proc() local OK, AA, NAME, INP, N, I, J, A, K, Q, KK, S, RSQ, V, U, PROD, Z, L, FLAG, OUP;
> printf(`This is the Householder Method.\n`);
> OK := FALSE;
> printf(`The symmetric array A will be input from a text file\n`);
> printf(`in the order:\n`);
> printf(`              A(1,1), A(1,2), A(1,3), ..., A(1,n),\n`);
> printf(`                      A(2,2), A(2,3), ..., A(2,n),\n`);
> printf(`                              A(3,3), ..., A(3,n),\n`);
> printf(`                                      ..., 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.\n`);
> N := scanf(`%d`)[1];
> if N > 1 then
> for I from 1 to N do
> for J from I to N do
> A[I-1,J-1] := fscanf(INP, `%f`)[1];
> A[J-1,I-1] := A[I-1,J-1];
> od;
> od;
> fclose(INP);
> OK := TRUE;
> else
> printf(`Dimension must be greater than 1.\n`);
> fi;
> od;
> else
> printf(`The program will end so the input file can be created.\n`);
> fi;
> if OK = TRUE then
> # Step 1
> for K from 1 to N-2 do
> Q := 0;
> KK := K+1;
> # Step 2
> for I from KK to N do
> Q := Q+A[I-1,K-1]*A[I-1,K-1];
> od;
> # Step 3
> if abs(A[K,K-1]) <= 1.0e-20 then
> S := sqrt(Q);
> else
> S := A[K,K-1]/abs(A[K,K-1])*sqrt(Q);
> fi;
> # Step 4
> RSQ := (S+A[K,K-1])*S;
> # Step 5
> V[K-1] := 0;
> V[K] := A[K,K-1]+S;
> for J from K+2 to N do
> V[J-1] := A[J-1,K-1];
> od;
> # Step 6
> for J from K to N do
> U[J-1] := 0;
> for I from KK to N do
> U[J-1] := U[J-1]+A[J-1,I-1]*V[I-1];
> od;
> U[J-1] := U[J-1]/RSQ;
> od;
> # Step 7
> PROD := 0;
> for I from K+1 to N do
> PROD := PROD + V[I-1]*U[I-1];
> od;
> # Step 8
> for J from K to N do
> Z[J-1] := U[J-1] - 0.5*PROD*V[J-1]/RSQ;
> od;
> # Step 9
> for L from K+1 to N-1 do
> # Step 10
> for J from L+1 to N do
> A[J-1,L-1] := A[J-1,L-1]-V[L-1]*Z[J-1]-V[J-1]*Z[L-1];
> A[L-1,J-1] := A[J-1,L-1];
> od;
> # Step 11
> A[L-1,L-1] := A[L-1,L-1] - 2*V[L-1]*Z[L-1];
> od;
> # Step 12
> A[N-1,N-1] := A[N-1,N-1]-2*V[N-1]*Z[N-1];
> # Step 13
> for J from K+2 to N do
> A[K-1,J-1] := 0;
> A[J-1,K-1] := 0;
> od;
> # Step 14
> A[K,K-1] := A[K,K-1]-V[K]*Z[K-1];
> A[K-1,K] := A[K,K-1];
> od;
> # Step 15
> 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, `HOUSEHOLDER METHOD\n\n`);
> fprintf(OUP, `The similar tridiagonal matrix follows - output by rows\n\n`);
> for I from 1 to N do
> for J from 1 to N do
> fprintf(OUP, ` %11.8f`, A[I-1,J-1]);
> od;
> fprintf(OUP, `\n\n`);
> od;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg095();