(* HOUSEHOLDER'S ALGORITHM 9.5
*
*  To obtain a symmetrical tridiagonal matrix A(n-1) similar
*   to the symmertic matrix A=A(1), construct the following
*   matricies 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.)
*)
Print["\n"];
Print["The symmetric array A will be input from a text\n"];
Print[" file in the order:\n"];
Print["\n"];
Print["      A(1,1), A(1,2), A(1,3), ..., A(1,n),\n"];
Print["              A(2,2), A(2,3), ..., A(2,n),\n"];
Print["                      A(3,3), ..., A(3,n),\n"];
Print["                              ..., A(n,n)\n"];
Print["\n"];
Print["\n"];
AA = InputString["This is the Householder Method\n
   The array will be input from a text file in the order:(see screen)\n
   Has the input file been created?\n
   Enter 'yes' or 'no'\n"];
OK = 0;
If[AA == "yes" || AA == "y" || AA == "Y",
   NAME=InputString["Input the file name in the form - \n
      drive:\ name.ext      for example\n
         A:\\DATA\\ALG095.DTA\n"];
   INP = OpenRead[NAME];
   OK = 0;
   While[OK == 0,
      n = Input["Input the dimension n \n"];
      If[n > 1,
	 For[i = 1, i <= n , i++,
	    For[J = i, J <= n, J++,
	       A[i-1,J-1] = Read[INP,Number];
	       A[J-1,i-1] = A[i-1,J-1];
   	    ];
	 ];
	 Close[INP];
         OK = 1,      	 
	 Input["The dimension must be greater than one\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ],  
   Input["This program will end so the input file\n
   can be created.\n
   \n
   Press 1 [enter] to continue\n"];
];
If[OK == 1,
   (* Step 1 *)
   For[K = 1, K <= n-2, K++,
      Q = 0;
      KK = K+1;
      (* Step 2 *)
      For[i = KK, i <= n, i++,
	 Q = Q+A[i-1,K-1]*A[i-1,K-1];
      ];
      (* Step 3 - S is used in place of alpha *)
      If[Abs[A[K,K-1]] <= 0,
	 S = Sqrt[Q],
	 S = A[K,K-1]/Abs[A[K,K-1]]*Sqrt[Q];
      ];
      (* Step 4 *)
      RSQ = (S+A[K,K-1])*S;
      (* Step 5 *)
      V[K-1] = 0;
      V[K] = A[K,K-1]+S;
      For[J = K+2, J <= n, J++,
	 V[J-1] = A[J-1,K-1];
      ];
      (* Step 6 *)
      For[J = K, J <= n, J++,
	 U[J-1] = 0;
	 For[i = KK, i <= n, i++,
	    U[J-1] = U[J-1]+A[J-1,i-1]*V[i-1];
	 ];
	 U[J-1] = U[J-1]/RSQ;
      ];
      (* Step 7 *)
      PROD = 0;
      For[i = K+1, i <= n, i++,
	 PROD = PROD+V[i-1]*U[i-1];
      ];
      (* Step 8 *)
      For[J = K, J <= n, J++,
	 Z[J-1] = U[J-1]-0.5*PROD*V[J-1]/RSQ;
      ];
      (* Step 9 *)
      For[L = K+1, L <= n-1, L++,
         (* Step 10 *)
	 For[J = L+1, J <= n, J++,
	    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];
	 ];
         (* Step 11 *)
	 A[L-1,L-1] = A[L-1,L-1]-2*V[L-1]*Z[L-1];
      ];
      (* Step 12 *)
      A[n-1,n-1] = A[n-1,n-1]-2*V[n-1]*Z[n-1];
      (* Step 13 *)
      For[J = K+2, J <= n, J++,
 	 A[K-1,J-1] = 0;
	 A[J-1,K-1] = 0;
      ];
      (* Step 14 *)
      A[K,K-1] = A[K,K-1]-V[K]*Z[K-1];
      A[K-1,K] = A[K,K-1];
   ];   
   (* Step 15 *)
   FLAG = Input["Select output destination\n
      1. Screen\n
      2. Text file\n
    Enter 1 or 2\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name\n
         For example:   output.dta\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP," HOUSEHOLDER METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"The similar tridiagonal matrix follows \n
                      - out by rows\n"];
   Write[OUP,"\n"];
   For[i = 1, i <= n, i++,
      For[J = 1, J <= n, J++,
          Write[OUP,SetPrecision[A[i-1,J-1],9]];	       
      ];
      Write[OUP,"\n"];
      Write[OUP,"\n"];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];
Quit[];