(* WIELANDT'S DEFLATION ALGORITHM 9.4
 *
 * To approximate the second most dominant eigenvalue and an
 * associated eigenvector of the n by n matrix A given an
 * approximation LAMBDA to the dominant eigenvalue, an
 * approximation V to a corresponding eigenvector and
 * a vector X belonging to R^(n-1), tolerance TOL, maximum
 * number of iterations N.
 *
 * INPUT:    Dimension n; matrix A; approximate eigenvalue
 *           LAMBDA; approximate eigenvector V belonging
 *           to R^n; vector X belonging to R^(n-1).
 *
 * OUTPUT:   Approximate eigenvalue MU; approximate
 *           eigenvector U or a message that the 
 *           method fails.
*)
(* INPUT *)
AA = InputString["This is Wielandt Deflation.\n\n
Has the Input file been created? - enter y or n.\n"];
If[AA == "Y" || AA == "y",
   NAME = InputString["Input the file name in the form - drive:\\name.ext\n
   for example:   A:\\data\\alg094.dta\n"];
   INP = OpenRead[NAME];
   OK = 0;
   While[OK == 0,
      n = Input["Input the dimension n.\n"];
      If[n > 1,
         OK = 1,
         Input["Dimension must be greater than 1.\n
         Press 1 [enter] to continue\n"];
      ];
   ];
   OK = 0;
   While[OK == 0,
      TOL = Input["Input a positive tolerance for the power method.\n"];
      If[TOL > 0,
         OK = 1,
         Input["Tolerance must be a positive number\n
         Press 1 [enter] to continue\n"];
      ];
   ];
   OK = 0;
   While[OK == 0,
      NN = Input["Input the maximun number of iterations for the power method.\n"];
      If[NN > 0,
         OK = 1,
         Input["The number must be a positive integer.\n
         Press 1 [enter] to continue\n"];
      ];
   ];
   For[i = 1,i <= n,i++,
      For[J = 1 ,J <= n,J++,
         A[i-1,J-1] = Read[INP,Number];
      ];
   ];
   OK = 0;
   For[i = 1,i <=n,i++,
      V[i-1] = Read[INP,Number];
      If[Abs[V[i-1]] > 0,
         OK = 1;
      ];
   ];
   XMU = Read[INP,Number];
   M = n-1;
   If[OK == 1,
      OK = 0;
      For[i = 1,i <= M,i++,
         X[i-1] = Read[INP,Number];
         If[Abs[X[i-1]] > 0,
            OK = 1;
         ];
      ];
   ];
   If[OK == 0,
      Input["All vectors must be nonzero.\n
      \n
      Press 1 [enter] to continue.\n"];
   ];
   Close[INP];
   ,
   OK = 0;
   Input["The program will end so the input file can be created.
   \n
   Press 1 [enter] to continue.\n"];
];
If[OK == 1,
   FLAG = Input["Choice of output method:\n
   1. Output to screen\n
   2. Output to text file\n
   Please enter 1 or 2.\n"];
   If[FLAG == 2,
      NAME = InputString["Input the file name in the form - drive:\\name.ext\n
      for example   A:\\OUTPUT.DTA\n"];
      OUP = OpenWrite[NAME,FormatType->OutputForm],
      OUP = "stdout";
   ];
   Write[OUP,"\nWIELANDT DEFLATION\n\n"];
(*  Step 1  *)
   i = 1;
   AMAX = Abs[V[0]];
   For[J = 2,J <= N,J++,
      If[Abs[V[J-1]] > AMAX,
         i = j;
         AMAX = Abs[V[J-1]];
      ];
   ];
(*  Step 2  *)
   If[i != 1,
      For[K = 1,K <= i-1,K++,
         For[J = 1,J <= i-1,J++,
            B[K-1,J-1] = A[K-1,J-1]-V[K-1]*A[i-1,J-1]/V[i-1];
         ];
      ];
   ];
(* Step 3  *)
   If[i != 1 && i !=n,
      For[K = i,K <= n-1,K++,
         For[J = 1,J <= i-1,J++,
            B[K-1,J-1] = A[K,J-1]-V[K]*A[i-1,J-1]/V[i-1];
            B[J-1,K-1] = A[J-1,K]-V[J-1]*A[i-1,K]/V[i-1];
         ];
      ];
   ];
(* Step 4  *)
   If[i != n,
      For[K = i,K <= n-1,K++,
         For[J = i,J <= n-1,J++,
            B[K-1,J-1] = A[K,J]-V[K]*A[i-1,J]/V[i-1];
         ];
      ];
   ];
(* Step 5   POWER METHOD  *)
   K = 1;
   LP = 1;
   AMAX = Abs[X[0]];
   For[i1 = 2,i1 <= M,i1++,
      If[Abs[X[i1-1]] > AMAX,
         AMAX = Abs[X[i1-1]];
         LP = 1;
      ];
   ];
   DONE = 0;
   For[i1 = 1,i1 <= M,i1++,
       X[i1-1] = X[i1-1]/AMAX
   ];
   While[K <= NN && DONE == 0,
      For[i1 = 1,i1 <=M,i1++,
         Y[i1-1] = 0.0;
         For[J = 1,J <= M,J++,
           Y[i1-1] = Y[i1-1]+B[i1-1,J-1]*X[J-1]
         ];
      ];
      YMU = Y[LP-1];
      LP = 1;
      AMAX = Abs[Y[0]];
      For[i1 = 2,i1 <= M,i1++,
         If[Abs[Y[i1-1]] > AMAX,
            AMAX = Abs[Y[i1-1]];
            LP = i1;
         ];
      ];
      If[AMAX <= 0.000000000001,
         Write[OUP,"Zero eigenvalue - B is singular\n"];
         OK = 0,
         ERR = 0.0;
         For[i1 = 1,i1 <= M,i1++,
            T = Y[i1-1]/Y[LP-1];
            If[Abs[X[i1-1]-T] > ERR,
               ERR = Abs[X[i1-1]-T];99999999999999999999999999
            ];
            X[i1-1] = T;
         ];
         If[ERR < TOL,
            For[i1 = 1,i1 <= M,i1++,
                Y[i1-1] = X[i1-1]
            ];
            DONE = 1,
            K++;
         ];
      ];
   ];
   If[K > NN && OK == 1,
      Write[OUP,"Power method did not converge in ",NN," iterations.\n"];
      OK = 0,
      Write[OUP,"Number of iterations for Power Method = "
            ,K];
   ];
(* Step 6  *)
   If[OK == 1,
      If[i != 1,
         For[K = 1,K <= i-1,K++,
                W[K-1] = Y[K-1]
         ];
      ];
(* Step 7  *)
      W[i-1] = 0.0;
(* Step 8  *)
      If[i != n,
         For[K = i+1,K <= n,K++,
            W[K-1] = Y[K-2]
         ];
      ];
(* Step 9  *)
      S = 0.0;
      For[J = 1,J <= n,J++,
         S = S+A[i-1,J-1]*W[J-1]
      ];
      S = S/V[i-1];
      For[K = 1,K <= n,K++,
(* Compute eigenvector 
   VV is used in place of u here.  *)
         VV[K-1] = (YMU-XMU)*W[K-1]+S*V[K-1];
      ];
      Write[OUP,"The reduced matrix B :\n"];
      For[L1 = i,L1 <= M,L1++,
         For[L2 = 1,L2 <= M,L2++,
            Write[OUP,SetPrecision[B[L1-1,L2-1],9]];
         ];
         Write[OUP,"\n"];
      ];
      Write[OUP,"The Eigenvalue = ",SetPrecision[YMU,10]];
      Write[OUP,"to Tolerance = ",SetPrecision[TOL,9],"\n\n"];
      Write[OUP,"Eigenvector is:"];
      For[i = 1,i <= n,i++,
         Write[OUP,SetPrecision[VV[i-1],10]];
      ];
      Write[OUP,"\n"];
  ];
];
Quit[];