> 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.)
> printf(`This is the Householder Method.`):
> OK := FALSE:
> 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 symmetric array A will be input from a text file`):
> print(`in the order:`):
> print(`              A(1,1), A(1,2), A(1,3), ..., A(1,n),`):
> print(`                      A(2,2), A(2,3), ..., A(2,n),`):
> print(`                              A(3,3), ..., A(3,n),`):
> print(`                                      ..., A(n,n)`):
> print(`Place as many entries as desired on each line, but separate `):
> print(`entries with 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.`):
> N := scanf(`%d`)[1]:print(`n is`): print(N):
> if N > 1 then
> for I1 from 1 to N do
> for J from I1 to N do
> A[I1-1,J-1] := fscanf(INP, `%f`)[1]:
> A[J-1,I1-1] := A[I1-1,J-1]:
> od:
> od:
> fclose(INP):
> OK := TRUE:
> else
> print(`Dimension must be greater than 1.`):
> fi:
> od:
> else
> print(`The program will end so the input file can be created.`):
> 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:
> fi:
> if OK = TRUE then
> # Step 1
> for K from 1 to N-2 do
> Q := 0:
> KK := K+1:
> # Step 2
> for I1 from KK to N do
> Q := Q+A[I1-1,K-1]*A[I1-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 I1 from KK to N do
> U[J-1] := U[J-1]+A[J-1,I1-1]*V[I1-1]:
> od:
> U[J-1] := U[J-1]/RSQ:
> od:
> # Step 7
> PROD := 0:
> for I1 from K+1 to N do
> PROD := PROD + V[I1-1]*U[I1-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
> 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]:
> 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]:
> OUP := fopen(NAME,WRITE,TEXT):
> else
> OUP := default:
> fi:
> fprintf(OUP, `HOUSEHOLDER METHOD\n`):
> fprintf(OUP, `The similar tridiagonal matrix follows - 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:
> if OUP <> default then
> fclose(OUP):
> print(`Output file `,NAME,`  created successfully`):
> fi:
> fi: