! Jacobi Method
! DOUBLE PRECISION
! Real Symmetric matrix. Input = Upper triangle
! Matrix A = stored in a(n*(n+1)/2), columnwise.

! Data File:
! 1) Title
! 2) n, kodA, kodPrintA
! 3) A - upper triangle on ROWS, for kodA =1
!	3-1) A_factor: only for kodA >=1 (inputted matrix)
! 4) Tol, kodPrintJ
! 5) kodVP, kod_writeVP
! 6) check
!
! Title: text <= 80 chars: Problem description
!
! n = matrix order
!
! 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).
!
! kodPrintA = Printing code for A: 0/#0 = No/Yes
!
!	A_factor = factor to multiply the inputed matrix: kodA >=1
!
! Tol = tolerance for off-diagonals
!
! kodPrintJ = Printing code for succesive transformations of A: 0 =No; #0 =Yes
!
! kodVP = Code for eigenvectors computation: 0 =No; #0 =Yes
! kod_writeVP =<=0/>0: do/do not: write eigenvectors in input file
!
! check: Only for kodVP #0.
!		 Text <=20 chars. First character is significant:
!		 C,c,P,p, F,f / other chars: Do/Do not check eigensystem. 
!		F,f: Full-Check: writes check for each eigenvector

Program Jacobi_Method
use work
use dflib

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

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

real(8), allocatable :: a0(:), y(:), a_i(:)
real(8) lambda, ay, ly

real(8) a_factor

character(80) dummy

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



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


open(1, file =fname, status ='old', form ='formatted')
read(1, '(a)') titlu
read(1,*) n, kodA, kodPrintA

n1 =n*(n+1)/2
print_it =kodPrintA /=0

! Allocate:
allocate(a(n1))
a =0.

! Read/Generate A:
if (kodA <=0) then
call Agen(a)

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

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

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

	close(2)	
	endif



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

read(1,*) tol, kodPrintJ
read(1,*) kodVP, kod_writeVP

if(kodVP /=0) then
read(1,*) check

full_check =check(1:1) =='F' .or. check(1:1) =='f' 

do_check =check(1:1) =='C' .or. check(1:1) =='c' .or. check(1:1) =='P' .or. check(1:1) =='p' .or. full_check
endif

if(do_check) then
allocate(a0(n1))
a0 =a
endif

! Print A:
IF(print_it) THEN
write(1, '(/,a)') ' Matrix A:'

do i =1,n
j1 =i; j2 =min(i+n-1, n); nj =j2-j1+1
write(1, '(<nj>g24.16)') (a(loca(i,j)), j= j1,j2)
enddo

ENDIF


! Eigenvectors:
! Allocate & Initialize eigenvectors' matrix:
if(kodVP /=0) then
allocate(X(n,n))

! Initialize X =I:
X =0.
do i =1,n
X(i,i) =1.
enddo
endif

call Jacobi(kodVP, kodPrintJ, iter, amax)

! Print the no. of iterations:
iwi =width(iter)
write(1, '(/, a, I<iwi>, a)') ' Iteration #', iter, ':'
write(1, '(a, 1pg24.16)') ' Maximum off-diagonal:', amax

! Print A:
IF(print_it) THEN

write(1, '(/,a)') ' Matrix A - diagonalized:'

do i =1,n
j1 =i; j2 =min(i+n-1, n); nj =j2-j1+1
write(1, '(<nj>g24.16)') (a(loca(i,j)), j= j1,j2)
enddo

ENDIF

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

do i=1,n
iw =width(i)+1
k =loca(i,i)		!k =(i-1)*i/2
write(1, '(I<iw>, a, 1pg24.16)') i, ')', a(k)
write(*, '(I<iw>, a, 1pg24.16)') i, ':', a(k)
enddo

! Write Eigenvectors:
IF(kodVP /=0 .and. kod_writeVP >0) THEN

if(n<=6) then		! Matrix
write(1, '(/, a)') ' Eigenvectors'' Matrix'
n6 =min(n,6)
do i=1,n
write(1, '( <n6>(g24.16, 1x) )') (X(i,j), j =1,n)
enddo

else	! Individual vectors:
write(1, '(/, a, \)') ' Eigenvectors'

do ivp =1,n
write(1, '(/, I4)') ivp
	do l=1,n
	write(1, '(g24.16)') X(l,ivp)
	enddo
enddo

endif
ENDIF


! 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')


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

do iv =1,n
iwi =width(iv)
k =loca(iv,iv)
write(2, '(1x, I<iwi>, 1pg24.16)') iv, a(k)
enddo

close(2)

write(3, '(a, \)') ' Eigenvectors'

do ivp =1,n
write(3, '(/, I8)') ivp

do l =1,n
write(3, '(g24.16)') X(l, ivp)
enddo

enddo

close(3)


! Check:
Check_If: IF(do_check) THEN
write(1,'(//, a)') ' Check: A*y =lambda*y'

allocate (y(n), a_i(n))
y =0.; a_i =0.;

ivp_max =0
dif_max =0.
i_max =0

IVP_DO: DO ivp =1,n
iwi =width(ivp)

if(full_check) write(1, '(/, a, I<iwi>)') ' Eigenvector #', ivp

ki =loca(ivp,ivp)
lambda =a(ki)
y =X(:,ivp)

dif =0.
idif =0

DO i =1,n

! Set a_i: a_i(j) =A0(i,j)
	do j =1,n

	if(j >=i) then
	a_i(j) =a0(loca(i,j))
	else
	a_i(j) =a0(loca(j,i))
	endif

	enddo

ay =dot_product(a_i(:), y)
ly =lambda*y(i)

dif_i =ay -ly
if(full_check) write(1,*) dif_i		!ay, ly

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

ENDDO


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

ENDDO IVP_DO

write(1, '(/, a, 1pg10.3)') ' Maximum difference in modulus (Ay-lambda*y)  =', dif_max

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

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

ENDIF Check_If

if(allocated(X)) deallocate(X)

deallocate(a)
END Program Jacobi_Method
