! Main for Shifted Inverse Power method
! Real matrix A: Real eigenvalues
! Input File:
! - Title
! - Matrix order, kodA, kodw0
! - Matrix A (on rows)
! - Vector w0
! - TOL, lnit
! - Number of eigenvalues (aproximations) "ns"
! - Eigenvalues "s": "ns" values (s =shift: A -> A-bI),
! - kod_met
! - kodTest, kodPrint

! kodA =0: read A; #0: generate A by subruutine Agen(A)
! kodw0 =0: read w0; #0: generate w0 by subroutine W0gen(w0)
!
! TOL = Tolerance  on "lambda" or for z||z0; lnit = Limit number of iterations.
!
! kod_met =0/#0: iterate with B=A^(-1) / use LU solver.
! kodTest =0/#0: Test on |lambda -lambda0| / Test on z || z0.
! kodPrint =0/1/>=2: Do Print: Nothing/Eigenvectors' check/Iterations & Eigenvectors' check
!

Program PowerShiftMethod
use work

character(260) fname, titlu *80
character(10) ch; logical do_check;

integer width
real az, lz
integer status


call GetFile(fname, status, 'Select Input file'C)
open(1, file =fname, status ='old', form ='formatted', iostat =io)

Err_If: if(io /=0) then
write(*, '(/,a,/)') ' Error in opening Input File !'
goto 600
endif Err_If

read(1, '(a)') titlu
read(1,*) n, kodA, kodw0

! Allocate:
allocate(A(n,n), A0(n,n), w0(n), z0(n), z(n))

A =0.; A0 =0.;
w0 =1.;
z0 =0.; z=0.

! Read/Generate:
IF(kodA ==0) THEN
do i =1,n
read(1,*) (a0(i,j), j =1,n)
enddo

ELSE	! Generate:
call agen(A0)
ENDIF

IF(kodw0 ==0) THEN
read(1,*) (w0(i), i =1,n)

ELSE	! Generate:
call w0gen()
ENDIF

read(1,*) TOL, lnit

read(1,*) ns
allocate(shift(ns))

read(1,*) shift

read(1,*) kod_met
read(1,*) kodTest, kodPrint

read(1,'(a)') ch
do_check =ch(1:1) =='C' .or. ch(1:1) =='c' .or. ch(1:1) =='P' .or. ch(1:1) =='p'

allocate(val(n), vec(n,n))
val =0.
vec =0.

! Build A-bI:
S_DO: DO is =1, ns
A =A0

s =shift(is)
if(s /=0) then
do i=1,n
a(i,i) =a(i,i) -s
enddo
endif

iwis =width(is)
if(kodPrint ==0) then
write(1,'(//, a, I<iwis>, a, \)') ' Shift #', is, ':'
else
write(1,'(//, a, I<iwis>)') ' Shift #', is
endif


if(kod_met ==0) then
call Inverse_Shift(iter, kodTest, kodPrint)
else
call LU_Shift(iter, kodTest, kodPrint)
endif

val(is) =1./lambda +s

if(kodPrint ==0) then
iw =width(iter)
write(1, '(a, I<iw>)') ' Iteration No. ', iter
endif

write(1, '(/, a, g15.7)') ' Eigenvalue:', val(is)

write(1, '(/, a)')  ' Eigenvector:'

! z is Normalized:
vec(:,is) =z

do i=1,n
iwi =width(i)
write(1,'(a, I<iwi>, a, g15.7)') ' (', i, ')', z(i)
enddo

ENDDO S_DO

IF(do_check) THEN
write(1,'(//, a)') ' Check: A*z -lambda*z ~ 0'
if(kodPrint >=1) write(1, '(/,a,\)') ' Eigenvectors:'

dif_max =0; is_max =0;

DO is =1,ns
iwis =width(is)+1
if(kodPrint >=1) write(1, '(/, I<iwis>)') is

lambda =val(is)
z =vec(:,is)

dif =0.
idif =0

do i=1,n
az =dot_product(a0(i,:), z)
lz =lambda*z(i)
!write(1,*) az ,lz

dif_i =az -lz

	if(kodPrint >=1) then
	write(1,*) dif_i
	endif

if(abs(dif_i) >abs(dif)) then
dif =dif_i
idif =i
endif

if(abs(dif) >abs(dif_max)) then
dif_max =dif
idif_max =idif
is_max =is
endif
enddo

!enddo


!write(1, '(/, a, 1pg10.3)') ' Maximum difference (in modulus) =', dif

!if(idif >0) then
!iwi =width(idif)
!write(1, '(a, I<iwi>)') ' Coordinate #', idif
!endif

ENDDO

write(1, '(/, a, (1pg10.3, :, 3x, 1pg10.3))') ' Maximum difference (in modulus) =', dif_max
iws =width(is_max); iwi =width(idif_max);
write(1,'(a, I<iws>, a, I<iwi>, a)') ' Eigenvector # ', is_max,';  Coordinate # ', idif_max, '.'

ENDIF

600 continue

END program PowerShiftMethod



