! MultiVector_Power method
! A, W and Z - defined in "work"
! ne <=n eigenvectors are required
! Test on |w-w0|
!
subroutine Sim_Iter(iter, kodTest, kodPrint, iter_flag, test_val)
use work

implicit real(8) (a-h, o-z)

!real lambda0(ne)	! Previous value - in work

!real colin(ne), test_val
logical herm

complex(8) colin(n)
real(8) test(n), test_val
complex(8) z(n), z0(n)
integer width


n3 =min(ne,3)

! Printing:
IF(kodPrint >=2) THEN
write(1, '(/,a)' ) ' Guess Vector''s Matrix:'

!do i=1,n
!write(1, '(t2, <n3>(g14.7, g14.7, 1x))') (W0(i,j), j=1,ne)
!enddo

DO J =1,ne,n3
NJ =J+n3-1
if(NJ >ne) NJ =ne

do i =1,n
write(1, '(<n3>(g14.7, g14.7, 1x) )') (W0(i,jw), jw =J,NJ)
enddo
write(1,*)
ENDDO


itl =n3*28 +11
write(1, '(/, a, t8, a, t<itl>, a)') ' Iter', 'Current Vector''s Matrix', 'Current Eigenvalues'
ENDIF

! Perform Gramm-Schmidt Normalization of W0:
call Gram_Schmidt(n, ne, W0, Tol)

lambda =(0.,0.)
W =W0
colin =0.

! Initialize counter:
iter =0
iter_flag =-1

100 iter =iter+1

if(n >=20 .and. lnit >=100 .and. mod(iter, 100)==0) then
iw =width(iter)
write(*, '(/, a, I<iw>)') ' Iteration: ', iter
endif


lambda0 =lambda
W0 =W

! Calculate z =Aw
W =matmul(A, W0)

! Calculate lambdas by Rayleigh quotient:
do j=1,ne
lambda(j) =dot_product(W0(:,j), W(:,j))
enddo


! Gramm-Schmidt
call Gram_Schmidt(n, ne, W, Tol)

! Iterations printing:
IF(kodPrint >=2) THEN
iwi =width(iter)+1
write(1, '(I<iwi>)') iter

!do i=1,n
!if(i<=ne) then
!write(1, '(t5, <ne>(g16.7, g16.7), g16.7, g16.7)') (W(i,j), j=1,ne), lambda(i)

!else
!write(1, '(t5, <n3>(g16.7, g16.7))' ) (W(i,j), j=1,ne)
!endif
!enddo
!ENDIF


DO J =1,ne,n3
NJ =J+n3-1
if(NJ >ne) NJ =ne

do i =1,n
if(j<=n3) then
write(1, '(<n3>(g15.7, g15.7), g16.7,g15.7 )') (W(i,jw), jw =J,NJ), lambda(i)

else
write(1, '(<n3>(g15.7, g15.7), g16.7,g15.7 )') (W(i,jw), jw =J,NJ)
endif

enddo
write(1,*)
ENDDO

ENDIF

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

! Test |lambda -lambda0| <= TOL, or z(k+1) || z(k):
SELECT CASE(kodTest)

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

CASE(1:)	! Test z||z0

DO je =1,ne
z =W(:,je); z0 =W0(:,je);

imax =maxindex_z(z)

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

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

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

colin(i) =colin(imax)
endif
enddo

colin =colin -colin(imax)
test(je) =vNorm2(n, colin)

ENDDO

test_val =maxval(test)

END SELECT


IF(test_val <=TOL) THEN
iter_flag =0

ELSEIF(iter <lnit) THEN

GOTO 100

ELSE	! iter =lnit, and test_val >TOL
iter_flag =1

ENDIF

!*
!write(1,*) 'Colin:'
!do ie =1,ne
!write(1,*) colin(ie)
!enddo
!*

! Check:
!write(1,*) 'WH*W:'
!write(1,*) matmul(transpose(conjg(W)), W);

end subroutine Sim_Iter

