! Main - for MultiVector Power Iteration
! Complex Hermitian matrix
! ne <=n Eigenvectors required
! Data File:
! - Title
! - n, ne, kodA, kodW0
! - Matrix A (on rows) - if kodA =0
! - Initial Vector_Matrix W0	- if kodW0 =0
! - TOL, lnit
! - kodTest, kodPrint
! - check_text
!
! TOL: Tolerance for "lambda" - or for ||z-z0||
! lnit = Limit number of iterations
!
! kodA:
!		<0: read upper triangle of A;
!		=0: read full A; 
!		#0: generate A bt subroutine Agen(A)
! kodw0 =0: read W0; #0: generate W0 by subroutine W0gen(W0)
!
! kodTest =0 /#0: Test on "lambda"/ on z||z0
! kodPrint =0/1/2: Do not print/Print Eigenvectors check/ Print iterations.& Eigenvectors check
!
! check_test: text <=20 chars. First char is significant:
!			C, c, P, p / other char: do / do not - test eigensystem.

Program Multivector_Power_Method
use work

character(260) fname, titlu *80, check*20
integer width

logical do_check

logical IsHermitian, IsSimetric

complex, allocatable :: z(:)
complex lambda1, az, lz
complex dif_i, dif, dif_max

logical real_matrix, real_values

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

read(1, '(a)') titlu
read(1,*) n, ne, kodA, kodw0
if(ne ==0) ne =n

! Allocations:
allocate(A(n,n), W0(n,ne), W(n,ne))
allocate(lambda(ne))

A =0.
w0 =1.
w =0.

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

! Fill hermitian:
do i=1,n
do j=1,i-1
a(i,j) =conjg(a(j,i))
enddo
enddo

ELSEIF(kodA ==0) THEN	! Read full
do i =1,n
read(1,*) (a(i,j), j =1,n)
enddo

ELSE	! Generate:
call agen(A)
ENDIF


IF(kodw0 ==0) THEN
do i=1,ne
read(1,*) (w0(i,j), j =1,n)
enddo

ELSE	! Generate:
call W0gen(w0)
ENDIF

read(1,*) TOL, lnit
read(1,*) kodTest, kodPrint

if(kodTest ==0) then
allocate(lambda0(ne))
endif

read(1,*) check
do_check =check(1:1) =='C' .or. check(1:1) =='c' .or. check(1:1) =='P' .or. check(1:1) =='p'


real_matrix =.true.
i_do: do i=1,n
do j=1,n
	if(imag(a(i,j)) /=0) then
	real_matrix =.false.
	exit i_do
	endif

enddo
enddo i_do

write(1,*)

if(IsSimetric(n,a)) then
	if(real_matrix) then
	write(1,'(a)') ' Matrice simetrica'
	else
	write(1,'(a)') ' Matrice simetrica (complexa)'
	endif
elseif(IsHermitian(n,a)) then
write(1,'(a)') ' Matrice hermitiana'

else
write(1,'(a)') ' Matrice generala'
endif

! Print A:
IF(kodA /=0 .and. kodPrint /=0) THEN

n8 =min(n,8)
do i=1,n
	if(real_matrix) then
	write(1, '(<n8>g15.7)') real(a(i,:))
	else
	write(1, '(<n8>(g15.7,g15.7))') real(a(i,:))
	endif
if(n >8) write(1,*)
enddo

ENDIF


iter =0
!call Sim_Iter(herm, iter, kodTest, kodPrint, iter_flag, test_val)
call Sim_Iter(iter, kodTest, kodPrint, iter_flag, test_val)

iw =width(iter)
IF(iter_flag ==0) THEN
	if(kodTest ==0) then
	write(1, '(/, a, I<iw>, a)') ' Iteration No. ', iter, ': Tolerance on "lambda" met.'
	else
	write(1, '(/, a, I<iw>, a)') ' Iteration No. ', iter, ': Tolerance on "W" met.'
	endif

ELSE
	write(1,'(/,a, I<iw>, a)') ' * Tolerance not met in ', iter, ' iterations'
ENDIF

write(1,'(a, 1pg10.3)') ' Test Value is:', test_val


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

new =width(ne)+1
do je =1,ne
	if(imag(lambda(je)) ==0.) then
	write(1, '(I<new>, g15.7)') je, real(lambda(je))
	else
	write(1, '(I<new>, g15.7, g15.7)') je, lambda(je)
	endif
enddo

write(1, '(/,a,\)') ' Eigenvectors:'

real_values =.true.
do je =1,ne
	if(imag(lambda(je)) /=0.) then
	real_values =.false.
	exit
	endif
enddo


IF(real_matrix .and. real_values) THEN
n4 =min(ne,4)

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

write(1, '(/, <n4>(I3, 13x))') (jw, jw =J,NJ)

do i =1,n
write(1, '(<n4>(g15.7, 1x) )') (real(W(i,jw)), jw =J,NJ)
enddo
ENDDO

ELSE
n2 =min(ne,2)

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

write(1, '(/, <n2>(I3, 28x))') (jw, jw =J,NJ)

do i =1,n
write(1, '(<n2>(g15.7, g15.7, 1x) )') (W(i,jw), jw =J,NJ)
enddo

ENDDO

ENDIF

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

allocate (z(n))
z =(0.,0.)

ivp_max =0
dif_max =0.
i_max =0

IVP_DO: DO ivp =1,ne
iwi =width(ivp) +1
if(kodPrint >=1) write(1, '(/, I<iwi>)') ivp

lambda1 =lambda(ivp)
z =W(:,ivp)

dif =0.
idif =0

DO i =1,n

az =dot_product(conjg(a(i,:)), z)
lz =lambda1*z(i)
dif_i =az -lz

if(kodPrint >=1) then
if(imag(dif_i) /=0) then
write(1,*) dif_i		!az, lz
else
write(1,*) real(dif_i)
endif
endif

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

ENDDO


if(abs(dif) >abs(dif_max)) then
dif_max =dif
ivp_max =ivp
i_max =idif
endif

ENDDO IVP_DO

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

if(ivp_max >0) then
iwp =width(ivp_max); iwp =width(i_max);

write(1, '(a, I<iwp>, a, I<iwi>, a)') ' Eigenvector #', ivp_max, ';  Coordinate #', i_max, '.'
endif

ENDIF Check_If

end Program Multivector_Power_Method




