! Main for Shift Inverse Power method
! Complex/Real matrix A
!
! Input File:
! - Title
! - Matrix order, kodA, kodw0
! - If kodA =0: Matrix A (on rows)
! - If kodw0 =0: Vector w0
! - TOL, lnit
! - ns
! - Shifts "s" (Eigenvalue approximations) - "ns" values (s =shift: A -> A-sI),
! - kod_met
! - kodTest, kodPrint
! - Check_text
!
! kodA =0: read A; #0: generate A by subroutine Agen(A)
! kodw0 =0: read w0; #0: generate w0 by subroutine w0gen(reals)
!
! TOL =Tolerance TOL on "lambda" or for z||z0
! lnit =limit number of ietrations
!
! ns =Number of eigenvalues (aproximations)
!
! kod_met =0: iterate with B=A^(-1)
! kod_met #0: use LU solver.
!
! kodTest = 0/#0: Test on |lambda -lambda0| / Test on z || z0.
! kodPrint =0/#0: Do/Do not print iterations.
!
! Check_text (10 chars; First char is significant): "C","c", "P", "p" = Do check
!
Program InversePower_Complex
use work

character(260) fname, titlu *80

character(10) check_text; logical do_check;

integer width
complex az, lz;
complex dif_zi, dif, dif_max;

complex s;
integer, allocatable :: Skip_Check(:)
logical hermitian


print '(a, \)', ' Input File: '
read *, fname
open(1, file =fname, status ='old', form ='formatted')

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

! Alocari:
allocate(A(n,n), A0(n,n))
A =0.; A0 =0.;

! Citiri/Generare:
IF(kodA ==0) THEN
do i =1,n
read(1,*) (a0(i,j), j =1,n)
enddo
ELSE
! Generate:
call agen(A0)
ENDIF

allocate(w0(n), z0(n), z(n))
w0 =(1.,0.); z0 =(1.,0.) ; z =0.


IF(kodw0 ==0) THEN
read(1,*) (w0(i), i =1,n)
ELSE
! Generate:
call w0gen()
ENDIF

read(1,*) ns

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


! Coduri:
read(1,*) TOL, lnit
read(1,*) kod_met
read(1,*) kodTest, kodPrint

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

if(do_check) then
allocate(Skip_Check(n))
Skip_Check =0
endif

! Verifica hermitiana:
if(hermitian(n,a0)) then
write(1,'(/,a)') ' Matrice hermitiana.'
else
write(1,'(/,a)') ' Matrice generala.'
endif


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


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

a =a0
s =shift(is)

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


if(kod_met ==0) then
call Inverse_Shift_Complex(iter, kodTest, kodPrint, kodLU, test_val)
else
call LU_Shift_Complex(iter, kodTest, kodPrint, kodLU, test_val)
endif

if(kodLU <0) then
	if(kodLU ==-1) then
	write(1,'(a)') ' Pivot < Threshold. Step left.'
	else
	write(1,'(a)') ' Singular matrix. Step left.'
	endif
Skip_Check(is) =1
cycle S_DO

else
val(is) =(1.,0.)/lambda +s
endif

! Write:
iwis =width(is)
write(1,'(/, a, I<iwis>, a)') ' Shift #', is, ':'

iw =width(iter)
write(1, '(a, I<iw>)') ' Iteration No. ', iter

if(kodPrint >=1) write(1,*) 'test_val:', test_val

	if(imag(val(is)) /=0) then
	write(1, '(/, a, g15.7, g15.7)') ' Eigenvalue:', val(is)

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

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

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


do i=1,n
iwi =width(i)
	if(imag(z(i)) /=0) then
	write(1,'(a, I<iwi>, a, g15.7, g15.7)') ' (', i, ')', z(i)

	else
	write(1,'(a, I<iwi>, a, g15.7)') ' (', i, ')', real(z(i))
	endif
enddo

ENDDO S_DO

IF(do_check) THEN
write(1,'(//, a)') ' Check: A*z -lambda*z'

dif_max =(0.,0.); idif_max =0;

S2_DO: DO is =1,ns

if(kodPrint >=1) then
iwis =width(is)
write(1, '(/, a, I<iwis>)') ' Eigenvector #', is
endif

if(Skip_check(is) /=0) then
write(1,*) 'Skipped.'
cycle S2_DO
endif

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

dif =(0.,0.); idif =0; 

do i=1,n
az =dot_product(conjg(a0(i,:)), z)
lz =lambda*z(i)
dif_zi =az -lz

if(kodPrint >=1) then
	if(imag(z(i)) /=0) then
	write(1,*) dif_zi
	else
	write(1,*) real(dif_zi)
	endif
endif

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

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

ENDDO S2_DO


if(imag(dif_max) /=0) then
write(1, '(/, a, (1pg10.3, :, 3x, 1pg10.3))') ' Maximum difference (in modulus) =', dif_max
else
write(1, '(/, a, (1pg10.3))') ' Maximum difference (in modulus) =', real(dif_max)
endif

iws =width(is_max); iwi =width(idif_max);
write(1,'(a, I<iws>, a, I<iwi>)') ' Eigenvector # ', is_max,';  Coordinate ', idif_max
ENDIF

END program InversePower_Complex



