! EigenSystem REFINEMENT
!
! Input Files:
! - Matrix A
! - Eigenvalues File
! - Eigenvectors File
!
! DOUBLE Precision.
! Eigensystem is printed in 'g24.16' format.
!
! Main - Refine
! Symmetric or Not_Symmetric matrices
!
! Data File (Input):
! -----------------
! - Title
! - n, kodA, kodPrintA
! If kodA <0:
!	- job: parameter for "agen";
!
! If kodA =0, 1
!   - Matrix A (on rows)
!
! [If kodA >=2: matrix read from file "Fname_2". No data input. Fname_2 is selected via the standard dialog window.]
!
! - real_val (T/F)
! - Eps_lambda, Lnit, kodIter
! - kod_PrintVP
! - is (integer)
! - kodNorm
! - Check_text (C,c,P,p, F,f; or, other chars)
! - Do_Gen_Files (logical:T,F)
!
! Data meanning:
! -------------
! Title: text data (1-80 chars). Problem description.
!
! n: matrix order
! kodA >=0: matrix read from data file; kodA  <0: generated matrix (by Agen);
!	Do Input:
!	kodA =0:  upper triangle;
!	kodA =1:  full input
!
!   kodA >=2: matrix read from file "fname_2", namely:
!		kodA =2:  upper triangle.
!		kodA =3:  full matrix
!		kodA >=4: non-zero matrix elements: i, j, a(i,j).
!	 
!	Note: "Fname_2" structure:
!         - The first line is a text-line (max. 80 chars; e.g.: problem title).
!         - Matrix entries, according to kodA.
!
! kodA<0: job: Character(20): Parameter to be transmitted to "agen".
!
! kodPrintA:
!		>0: print A (given & block-tri reduced)
!		=0: print A block-tri reduced
!		<0: print nothing
!
! Eps_lambda = Tolerance to two succesive lambdas
!
! Lnit = Limiting number of QR iterations
!
! kodIter =0/#0: do / do not: print iterations
!
! kod_PrintVP =0/#0:
!        kod_PrintVP #0: Do Write eigenvectors in the input file
!        kod_PrintVP <=0: Do Not Write eigenvectors
!
! Note: eigenvectors are already normalized by the 2-norm.
!
! is: Code that defines the shift for Inverse Power:
!			is >0: shift = lambda -one unit in the "is"th significant digit of lambda.
!			Recommendation: is =4; 5,6.
!			is <0: shift =lambda -1D-4
!
! kodNorm =Norm code:
!		<0: do not normalize
!		=0: maximum norm (Infinity-Norm)
!		=1: 1-norm
!		=2: euclidian norm (2-Norm)
!		=3: "eng. norm": First coordinate =1.
!
! Check_text: text data (1-20 chars). First character is significant:
!		- "C","c","P","p", "F", "f"/Other chars: with/without check on eigensystem;
!		  If "F" or "f": Full printing: Differencies (A*z -lambda*z) are printed, too.
!
! Do_Gen_Files: Generate or no (T/F): Eigenvalues File & Eigenvectors File.
!
!
! Data Files:
! -----------
! - stored in "Dat" Folder
! - selected via OpenFileName
!
Program Refine_Eigen
use work
use dflib

character($maxpath) fname

character(80) titlu
character($maxpath) fname_2, title_2*80

character($maxpath) fname_VAL, fname_VEC

character(20) shift, check_text, ch1*1
logical do_check, full_print_check

integer width

logical symm

logical real_val

! In "work":
!real(8), allocatable :: VEC(:,:)
!real(8), allocatable :: y(:)

!complex(8), allocatable :: VEC_Z(:,:)
!complex(8), allocatable :: y_z(:)     ! in "work"

real(8) value; complex(8) value_z;

real(8) s, lambda1
real(8), allocatable :: as(:,:), w0(:)

character(20) job

integer status
real(8) ulp

logical Do_Gen_Files

character(80) dummy



Do_Gen_Files =.FALSE.

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

if(status <=0) then
write(*,*) 'Error in selecting file !'
GOTO 600
endif

! Open:
open(1, file =fname, status ='old', form ='formatted')

! READ:
read(1, '(a)') titlu
read(1,*) n, kodA, kodPrintA

! Allocations:
allocate(A(n,n))

A =0.

! Read/Generate:
kodA_Select: SELECT CASE(kodA)

CASE(:-1)		! kodA <0: Generate
job(1:20) =' '
read(1,*) job
call agen(A, job)

CASE(0)			! kodA =0: Read from input file
do i =1,n
read(1,*) (a(i,j), j =i,n)	! upper triangle on rows
enddo
! Fill symmetric:
do i=2,n
do j=1, i-1
a(i,j) =a(j,i)
enddo
enddo

CASE(1:)		! kodA >0: Read from file
	if(kodA ==1) then	! Input File
	do i =1,n
	read(1,*) (a(i,j), j =1,n)	! full input on rows
	enddo

	else	! Specified File "fname_2"
	call GetFile(fname_2, status, 'Select Matrix File'C)

		if(status <=0) then
		write(*,*) 'Error in selecting file !'
		GOTO 600
		endif

	open(2, file =fname_2, status ='old', form ='formatted')
	read(2,*) title_2

		if(kodA ==2) then	! read upper triangle
		do i =1,n
		read(2,*) (a(i,j), j =i,n)	! upper triangle on rows
		enddo

		! Fill symmetric:
		do i=2,n
		do j=1, i-1
		a(i,j) =a(j,i)
		enddo
		enddo

		elseif(kodA ==3) then	! full input
		do i =1,n
		read(1,*) (a(i,j), j =1,n)	! full input on rows
		enddo

		else	! kodA.=4

		40 read(2,*, err =50, end =50) i, j, a(i,j)		! non-zero entries
		goto 40
		endif

	close(2)
	endif

END SELECT kodA_Select

50 continue
! Check for symmetric matrix:
if(kodA ==0) then
simetric =.true.

else
simetric =symm(n,a)
endif

read(1,*) real_val
read(1,*) eps_lambda, lnit, kodIter

! Store
lnit0 =lnit


read(1,*) kod_PrintVP

read(1,*) is

read(1,*) kodNorm

read(1,*) check_text
ch1 =check_text(1:1)
do_check = ch1 =='C' .or. ch1 =='c' .or. ch1 =='P' .or. ch1 =='p' &
			.or. ch1 =='F' .or. ch1 =='f'
full_print_check =ch1 =='F' .or. ch1 =='f'


! Generate Eigen Files
read(1,*, err =130) Do_Gen_Files

130 continue

! For printing:
n8 =min(n,8)

! Printing:
Print_A1: IF(kodPrintA >0) THEN

if(simetric) then
write(1, '(/,a)') ' Matrix A (Symmetric)'
else
write(1, '(/,a)') ' Matrix A (Not Symmetric)'
endif

call PrintA(n,A)

ENDIF Print_A1


! Save A - for eigenvectors check & shift:
allocate(a0(n,n))
a0 =a

!if(real_val) then
allocate(lambda(n))
allocate(lambda_z(n))

lambda =0
lambda_z =0

! Select Eigenvalues File & eigenvectors file:
call GetFile(fname_Val, status, 'Select EigenValues File'C)
call GetFile(fname_Vec, status, 'Select EigenVectors File'C)


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

! Read lambdas:
! Lambda:
read(2,*) dummy
do ii =1,n
	if(real_val) then
	read(2,*) i, lambda(i)
	else
	read(2,*) i, lambda_z(i)
	endif

enddo

close(2)

! Write eigenvalues:
write(1, '(/,a)') ' Eigenvalues (given)'
write(*, '(/,a)') ' Eigenvalues (given)'

IF(real_val) THEN	

do i=1,n
iw =width(i)+1

write(1, '(I<iw>, 1pG24.16)') i, lambda(i)
write(*, '(I<iw>, 1pG20.10)') i, lambda(i)

enddo

ELSE
! Write lambda_z
do i =1,n
iw =width(i)+1

if(dimag(lambda_z(i)) /=0) then

write(1, '(I<iw>, 1pg24.16, 1pg24.16)') i, lambda_z(i)
write(*, '(I<iw>, 1pg20.10, 1pg20.10)') i, lambda_z(i)

else
write(1, '(I<iw>, 1pG24.16)') i, real(lambda_z(i))
write(*, '(I<iw>, 1pG20.10)') i, real(lambda_z(i))
endif

enddo

ENDIF



REAL_VAL_IF: IF( real_val ) THEN
!lambda =real(lambda_z)

! Check for equal eigenvalues:
allocate(multiple(n), done(n))
multiple =1

ir_do: do i =1,n
value =lambda(i)

	kr_do: do k =1,n
	if(k ==i) cycle kr_do

	if( abs(lambda(k) -value) <=EPS_lambda ) then
	multiple(i) =multiple(i) +1
	endif
	enddo kr_do

enddo ir_do


max_multiple =maxval(multiple)

! Write Multiple:
Write_Multiple: IF(max_multiple >1) THEN
write(1, '(/, a)') ' Multiplicity (algebraic)'

done =0

ir2_do: do i =1,n
value =lambda(i)

if(done(i) >1) then
cycle ir2_do

else
m_i =multiple(i)


	if(m_i >1) then
	iwi =width(i)+1; iwmi =width(m_i);

	write(1, '(I<iwi>, 1pg24.16, a, I<iwmi>)'), i, value, ': ', m_i
	endif

done(i) =1

	kr2_do: do k =i+1,n

	if( abs(lambda(k) -value) <=EPS_lambda ) then
	done(k) =2
	endif
	enddo kr2_do

endif

enddo ir2_do

ENDIF Write_Multiple


allocate(t(n,n))
t =a

allocate(y(n))
y =0


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

! Read VEC:
! Vec_y:
read(3,*) dummy
do ii =1,n
read(3,*) i

do j =1,n
read(3,*) Vec(i, j)
enddo

enddo

close(3)


allocate(iter_refine(n))
iter_refine =0

allocate(p(n))
p =0.	! To store test_val
! Restore
lnit =lnit0

if(n >50) write(*, '(/, a)') ' * Refining ...'
write(1, '(//, a)') ' * Refinement:'

!!TOL =1D-15
!!TOL =2.220446049250313E-016
!TOL =5D-16

allocate(TOL(n))

do ivp =1,n
TOL(ivp) =ulp(lambda(ivp))

! Skip if lambda =0:
if(lambda(ivp) ==0) cycle

call Power_Refine(ivp)
enddo


iwlnit =width(lnit)+2; iwn =width(n)+1;
it1 =iwn+4; it2 =it1+16; it3 =it1+12; it4 =it3+23;
write(1,'(t<it1>, a, t<it2>, a, t<it4>, a)') 'Iterations', 'Test Value', ' Tolerance'

do ivp =1,n
if(iter_refine(ivp) >=0) then
write(1,'(I<iwn>, I<iwlnit>, t<it3>, 1pg20.10, t<it4>, 1pg10.3)') ivp, iter_Refine(ivp), p(ivp), TOL(ivp)
else
write(1,'(I<iwn>, I<iwlnit>, t<it3>, 1pg20.10, a)') ivp, -iter_Refine(ivp), p(ivp), '(stationary)'
endif
enddo


! Write eigenvalues:
write(1, '(/,a)') ' Eigenvalues'
write(*, '(/,a)') ' Eigenvalues'

do ivp =1,n
iwi =width(ivp)

write(1, '(1x, I<iwi>, 1pg24.16)') ivp, lambda(ivp)
write(*, '(1x, I<iwi>, 1pg20.10)') ivp, lambda(ivp)
enddo



IF(kodNorm >=0) THEN
call Normalize_r(VEC, kodNorm)
ENDIF

n6 =min(n,6)

! Write eigenvectors:
write_vectors: IF(kod_PrintVP >0) THEN
write(1, '(/, a, \)') ' Eigenvectors'

ii_DO: DO ii =1,n,6
n1 =ii
n2 =ii+5
if(n2>n) n2 =n

iw2 =width(n2)
write(1,'(/, I11,\)') n1
write(1,'(<n6-1>(23x,I<iw2>))') (i, i =n1+1, n2)

do i =1,n
write(1, '( <n6>(g24.16) )') VEC(i, n1:n2)
enddo

ENDDO ii_DO
ENDIF write_vectors
!kodCheck =0


ELSE		! REAL_VAL_IF		! Work in Complex. With real lambdas, too.
allocate(VEC_Z(n,n))
vec_z=(0.,0.)

allocate(t_z(n,n))
t_z =a

allocate(y_z(n))
y_z =0

! Vec_Z:
read(3,*) dummy
do ii =1,n
read(3,*) i

do j =1,n
read(3,*) Vec_Z(i, j)
enddo

enddo

close(3)


! Establish multiplicity:
allocate(multiple(n), done(n))
multiple =1

i_do: do i =1,n
value_z =lambda_z(i)

	k_do: do k =1,n
	if(k ==i) cycle k_do

	if( abs(lambda_z(k) -value_z) <=EPS_lambda) then
	multiple(i) =multiple(i) +1
	endif
	enddo k_do

enddo i_do

max_multiple =maxval(multiple)

! Write Multiple:
Write_Multiple_Z: IF(max_multiple >1) THEN
write(1, '(/, a)') ' Multiplicity (algebraic)'

done =0

i2_do: do i=1,n
value_z =lambda_z(i)

if(done(i) >1) then
cycle i2_do

else
m_i =multiple(i)
iwmi =width(m_i)

if(m_i >1) then
iwi =width(i)+1; iwmi =width(m_i);

	if(dimag(value_z) ==0) then
	write(1, '(I<iwi>, 1pg20.10, a, I<iwmi>)') i, real(value_z), ': ', m_i

	else
	write(1,'(I<iwi>, 1pg20.10, 1pg20.10, a, I<iwmi>)') value_z, ': ', m_i
	endif
endif

done(i) =1

	k2_do: do k =i+1,n

	if( abs(lambda_z(k) -value_z) <=EPS_lambda ) then
	done(k) =2
	endif
	enddo k2_do
endif

enddo i2_do

ENDIF Write_Multiple_Z



! Refine:
allocate(iter_refine(n))
iter_refine =0
p =0.	! To store test_val
! Restore
lnit =lnit0

if(n >50) write(*, '(/, a)') ' * Refining ...'
write(1, '(//, a)') ' * Refinement:'

!TOL =2.220446049250313D-016
!TOL =5D-16

!write(1, '(/, a, 1pg20.10)') ' TOL =', tol
allocate(TOL(n))

do ivp =1,n
TOL(ivp) =ulp(abs(lambda_z(ivp)))

! Skip if lambda =0:
if(lambda_z(ivp) ==0) cycle

call Power_Refine_Z(ivp)
enddo


iwlnit =width(lnit)+2; iwn =width(n)+1;
it1 =iwn+4; it2 =it1+16; it3 =it1+12; it4 =it3+24;
write(1,'(t<it1>, a, t<it2>, a, t<it4>, a)') 'Iterations', 'Test Value', ' Tolerance'

do ivp =1,n
if(iter_refine(ivp) >=0) then
write(1,'(I<iwn>, I<iwlnit>, t<it3>, 1pg20.10, t<it4>, 1pg10.3)') ivp, iter_Refine(ivp), p(ivp), TOL(ivp)
else
write(1,'(I<iwn>, I<iwlnit>, t<it3>, 1pg20.10, a)') ivp, -iter_Refine(ivp), p(ivp), '(stationary)'
endif
enddo


! Write eigenvalues:
write(1, '(/,a)') ' Eigenvalues'
write(*, '(/,a)') ' Eigenvalues'

do ivp =1,n
iwi =width(ivp)
if(imag(lambda_z(ivp)) /=0) then
write(1, '(1x, I<iwi>, 1pg24.16, 1pg24.16)') ivp, lambda_z(ivp)
write(*, '(1x, I<iwi>, 1pg20.10, 1pg20.10)') ivp, lambda_z(ivp)

else
write(1, '(1x, I<iwi>, 1pg24.16, 1pg24.16)') ivp, real(lambda_z(ivp))
write(*, '(1x, I<iwi>, 1pg20.10, 1pg20.10)') ivp, real(lambda_z(ivp))
endif

enddo



IF(kodNorm >=0) THEN
call Normalize_Z(VEC_Z, kodNorm)
ENDIF

! Write eigenvectors:
write_vectors_Z: IF(kod_PrintVP >0) THEN
write(1, '(/, a, \)') ' Eigenvectors'

! Real lambdas: Write groups of six vectors
	Real_If: IF(real_val) THEN

ii2_DO: DO ii =1,n,6
n1 =ii
n2 =ii+5
if(n2>n) n2 =n

iw2 =width(n2)
write(1,'(/, I11,\)') n1
write(1,'(<n6-1>(23x,I<iw2>))') (i, i =n1+1, n2)

do i =1,n
write(1, '( <n6>(g24.16) )') VEC(i,n1:n2)
enddo

ENDDO ii2_DO

	ELSE	! Complex: write vectors one by one
wr_ivp_do: do ivp =1,n
write(1, '(/, I8)') ivp

do l=1,n
	if(imag(vec_z(l, ivp)) /=0) then
	write(1,'(g24.16, g24.16)') vec_z(l, ivp)

	else
	write(1,'(g24.16)') real(vec_z(l,ivp))
	endif
enddo

enddo wr_ivp_do

	ENDIF Real_If
ENDIF write_vectors_Z

!kodCheck =1

ENDIF REAL_VAL_IF




! Generate eigen_values file & eigen_vectors file:
Gen_Eigen_Files_If: IF(Do_Gen_Files) THEN

call Gen_Eigen_Files(fname, real_val)
ENDIF Gen_Eigen_Files_If



! CHECK:
CHECK_IF: IF(do_check) THEN

call Check(real_val, full_print_check)
ENDIF CHECK_IF


500 continue


if(allocated(A)) deallocate(A)
if(allocated(a0)) deallocate(a0)

if(allocated(t)) deallocate(t)
if(allocated(t_z)) deallocate(t_z)


if(allocated(y)) deallocate(y)
if(allocated(y_z)) deallocate(y_z)

if(allocated(VEC)) deallocate(VEC)
if(allocated(VEC_Z)) deallocate(VEC_Z)

if(allocated(iter_refine)) deallocate(iter_refine)

600 continue

end Program Refine_Eigen



