! DUBLA Precizie
! Val. & Vect. proprii ale unei matrici A = sub forma Generala/Hessenberg
!
! Fisierul de intrare:
! 1) Titlul problemei
! 2) n, kodA, kodPrintA
! 3) Tip matrice
! 4) [job - pentru kodA <=0]
! 5) [Matricea A - pentru kodA =0,1]
! 6) kod_write_VP
! 7) Proba
! 8) Do_Gen_Files (logic)
!
! Semnificatia datelor:
! 1) Titlul: text <=80 caractere
! 2) n: ordinul matricii
!	kodA: cod de introducere a matricii A:
!   kodA <0: matrice generata (de Agen)
!   kodA >=0: matrice citita din fisiser:
!
!		kodA =0,1: citita din fisierul de intrare:
!			- kodA =0: triunghiul superior;
!			- kodA =1: intreaga matrice
!
!	kodA >=2: matrice citita din fisier specificat = "fname_2":
!			kodA =2:  triunghiul superior
!			kodA =3:  intreaga matrice
!			kodA >=4: elementele non-zero ale matricii: i, j, a(i,j)
!
!		Nota: Structura "Fname_2":
!         - Prima linie este o linie text (max. 80 caractere; ex.: titlul problemei)
!         - Elementele matricii - cf. kodA.

! 3) Tip: text, primul caracter este semnificativ
!		"H", "h": Hessenberg
!		"S", "s": simetrica
!		alt caracter: Generala
! 4) Job: text <=20 caractere: identifica secventa de generare din Agen
! 5) Matricea A:
!    Tip ="S, s": triungiul superior, pe linii
!         Altfel: matricea A completa, pe linii
! 6) kod_write_VP: 
!       >0 / <=0: Scrie / Nu scrie - vectorii proprii in fisierul de date (ca iesire)
! 7) Proba: text, primul caracter este semnificativ
!		 "C","c","P","p", "F", "f"/Alt caracter: Cu/Fara proba sistemului propriu;
!		 "F" or "f": Tiparire diferentelor (A*z -lambda*z), pentru fiecare pereche (valoare proprie-vector propriu)
! 8) Genereaza sau nu (T/F): Fisierul de valori proprii & Fisierul de vectori proprii.

Program Main_QR_IMSL
use ifwin

!include 'link_fnl_shared.h'


character(max_path) fname, fname_2
character(80) titlu, titlu_2
real(8), allocatable :: a(:,:), rQ(:,:)
real(8), allocatable :: reval(:), revec(:,:)
complex(8), allocatable :: eval(:), evec(:,:), Q(:,:), evec_conj(:,:)

character(1) tip
logical hess/.false./, sim/.false./

integer width
logical all_real

complex(8), allocatable :: z1(:), z2(:), za(:,:)
complex(8) lambda
real(8), allocatable :: y1(:), y2(:)

character(1) w
logical proba, full_proba
complex(8) dif_z

logical symm
real(8) dif_r, dif_max, dify, dif
real(8), allocatable ::a0(:,:)

character(20) job
character(80) dummy

logical Do_Gen_Files
!integer status

character(max_path) CurrentDirectory, dat_dir
integer istat

integer status
integer lenCurr

integer Open_File

character(80) w_title
real(8) ereal, eimag


istat =GetCurrentDirectory(max_path, CurrentDirectory)

! Form "dat" dir
if(istat >0) then
dat_dir =CurrentDirectory(1:istat)//'\Dat\'//char(0)
!lpstrInitialDir =dat_dir
endif


!print '(a, \)', 'Input File: '
!read *, fname

fname ="*.dat"C	!""C
w_title ='Select Input File'C
status =Open_File(fname, w_title, dat_dir)

if(status ==0) then
print '(a, \)', ' Input File: '
read *, fname
endif

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

read(1, '(a)') titlu

print '(a, a)', ' Titlu: ', trim(titlu)

read (1,*) n, kodA, kodPrintA

!if(kodA <0) read(1,*) job

read(1, '(a1)') tip
hess =tip =='H' .or. tip =='h'
sim =tip =='S' .or. tip =='s'

allocate(a(n,n))
allocate (eval(n), evec(n,n))
allocate(reval(n),revec(n,n))

a =0.
eval =0.; evec =0.;
reval =0.; revec =0.;

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

CASE(:-1)		! kodA <0: Generate
job(1:20) =' '
read(1,*) job
call agen(n, 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"

istat =GetCurrentDirectory(max_path, CurrentDirectory)

fname_2 ="*.dat"C
w_title ='Select Matrix File'C
status =Open_File(fname_2, w_title, CurrentDirectory)

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

	open(2, file =fname_2, status ='old', form ='formatted')
	read(2,*) titlu_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(2,*) (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

! cod scriere vect. proprii in fisierul de date:
read(1,*) kod_write_VP

! Store:
allocate(a0(n,n))
a0 =a

read(1,'(a1)') w
proba = w =='P' .or. w =='p' .or. w =='C' .or. w =='c' &
.or. w =='F' .or. w =='f'

full_proba = w =='F' .or. w =='f'

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

! write A
! Check for symmetric:
sim =symm(n,A)

PrintA: IF(kodPrintA >0) THEN
IF(sim .or. kodA <0) THEN

if(.not. sim) then
write(1,'(/,a)') 'Matrix A:'
else
write(1,'(/,a)') 'Matrix A (Symmetric):'
endif

do i=1,n
write(1, '(8g15.7)') a(i, :)
if(n >=16) write(1,*)
enddo

ENDIF

ENDIF PrintA

!IF(sim) THEN		! Simetrica
!call devcsf(n, a, n, reval, revec, n)	! ???
!CALL EVCSF (A, EVAL, EVEC)
!ELSE
if(hess) then	! Hessenberg
call devcrh(n, a, n, eval, evec, n) 

else		! Generala
call devcrg(n, a, n, eval, evec, n)
endif

!ENDIF

! Reals:
all_real =.true.

IF(.not. sim) THEN
do i=1,n
all_real =all_real .and. (imag(eval(i)) ==0)
enddo
ENDIF

write(*, '(/, a, L)') ' * All Real:', all_real

if(all_real) then

!if(.not. sim) allocate(reval(n),revec(n,n))
reval =real(eval)
revec =real(evec)

write(1, '(/, a,\)') ' * All Real'
endif

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

IF(sim .or. all_real) THEN
do i=1,n
iw =width(i)
write(1, '(a, I<iw>, a, 1pg24.16)') '(',i,')', reval(i)
write(*, '(a, I<iw>, a, 1pg24.16)') '(',i,')', reval(i)
enddo

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

ereal =dreal(eval(i))
eimag =dimag(eval(i))

if(eimag /=0) then
write(1, '(I<iw>, 2x, 1pg24.16, 1pg24.16)')  i, eval(i)
write(*, '(I<iw>, 2x, 1pg24.10, 1pg24.10)')  i, eval(i)

else
write(1, '(I<iw>, 2x, 1pg24.16)') i, ereal
write(*, '(I<iw>, 2x, 1pg24.10)') i, ereal
endif

enddo

ENDIF


! Write each vector separately:
write_VP: IF(kod_write_VP >0) THEN
write(1, '(/, a)') ' Eigenvectors:'

jdo: do j =1,n
write(1, *)
write(1,*) j

ido: do i =1,n
ereal =dreal(evec(i,j))
eimag =dimag(evec(i,j))

if(eimag /=0) then
!write(1, '(a, g24.16, a, g24.16, a)') '(', ereal, ',', eimag, ')'
write(1, '(1pg24.16, 1pg24.16)') evec(i,j)

else
write(1, '(1pg24.16)') ereal
endif
enddo ido

enddo jdo
ENDIF write_VP

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

call Gen_Eigen_Files(fname, n, eval, evec, reval, revec, all_real)
ENDIF Gen_Eigen_Files_If


! Proba - Complex:
P_IF: IF(proba) THEN

Real_IF: IF(.not. all_real) THEN

write(1,'(/,a)') ' Proba - Complex:'
allocate(z1(n), z2(n))
allocate(za(n,n))

za =cmplx(a0)

!write(1,'(a,19x,a,\)') ' A*z', 'lambda*z:'
write(1,'(a,\)') ' A*z -lambda*z'

dif =0

DO ivp =1,n
IF(full_proba) write(1, '(/, I3)') ivp

z1 =matmul(za, evec(:,ivp))
z2 =eval(ivp)*evec(:,ivp)

eimag =imag(eval(ivp))


do i=1,n

IF(full_proba) THEN
if(eimag /=0) then
write(1,*)  z1(i) -z2(i)
else
write(1,*)  real(z1(i)) -real(z2(i))
endif
ENDIF

! diferenta:
dif_z =z1(i) -z2(i)
dif_r =max( abs(dreal(dif_z)), abs(dimag(dif_z)))
if( dif_r >dif ) then
dif =dif_r;		! write(1,*) 'dif =', dif;

ivmax =ivp
icoord =i
endif
enddo

ENDDO

ELSE
write(1,'(/,a)') ' Proba - Real:'
allocate(y1(n), y2(n))

!write(1,'(a,13x,a,\)') ' A*x', 'lambda*x:'
write(1,'(a,13x,a,\)') ' A*x -lambda*x:'

dif =0

DO ivp =1,n
IF(full_proba) write(1, '(/, I3)') ivp

y1 =matmul(a0, revec(:,ivp))
y2 =reval(ivp)*revec(:,ivp)

do i=1,n
if(full_proba) write(1,*)  y1(i) - y2(i)
! Diferenta:
dify =dabs(y1(i)-y2(i))
if( dify >dif ) then
dif =dify;			!write(1,*) 'dif =', dif;

ivmax =ivp
icoord =i
endif
enddo

ENDDO

ENDIF Real_IF

write(1, '(/, a)') ' Check result:'
write(1,'(a, 1pg10.3)')  ' Maximum absolute difference (A*y -lambda*y) =', dif

if(ivmax /=0) then
iwv =width(ivmax); iwc =width(icoord);
write(1, '(a, I<iwv>, a, I<iwc>, a)') ' Eigenvector #', ivmax, '; Coordinate #', icoord, '.'
else
write(1,'(a)') ' All vectors.'
endif

ENDIF P_IF

600 continue

END

