(* POWER METHOD ALGORITHM 9.1
*
*  To approximate the dominant 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 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,
            (* Use NN for n *)
            NN = Input["Input maximum number of iterations\n
               no decimal points\n"];
            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,"POWER METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"iter    approx         approx eigenvector\n"];
   Write[OUP,"        eigenvalue\n"];
   (* Step 1 *)
   K=1;
   (* Step 2 *)
   LP=1;
   AMAX=Abs[X[0]];
   For[i = 2, i <= n, i++,
      If[Abs[X[i-1]] > AMAX,
	 AMAX=Abs[X[i-1]];
	 LP=i;
      ];
   ];
   (* Step 3 *)
   For[i = 1, i <= n, i++,
      X[i-1]=X[i-1]/AMAX;
   ];
   (* Step 4 *)
   While[K <= NN && OK == 1,
      (* Step 5 *)
      For[i = 1, i <= n, i++,
	 Y[i-1]=0;
         For[ J = 1, J <= n, J++, 
	    Y[i-1]=Y[i-1]+A[i-1,J-1]*X[J-1];
    	 ];
      ];
      (* Step 6 *)
      YMU=Y[LP-1];
      (* Step 7 *)
      LP=1;
      AMAX=Abs[Y[0]];
      For[i = 2, i <= n, i++,
  	 If[Abs[Y[i-1]] > AMAX,
	    AMAX = Abs[Y[i-1]];
	    LP=i;
	 ];
      ];
   (* Step 8 *)
   If[AMAX <= 0,
      Write[OUP,"0 eigenvalue - select another
 	 initial vector and begin again\n"];
      OK = 0,
   ];
   (* Step 9 *)
   ERR = 0;
   For[i = 1, i <= n, i++,
      T=Y[i-1]/Y[LP-1];
      If[Abs[X[i-1]-T] > ERR,
	 ERR = Abs[X[i-1]-T];
      ];
      X[i-1]=T;
   ];
   Write[OUP,K,"     ",SetPrecision[YMU,9]];
   For[i = 1, i <= n, i++,
      Write[OUP,SetPrecision[X[i-1],9]];
   ];
   Write[OUP,"\n"];
   (* Step 10 *)
   If[ERR <= TOL,
      Write[OUP,"\n"];
      Write[OUP,"The eigenvalue = \n",SetPrecision[YMU,9]];
      Write[OUP,"to tolerance \n",SetPrecision[TOL,9]];
      Write[OUP,"obtained on iteration number = \n",K];
      Write[OUP,"\n"];
      Write[OUP,"Unit eigenvector is :\n"];
      Write[OUP,"\n"];
      For[i = 1, i <= n, i++,
	 Write[OUP,SetPrecision[X[i-1],9]];
      ];
      Write[OUP,"\n"];
      OK = 0;
   ];
   (* Step 11 *)
   K=K+1;
];
(* Step 12 *)
If[K > NN,
   Write[OUP,"Method did not converge within ",NN," iterations.\n"];
];
If[OUP == "OutputStream[",NAME," 3]",
   Print["Output file: ",NAME," created successfully\n"];
   Close[OUP]
];
];
Quit[];