! Inverse Power method - with LU
! A, w and z - defined in "work"

subroutine LU_Shift(iter, kodTest, kodPrint)
use work
real vNorm2
real lambda0	! Previous value
integer order(n)

real colin(n)
logical do_pvt


! LU decomposition:
!call LUdecomp(A, n, order, kodLU)
kod =0

! Initializeaza ORDER:
do i =1, n
order(i) =i
end do

do_pvt =.true.
kodPivot =1
call LUdecomp(A, n, order, do_pvt, kodPivot, kodLU)
if(kodLU /=0) then
print*, 'Pivot < 1e-7. Stop'
call EXIT
endif


! Normalize guess vector:
z0 =w0/vNorm2(n,w0)

! Printing:
IF(kodPrint >=2) THEN
write(1, '(/, a)') ' Guess vector:'
write(1,*) (z0(j), j=1,n)

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


lambda0 =0

! Initialize counter:
iter =0

100 iter =iter+1

! Solve Aw =z
!call LUsolve(A, z0, z, n, order)!LUsolve(LU, b, n, x, order)
call LUsolve(A, z0, n, z, order)

do l=1,n
if(isnan( z(l) )) then
write(1,*) l,z(l)
endif
enddo

! Calculate "lambda": Seek the largest coordinate of vector w:
imax =maxindex(z)
lambda =z(imax)/z0(imax)

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

! Iterations printing:
IF(kodPrint >=2) THEN
write(1,'(I10, <n>G15.7, g15.7)') iter, (z(j), j=1,n), lambda
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(0)		! Test "lambda":
test_val =abs(lambda -lambda0)

CASE DEFAULT	! 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)
test_val =vNorm2(n, colin)

END SELECT

if(test_val <=TOL) then
iter_flag =0
write(1,*) 'test_val:', test_val
RETURN

elseif(iter ==lnit) then
iter_flag =-1
write(1,*) 'test_val:', test_val
RETURN

else
z0 =z
lambda0 =lambda
GOTO 100
endif


500 continue
end subroutine LU_Shift
