% FIGURE 3.5.2
      function [EIG,V] = DPOWER(A,N,EIG,V,IUPDAT)
%
%  FUNCTION DPOWER FINDS ONE EIGENVALUE OF A, AND A CORRESPONDING
%    EIGENVECTOR, USING THE SHIFTED INVERSE POWER METHOD.
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    A      - THE N BY N MATRIX.
%
%    N      - THE SIZE OF MATRIX A.
%
%    EIG    - A (COMPLEX) INITIAL GUESS AT      AN EIGENVALUE OF A,
%             AN EIGENVALUE.                    NORMALLY THE ONE CLOSEST
%                                               TO THE INITIAL GUESS.
%
%    V      - A (COMPLEX) STARTING VECTOR       AN EIGENVECTOR OF A,
%             FOR THE SHIFTED INVERSE POWER     CORRESPONDING TO THE
%             METHOD.  IF ALL COMPONENTS OF     COMPUTED EIGENVALUE.
%             V ARE ZERO ON INPUT, A RANDOM
%             STARTING VECTOR WILL BE USED.
%
%    IUPDAT - THE NUMBER OF SHIFTED INVERSE
%             POWER ITERATIONS TO BE DONE
%             BETWEEN UPDATES OF P.  IF
%             IUPDAT=1, P WILL BE UPDATED EVERY
%             ITERATION.  IF IUPDAT > 1000,
%             P WILL NEVER BE UPDATED.
%
%-----------------------------------------------------------------------
%                              EPS = MACHINE FLOATING POINT RELATIVE
%                                    PRECISION
% *****************************
      EPS = eps;
% *****************************
%                             IF V = 0, GENERATE A RANDOM STARTING
%                             VECTOR
      if (V == 0)
      SEED = N+10000;
      DEN = 2.0^31-1.0;
      for I=1:N
         SEED = mod(7^5*SEED,DEN);
         V(I) = SEED/(DEN+1.0);
      end
      end
%                              NORMALIZE V, AND SET VN=V
      VNORM = 0.0;
      for I=1:N
         VNORM = VNORM + abs(V(I))^2;
      end
      VNORM = sqrt(VNORM);
      for I=1:N
         V(I) = V(I)/VNORM;
         VN(I) = V(I);
      end
%                              BEGIN SHIFTED INVERSE POWER ITERATION
      NITER = 1000;
      for ITER=0:NITER
         if (mod(ITER,IUPDAT) == 0)
%                              EVERY IUPDAT ITERATIONS, UPDATE PN
%                              AND SOLVE (A-PN*I)*VNP1 = VN
            PN = EIG;
            for I=1:N
            for J=1:N
               if (I == J)
                  B(I,J) = A(I,J) - PN;
               else
                  B(I,J) = A(I,J);
               end
            end
            end
            [VNP1,IPERM,B] = DLINEQ(B,N,V);
         else
%                              BETWEEN UPDATES, WE CAN USE THE LU
%                              DECOMPOSITION OF B=A-PN*I CALCULATED
%                              EARLIER, TO SOLVE B*VNP1=VN FASTER
            VNP1 = DRESLV(B,N,V,IPERM);
         end
%                              CALCULATE NEW EIGENVALUE ESTIMATE,
%                              PN + (VN*VN)/(VN*VNP1)
         RNUM = 0.0;
         RDEN = 0.0;
         for I=1:N
            RNUM = RNUM + VN(I)*VN(I);
            RDEN = RDEN + VN(I)*VNP1(I);
         end
         R = RNUM/RDEN;
         EIG = PN + R;
%                              SET V = NORMALIZED VNP1
         VNORM = 0.0;
         for I=1:N
            VNORM = VNORM + abs(VNP1(I))^2;
         end
         VNORM = sqrt(VNORM);
         for I=1:N
            V(I) = VNP1(I)/VNORM;
         end
%                              IF R*VNP1 = VN  (R = (VN*VN)/(VN*VNP1) ),
%                              ITERATION HAS CONVERGED.
         ERRMAX = 0.0;
         for I=1:N
            ERRMAX = max(ERRMAX,abs(R*VNP1(I)-VN(I)));
         end
         if (ERRMAX <= sqrt(EPS))
            return
         end
%                              SET VN = V = NORMALIZED VNP1
         for I=1:N
            VN(I) = V(I);
         end
      end
      error('***** INVERSE POWER METHOD DOES NOT CONVERGE *****')