! Power method
! A, w and z - defined in "work"
! Test on "lambda"
! with shift two dominant eigenvalues = complex & of equal moduli

subroutine Power_Complex(n, A, lnit, TOL, lambda, iter, iter_flag, test_val)

complex A(n,n)
complex w0(n), w(n), z0(n), z(n)

complex lambda, lambda0	! Previous value

integer width

complex s


! Set:
s =(0.,0.)

! Define guess vector:
w0 =(1.,1.)

! Normalize guess vector - 2_norm:
z0 =w0/Vector_Norm2(n, w0)


lambda0 =0

! Initialize counter:
90 iter =0

100 iter =iter+1

! Calculate z =Aw
z =matmul(A,z0)


! Calculate "lambda":
imax =maxindex(n,z)
lambda =z(imax)/z0(imax)

! Normalize z:
z =z/Vector_Norm2(n,z)


! Do at least 2 iterations:
if(iter ==1) then
z0 =z
lambda0 =lambda
goto 100
endif


! Test on "lambda":
!test_val =abs(lambda -lambda0)
if(abs(lambda) >1E-6) then 
test_val =abs(lambda -lambda0)/abs(lambda)

else
test_val =abs(lambda -lambda0)
endif

!*
write(11,*) iter, test_val, lambda
!*

IF(test_val <=TOL) THEN
iter_flag =0

!*
!write(11,*) iter, test_val, lambda
!*

if(s /=(0.,0.)) lambda =lambda +s

RETURN	!ok

ELSEIF(iter ==lnit) THEN
	if(s ==(0,0)) then
! Aplly an imag.shift: s =i
	s =(0.,1.)

	do i=1,n
	A(i,i) =A(i,i) -s
	enddo

	goto 90		! iter =0

	else
	iter_flag =-1
	RETURN
	endif


ELSE
z0 =z
lambda0 =lambda
GOTO 100
ENDIF


end subroutine Power_Complex

