function [ x, r ] = inpows ( A, xold, rold )
%
% function [ x, r ] = inpows ( A, xold, rold )
%
%  Discussion:
%
%    INPOWS carries out one step of the shifted inverse power method.
%
%  Modified:
%
%    09 April 2000
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A(N,N), the matrix.
%
%    Input, real XOLD(N), an estimate for the eigenvector.
%
%    Input, real ROLD, an estimate for the eigenvalue.
%
%    Output, real X(N), the value of inv(A-rold*I)*xold, normalized to have unit L2 norm.
%
%    Output, real R, the Rayleigh estimate of the eigenvalue.
%
[n,n] = size ( A );

A = A - rold * eye(n);

[P, L, U ] = ge_pp ( A )
%
%  Multiply:                   x = inv ( A - rold * I ) * xold  
%  <=>  
%  Solve:   ( A - rold * I ) * x =                        xold
%
x = plu_solve ( P, L, U, xold );
%
%  The matrix here is really (A - rold * I).
%  We have to add ROLD back to the Rayleigh quotient estimate.
%
r = rold + ( x' * A * x ) / ( x' * x );

x = x / norm ( x );
%
%  X and -X are both eigenvectors.  To make it easier to detect
%  convergence, pick the one whose maximum is positive.
%
x_max = max ( x );
if ( x_max < 0 )
  x = - x;
end

