(* INVERSE POWER METHOD ALGORITHM 9.3
*
*  To approximate an eigenvalue and an associated 
*    eigenvector of the n by n matrix A given a nonzero vector x: 
*
*   Input:  dimension n; matrix A; vector x; tolerance TOL;
*	     maximum number of iterations NN.
*
*  Output:  approximate eigenvalue MU; approximate eigenvector x;
*            or a message that the maximum number of
*	     iterations was exceeded.
*)
Print["\n"];
Print["A(1,1), A(2,1), ..., A(1,n), A(2,1), A(2,2), ...\n"];
Print["A(2,n), ..., A(n,1), A(n,2), ..., A(n,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["The initial approximation should follow in the same format\n"];
Print["\n"];
Print["\n"];
AA = InputString["This is the Inverse Power 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 - an integer\n"];
      If[n > 0,
	 For[i = 1, i <= n , i++,
	    For[J = 1, J <= n, J++,
	       A[i-1,J-1] = Read[INP,Number];
   	    ];
	 ];
         For[i = 1, i <= n , i++,
	    X[i-1]=Read[INP,Number];
         ];
      	 Close[INP];
         While[OK == 0,
	    TOL=Input["Input the tolerance.\n"];
            If[TOL > 0,
	       OK = 1,
	       Input["Tolerance must be positive\n
	       \n
	       Press 1 [enter] to continue\n"];
            ];
         ];
         OK = 0;
         While[OK == 0,
            NN = Input["Input maximum number of iterations\n
               no decimal points\n"];
            (* Use NN intead of N - maximum number of iterations *)
            If[NN <= 0,
               Input["Must be a positive integer.\n
               Enter 1 [enter] to continue\n"],
               OK = 1;
            ];
         ],
	 Input["The dimension must be a positive integer\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,
   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," INVERSE POWER METHOD\n"];
   Write[OUP,"\n"];
   (* Step 1 - Q could be inputed instead of comuted by deleting 
      the next 10 steps *)
   K = 1;
   Q = 0;
   S = 0;
   For[i = 1, i <= n, i++,
      S=S+X[i-1]*X[i-1];
      For[J = 1, J <= n, J++,
	 Q=Q+A[i-1,J-1]*X[i-1]*X[J-1];
      ];
   ];
   Q=Q/S;
   Print["Q is \n",Q];
   AA = InputString["Input new Q?  Enter 'yes' or 'no'\n"];
   If[ AA=="y" || AA=="Y" || AA=="yes",
      Q=Input["Input new Q\n"];
   ];
   Write[OUP,"Q is \n",Q];
   Write[OUP,"Iteration Eigenvalue Eigenvector\n"];
   (* Step 2 *)
   For[i = 1, i <= n, i++,
      A[i-1,i-1]=A[i-1,i-1]-Q;
   ];
   For[i = 1, i <= n, i++,
      NROW[i-1]=i;
   ];
   OK = 1;
   i = 1;
   M = n-1;
   While[i <= M && OK == 1,
      IMAX = i;
      J = i+1;
      For[IP = J, IP <= n, IP++,
 	 L1=NROW[IMAX-1];
	 L2=NROW[IP-1];
	 If[Abs[A[L2-1,i-1]] > Abs[A[L1-1,i-1]],
	    IMAX=IP;
	 ];
      ];
      If[Abs[A[NROW[IMAX-1]-1,i-1]] <= 10^-20,
         OK = 0;
	 Print["A - Q * i is singular, Q = ",SetPrecision[Q,9]," is an eigenvalue\n"],
         JJ=NROW[i-1];
	 NROW[i-1]=NROW[IMAX-1];
         NROW[IMAX-1]=JJ;
 	 I1=NROW[i-1];
	 For[JJ = J, JJ <= n, JJ++,
	    J1=NROW[JJ-1];
	    A[J1-1,i-1]=A[J1-1,i-1]/A[I1-1,i-1];
	    For[KK = J, KK <= n, KK++,
	       A[J1-1,KK-1]=A[J1-1,KK-1]-A[J1-1,i-1]
                             *A[I1-1,KK-1];
	    ];
	 ];
      ];
      i=i+1;
   ];
   If[Abs[A[NROW[n-1]-1,n-1]] <= 10^-20,
      OK = 0;
      Print["A - Q * i is singular, Q = ",SetPrecision[Q,9]," is an eigenvalue\n"];
   ];
   If[OK == 1,
      (* Step 3 *)
      LP = 1;
      For[i = 2, i <= n, i++,
	 If[Abs[X[i-1]] > Abs[X[LP-1]],
	    LP = i;
	 ];
      ];
      (* Step 4 *)
      AMAX=X[LP-1];
      For[i = 1, i <= n, i++,
   	 X[i-1]=X[i-1]/AMAX;
      ];
      (* Step 5 *)
      While[K <= NN && OK == 1,
         (* Steps 6 & 7 *)
	 For[i = 1, i <= n, i++,
	    B[i-1]=X[i-1];
	 ];
	 M=n-1;
	 For[i = 1, i <= M, i++,
	    J=i+1;
	    I1=NROW[i-1];
	    For[JJ = J, JJ <= n, JJ++,
	       J1=NROW[JJ-1];
	       B[J1-1]=B[J1-1]-A[J1-1,i-1]*B[I1-1];
	    ];
	 ];
	 N1=NROW[n-1];
	 Y[n-1]=B[N1-1]/A[N1-1,n-1];
 	 L=n-1;
	 For[K1 = 1, K1 <= L, K1++,
	    J=L-K1+1;
	    JJ=J+1;
	    N2=NROW[J-1];
	    Y[J-1]=B[N2-1];
	    For[KK = JJ, KK <= n, KK++,
	       Y[J-1]=Y[J-1]-A[N2-1,KK-1]*Y[KK-1];
	    ];
	    Y[J-1]=Y[J-1]/A[N2-1,J-1];
         ];
         (* Step 8 *)
	 YMU=Y[LP-1];
         (* Steps 9 & 10 *)
	 LP=1;
	 For[i = 2, i <= n, i++,
	    If[Abs[Y[i-1]] > Abs[Y[LP-1]],
	       LP=i;
	    ];
	 ];
	 AMAX=Y[LP-1];
	 ERR=0;
         For[i = 1, i <= n, i++,  
	    T=Y[i-1]/AMAX;
	    If[Abs[X[i-1]-T] > ERR,
	       ERR = Abs[X[i-1]-T];
	    ];
	    X[i-1]=T;
	 ];
	 YMU=1/YMU+Q;
         (* step 11 *)
	 Write[OUP,K,"        ",SetPrecision[YMU,10]];
	 For[i = 1, i <= n, i++,
	    Write[OUP,"                      \n"
                                ,SetPrecision[X[i-1],10]];
	 ];
	 Write[OUP,"\n"];
	 If[ERR < TOL,
	    OK = 0;
	    Write[OUP,"Eigenvalue = \n",SetPrecision[YMU,10]];
	    Write[OUP,"to tolerance = \n",SetPrecision[TOL,9]];
	    Write[OUP,"obtained on iteration number ",K];
	    Write[OUP,"\n"];
	    Write[OUP,"Unit eigenvector is : \n"];
	    For[i = 1, i <= n, i++,
	       Write[OUP,SetPrecision[X[i-1],10]];
	    ];
	    Write[OUP,"\n"],
            (* Step 12 *)
	    K=K+1;
	 ];
      ];
      If[K > NN,
	 Write[OUP,"No convergence in ",NN," iterations\n"];
      ];
   ];
   If[OUP == "OutputStream[",NAME," 3]",
      Print["Output file: ",NAME," created successfully\n"];
      Close[OUP]
   ];
];
Quit[];