(*  QR ALGORITHM 9.6
*
*  To obtain the eigenvalues of a symmetric, tridiagonal
*    n by n matrix
*
*      a(1)  b(2)
*      b(2)  a(2)  b(3)
*         .     .     .
*           .     .     .
*             .     .     .
*           b(n-1) a(n-1) b(n)
*		      b(n)  a(n)
*
*   INPUT:  n; A(1), ..., A(n) (diagonal of A); B(2), ..., B(n)
*	     (off-diagonal of A); maximum number of iterations M,
*	     tolerance TOL;
*
*  OUTPUT:  Eigenvalues of A or recommended splitting of A,
*	     or a message that the maximum number of iterations
*	     was exceeded.
*)
Print["\n"];
Print["The tridiagonal symmetric array A will be input from\n"];
Print[" a text file in the order:\n"];
Print["\n"];
Print["(diagonal): A(1), A(2), ..., A(n)     \n"];
Print["(subdiagonal): B(2), B(3), ..., B(n)    \n"];
Print["\n"];
Print["Place as many entries as desired on each line, but\n"];
Print["separate entries with at least one blank\n"];
Print["\n"];
Print["\n"];
AA = InputString["This is the QR 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.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++,
	    A[i-1] = Read[INP,Number];
	 ];
	 For[i = 2, i <= n , i++,
	    B[i-1] = Read[INP,Number];
	 ];
      	 OK = 1,      	 
	 Input["The dimension must be greater than one\n
	 \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
   OK = 0;
   While[OK == 0,
      TOL = Input["Input the tolerance\n"];
      If[TOL > 0,
	 OK = 1,
	 Input["Tolerance must be a positive, real number\n
         \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
   OK = 0;
   While[OK == 0,
      L = Input["Input the the maximum number of iterations\n"];
      If[L > 0,
	 OK = 1,
	 Input["The number must be a positive integer\n
         \n
	 Press 1 [enter] to continue\n"];
      ];
   ];
   Close[INP],
   Input["This program will end so the input file\n
   can be created.\n
   \n
   Press 1 [enter] to continue\n"];
];
If[OK == 1,
   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," QR METHOD\n"];
   Write[OUP,"\n"];
   (* Step 1 *)
   SHIFT = 0;
   K = 1;
   (* Step 2 *)
   While[K <= L && OK == 1,
      Write[OUP,"Iteration number ",K," N = ",n];
      Write[OUP,"The array A is now as follows:\n"];
      Write[OUP," Diagonal:\n"];
      For[i = 1, i <= n, i++,
 	 Write[OUP,SetPrecision[A[i-1],9]];
      ];
      Write[OUP,"\n"];
      Write[OUP,"Off diagonal:\n"];
      For[i = 2, i <= n, i++,
	 Write[OUP,SetPrecision[B[i-1],9]];
      ];
      Write[OUP,"\n"]; 
      (* Steps 3-7 - Test for success *)
      (* Step 3 *)
      If[Abs[B[n-1]] <= TOL,
	 A[n-1] = A[n-1]+SHIFT;
	 Write[OUP,"Eigenvalue = \n",SetPrecision[A[n-1],9]];
	 n = n-1;
      ];
      (* Step 4 *)
      If[Abs[B[1]] <= TOL,
	 A[0] = A[0]+SHIFT;
 	 Write[OUP,"Eigenvalue = \n",SetPrecision[A[0],9]];
	 n = n-1;
	 A[0] = A[1];
	 For[J = 2, J <= n, J++,
	    A[J-1] = A[J];
	    B[J-1] = B[J];
	 ];
      ];
      (* Step 5 *)
      If[n == 0,
	 OK = 0;
      ];
      (* Step 6 *)
      If[n == 1,
 	 A[0] = A[0]+SHIFT;
	 Write[OUP,"Eigenvalue = ",SetPrecision[A[0],9]];
	 OK = 0;
      ];
      (* Step 7 *)
      If[OK == 1,
	 M = n-1;
	 If[M >= 2,
	    For[ i = 2, i <= M, i++,
	       If[Abs[B[i-1]] <= TOL,
		  OK = 0;
		  J = i;
	       ];
	    ];
	    If[OK == 0,
	       Write[OUP,"Split the matrix into\n"];
	       For[i = 1, i <= J-1, i++,
	   	  Write[OUP,SetPrecision[A[i-1],9]];
	       ];
	       Write[OUP,"\n"];	
	       For[i = 2, i <= J-1, i++,
		  Write[OUP,SetPrecision[B[i-1],9]];
	       ];
	       Write[OUP,"\n"];
               Write[OUP,"and \n"];
 	       For[i = J, i <= n, i++,
		  Write[OUP,SetPrecision[A[i-1],9]];
	       ];
	       Write[OUP,"\n"];
    	       For[i = J+1, i <= n, i++,
	 	  Write[OUP,SetPrecision[B[i-1],9]];
	       ];
	       Write[OUP,"\n"];
	    ];
   	 ];
      ];
      If[OK == 1,
         (* Step 8 - Compute shift *)
	 B1 = -(A[n-1]+A[n-2]);
	 C1 = A[n-1]*A[n-2]-B[n-1]*B[n-1];
	 D1 = B1*B1-4*C1;
	 If[D1 < 0,
	    Write[OUP,"Problem with complex roots, D1 = \n",D1];
	    OK = 0,
	    D1 = Sqrt[D1];
            (* Step 9 *)
	    If[ B1 > 0,
	       X1 = -2*C1/(B1+D1);
	       X2 = -(B1+D1)/2,
	       X1 = (D1-B1)/2;
	       X2 = 2*C1/(D1-B1);
	    ];
            (* If n=2, then the 2 eigenvalues have been computed *)
            (* Step 10 *)
	    If[n == 2,
	       X1 = N[X1+SHIFT];
	       X2 = N[X2+SHIFT];
	       Write[OUP,"The last two eigenvalues are: ",
                    SetPrecision[X1,9]," ",SetPrecision[X2,9]];
	       OK = 0,
               (* Step 11 *)
	       If[Abs[A[n-1]-X1] > Abs[A[n-1]-X2],
		  X1 = X2;
	       ];
               (* Step 12 - Accumulate shift *)
	       SHIFT = SHIFT+X1;
               (* Step 13 - Perform shift *)
	       For[i = 1, i <= n, i++,
	          d[i-1] = A[i-1]-X1;
	       ];
               (* Steps 14 & 15  - Compute R(K) *)
               (* Step 14 *)
	       X[0] = d[0];
	       Y[0] = B[1];
               (* Step 15 *)
	       For[J = 2, J <= n, J++,
	       	  Z[J-2] = Sqrt[(X[J-2]*X[J-2])+(B[J-1]*B[J-1])];
		  c[J-1] = X[J-2]/Z[J-2];
	  	  S[J-1] = B[J-1]/Z[J-2];
	    	  Q[J-2] = c[J-1]*Y[J-2]+S[J-1]*d[J-1];
		  X[J-1] = c[J-1]*d[J-1]-S[J-1]*Y[J-2];
	 	  If[ J != n,
		     R[J-2] = S[J-1]*B[J];
		     Y[J-1] = c[J-1]*B[J];
		  ];
	       ];
	       M = n-1;
	       MM = n-2;
               (* Steps 16-18 - Compute A(K+1) *)
               (* Step 16 *)
	       Z[n-1] = X[n-1];
	       A[0] = c[1]*Z[0]+S[1]*Q[0];
	       B[1] = S[1]*Z[1];
               (* Step 17 *)
	       If[n > 2,
	          For[J = 2, J <= M, J++,
		     A[J-1] = c[J]*c[J-1]*Z[J-1]+S[J]*Q[J-1];
		     B[J] = S[J]*Z[J];
	 	  ];
	       ];
               (* Step 18 *)
	       A[n-1] = c[n-1]*Z[n-1];
            ];
         ];
      ];
      (* Step 19 *)
      K = K+1;
   ];
   (* Step 20 - Process is complete *)
   If[OK == 1 && K > L,
      Write[OUP,"Maximum number of iterations exceeded.\n"];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
   Print["Output file: ",NAME," created successfully\n"];
   Close[OUP]
   ];
];
Quit[];