(* SYMMETRIC POWER METHOD ALGORITHM 9.2
*
*  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 Symmetric 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];
   	    ];
	 ];
         (* Because of the way in which the subroutine norm computes
            the I2 norm of a vector, the initial input is into 
            Y and X is initialized at the zero vector *)
         For[i = 1, i <= n , i++,
	    Y[i-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 for N the 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," SYMMETRIC POWER METHOD\n"];
   Write[OUP,"\n"];
   Write[OUP,"iter    approx         approx eigenvector\n"];
   Write[OUP,"        eigenvalue\n"];
   (* Step 1 *)
   K=1;
   XL=0;
   (* Norm computes the norm of the vector X and returns X 
      divided by its norm in the variable Y *)
   For[i = 1, i <= n, i++,
      XL = XL+Y[i-1]*Y[i-1];
   ];
   (* 2-norm of Y *)
   XL=Sqrt[XL];
   ERR = 0;
   If[XL > 0,
      For[i = 1, i <= n, i++,
         T=Y[i-1]/XL;
	 ERR=ERR+(X[i-1]-T)*(X[i-1]-T);
	 X[i-1]=T;
      ];
      (* X has a 2-norm of 1.0 *)
      ERR=Sqrt[ERR],
      Print["A has a zero eigenvalue - select new vector\n
	and begin again\n"];
      OK = 0;
   ];
   If[OK == 1,
      (* Step 2 *)
      While[K <= NN && OK == 1,
         (* Steps 3 & 4 *)
	 YMU=0;
         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];
	    ];
	    YMU=YMU+X[i-1]*Y[i-1];
	 ];
         (* Step 5 is accomplished in the subroutine norm *)
         (* Step 6 *)
	 XL=0;
	 For[i = 1, i <= n, i++,
	    XL = XL+Y[i-1]*Y[i-1];
	 ];
         (* 2-norm of Y *)
	 XL=Sqrt[XL];
	 ERR=0;
	 If[XL > 0,
	    For[i = 1, i <= n, i++,
	       T=Y[i-1]/XL;
	       ERR=ERR+(X[i-1]-T)*(X[i-1]-T);
	       X[i-1]=T;
	    ];
            (* X has a 2-norm of 1.0 *)
	    ERR=Sqrt[ERR],
	    Print["A has a zero eigenvalue - select new vector\n
	     and begin again\n"];
	    OK = 0;
	 ];
	 Write[OUP,K,"     ",SetPrecision[YMU,9]];
         For[i = 1, i <= n, i++,
            Write[OUP,SetPrecision[X[i-1],9]];
         ];
	 Write[OUP,"\n"];
	 If[OK == 1,
            (* Step 7 - Procedure completed successfully *)
	    If[ERR < TOL,      
	       Write[OUP,"\n"];
               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,
   	       K=K+1;
            ];
	 ];
      ];
      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[];