! Power method
! A, w and z - defined in "work"
! Test on:
! kodTest =0: lambda
! kodTest #0: z || z0

subroutine Power_Complex(iter, kodTest, ikod, kodNorm, iter_flag, test_val, test_val_min)
use work
complex lambda0	! Previous value

integer width
complex colin(n)

complex eps1

complex lambda1/(3.025604,0.)/
complex rap

colin =(0.,0.)
test_val_min =10.

! Printing:
IF(ikod /=0) THEN
write(1, '(/, a)') ' Guess vector:'
write(1,*) (w0(j), j=1,n)

itl =14+n*15
write(1, '(/, a, t14, a, t<itl>, a)') ' Iteration', 'Current_vector', 'Current_lambda'
ENDIF

! Normalize guess vector - infinity_norm:
if(kodNorm <=1) then
z0 =w0/vNorm$(n, w0)
else
z0 =w0/vNorm2(n, w0)
endif

lambda0 =0

! Initialize counter:
iter =0

100 iter =iter+1

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

! For kodTest =3:
!if(kodTest >=3) then
if(kodTest >=2) then
	if(kodNorm <=1) then
	test_val3 =vNorm$(n, z-lambda*z0)
	else
	test_val3 =vNorm2(n, z-lambda*z0)
	endif
endif

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

! Normalize z:
if(kodNorm <=1) then
z =z/vNorm$(n,z)

else
z =z/vNorm2(n,z)
endif


! Iterations printing:
IF(ikod /=0) THEN
write(1, '(I10, <n>(g15.7,g15.7), g15.7,g15.7)') iter, (z(j), j=1,n), lambda	

!*
if(ikod >=2) then
eps1 =lambda/abs(lambda)
write(1,*) eps1
write(1,*) 'eps1*z0:', eps1*z0
write(1,*)
endif
!*

ENDIF

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


! Test on |lambda-lambda0| or z||z0:
SELECT CASE(kodTest)

CASE(1)	! Test z||z0

do i=1,n
if(abs(z(i)) <TOL .and. abs(z0(i)) <TOL) then
cycle

else
colin(i) =z(i)/z0(i)
endif
enddo

do i=1,n
if(abs(z(i)) <TOL .and. abs(z0(i)) <TOL) then

colin(i) =colin(imax)
endif
enddo

!*
!write(1,*) colin
!*
colin =colin -colin(imax)

	if(kodNorm <=1) then
	test_val =vNorm$(n, colin)

	else
	test_val =vNorm2(n, colin)
	endif

CASE(:0)		! Test "lambda":
test_val =abs(lambda -lambda0)


CASE(2:)	! Test A*z-lambda*z

test_val =test_val3

END SELECT

if(test_val < test_val_min) test_val_min =test_val

if(test_val <=TOL) then
iter_flag =0

!* Test convergence for lambda:
!rap =(lambda -lambda1)/(lambda0 -lambda1)
!write(2,*) iter, rap
!*

RETURN

elseif(iter ==lnit) then
iter_flag =-1
RETURN

else
!* Test convergence for lambda:
!rap =(lambda -lambda1)/(lambda0 -lambda1)
!write(2,*) iter, rap
!*

z0 =z
lambda0 =lambda
GOTO 100
endif


end subroutine Power_Complex

