! Main - for MultiVector Power Iteration
! Double Precision
! Complex Hermitian matrix
!
! Data File:
! - Title
! - n, ne, kodA, kodW0
! - Matrix A (on rows) - if kodA =0
! - Initial Vector_Matrix W0	- if kodW0 >0
! - tol, lnit
! - kod_writeVP
! - kodTest, kodPrint
! - check_text
!
! Data meaning:
! n : Matrix order
! ne: Number of required eigenvectors (ne <=n)
!
! kodA: Input code for A:
!		kodA <=0: generated matrix - by subroutine Agen;
!		kodA =1: inputted matrix - read from input file;
!		kodA >=2 inputted matrix - read from specified file. 
!				The first line in the specified file is a text-line (dummy).
!
! kodw0 =0: read W0; #0: generate W0 by subroutine W0gen(W0)
!
! tol : Tolerance for "lambda" - or for ||z-z0||
! lnit: Limit number of iterations
!
! kod_writeVP: <=0/>0: do/do not: werite eigenvectors in input file
!
! 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

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

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

logical do_check

logical IsHermitian, IsSimetric

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

logical real_matrix, real_values

real(8) a_factor

character(80) dummy

 CHARACTER(3)        drive
 CHARACTER(256)      dir
 CHARACTER(256)      name
 CHARACTER(256)      ext


! Input File:
call GetFile(fname, status, 'Select Input File'C)

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:
! Read/Generate A:
if (kodA <=0) then
call Agen(a)

else
	if (kodA ==1) then
	do i =1,n
	read(1,*) (a(i,j), j =i,n)
	enddo

	else
	call GetFile(fname2, status, 'Select Matrix File'C)
	open(2, file =fname2, status ='old', form ='formatted')

	read(2,*) dummy
	do i =1,n
	read(2,*) (a(i,j), j =i,n)
	enddo

	close(2)	
	endif

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

read(1,*) a_factor
A =A*a_factor
endif



IF(kodw0 ==0) THEN
call W0gen(w0)

ELSE	! Generate:

do i=1,ne
read(1,*) (w0(i,j), j =1,n)
enddo
ENDIF

read(1,*) TOL, lnit
read(1,*) kod_writeVP
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>(1pg24.16))') real(a(i,:))
	else
	write(1, '(<n8>(1pg24.16,1pg24.16))') 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>, 1pg24.16)') je, real(lambda(je))
	else
	write(1, '(I<new>, 1pg24.16, 1pg24.16)') je, lambda(je)
	endif
enddo

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

Write_VP: IF(kod_writeVP >0) THEN

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

real_if: 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, 22x))') (jw, jw =J,NJ)

do i =1,n
write(1, '(<n4>(1pg24.16, 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, 46x))') (jw, jw =J,NJ)

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

ENDDO

ENDIF real_if

ENDIF Write_VP


! Write Output Files:
result = SPLITPATHQQ (fname, drive, dir, name, ext)
ln =len_trim(name)

fname_Val =trim(name)//'-VAL.dat'
fname_Vec =trim(name)//'-VEC.dat'

open(2, file =fname_Val, status ='unknown', form ='formatted')
open(3, file =fname_Vec, status ='unknown', form ='formatted')


! Eigenvalues:
write(2, '(a)') 'Eigenvalues:'

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

! Eigenvectors:
write(3, '(/,a,\)') ' Eigenvectors:'


DO ivp =1,ne
iwi =width(ivp)+1
write(3, '(/,I<iwi>)') ivp

do i =1,n

	if(real_matrix .and. real_values) then
	write(3, '(1pg24.16)') real(W(i,ivp))

	else
	write(3, '(1pg24.16)') W(i,ivp)
	endif

enddo
ENDDO

close(3)

! 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




