! DUBLA Precizie
! Val. & Vect. proprii ale unei problemei generalizate (A-lambda*B)x =0
! Cu: A, B = SIMETRICE & pozitiv definite.
! Valori si vectori proprii = reali
!
! A =K; B =M;
!
! Intrare - matrici A, B:
! - Citire din fiser: TRIUNGHIUL SUPERIOR - pe linii
! - Generare cu subrutina A_gen, B_gen: Matricea COMPLETA (n,n)
!
! Fisierul de intrare:
! 1) Titlul problemei
! 2) n, kod_A, kod_B, kodPrint
!
! 3_1) [job_A - pentru kod_A <=0]
! 3_2) [Matricea A - pentru kod_A =1]
!
! 4_1) [job_B - pentru kod_B <=0]
! 4_2) [Matricea B - pentru kod_B =1]
!
! 5) [A_factor: pentru kod_A >=2]
!
! 6) cod_Eigen
!
! 7) kodPrintVP
!
! 8) Proba
! 9) Do_Gen_Files (logic)
!
! Semnificatia datelor:
! 1) Titlul: text <=80 caractere
! 2) n: ordinul matricii
!	 kod_A, kod_B: cod de introducere a matricilor A, B.
!		kod_A, kod_B <=0: matrice generata de A_gen, respectiv B_gen (matricea completa)
!		kod_A, kod_B  =1: matrice citita din fisierul de intrare  (triungiul superior, pe linii)
!		kod_A, kod_B >=2: matrice citita din fisier specificat  (triungiul superior, pe linii)
!	kodPrint: cod tiparire matrici A, B: 0/#0
!       =0: nu se tiparesc A si B
!		=1: tipareste A
!		=2: tipareste B
!	   >=3: tipareste A si B
!
! 3_1, 4_1) Job_A, job_B: text <=20 caractere: identifica secventa de generare din Agen, Bgen.
! 3_2, 4_2) Matricile A, B.
!
! 5) A_factor: factor de multiplicare pentru matricea A: numai pt. matrice citita din fisier specificat (Kod_A >=2)
!
! 6) cod_Eigen: text 10 caractere; primul caracter este semnificativ:
!     p = pulsatii
!	  f = frecvente
!	  caracter #p,f: lambda
!
! 7) kodPrintVP: cod tiparire vectori proprii (in fisierul de date): 0/#0 = Nu/Da
!
! 8) 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)
! 9) Genereaza sau nu (T/F): Fisierul de valori proprii & Fisierul de vectori proprii.

use numerical_libraries


character(255) fname, fname2, titlu *80
real(8), allocatable :: a(:,:), b(:,:), Q(:,:)
real(8), allocatable :: eval(:), eval_pf(:), evec(:,:)

real(8) lambda
real(8), allocatable :: y1(:), y2(:)

character(1) w
logical proba, full_proba

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

character(20) job_A, job_B
character(80) dummy

logical Do_Gen_Files

integer width

logical do_PrintA, do_PrintB

character(10) cod_Eigen
logical pulsatii, frecv

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

! Retrieve working directory:
call GetFile(fname, status, 'Selecteaza Fisierul de Intrare'C)

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

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

read (1,*) n, kod_A, kod_B, kodPrint

allocate(a(n,n), b(n,n))
allocate (eval(n), evec(n,n))

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

kod_A_IF: IF(kod_A <=0) THEN
read(1,*) job_A

call A_gen(n, A, job_A)

ELSE	! Read A: kod_A >=1

	if(kod_A ==1) then
	iunit =1

	else
	iunit =2
	call GetFile(fname2, status, 'Selecteaza Fisierul matricii A  (matricea K)'C)

	open(2,file =fname2, status ='old', form ='formatted')
	read(2,*) dummy
	endif

do i =1,n
read(iunit,*) (a(i,j), j =i,n)	! Upper Triangle
enddo

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

if(iunit ==2) close(2)

ENDIF kod_A_IF


kod_B_IF: IF(kod_B <=0) THEN
read(1,*) job_B

call B_gen(n, B, job_B)

ELSE	! Read B: kod_B >=1

	if(kod_B ==1) then
	iunit =1

	else
	iunit =2
	call GetFile(fname2, status, 'Selecteaza Fisierul matricii B  (matricea M)'C)

	open(2,file =fname2, status ='old', form ='formatted')
	read(2,*) dummy
	endif

do i =1,n
read(iunit,*) (b(i,j), j =i,n)	! Upper Triangle
enddo

! Fill lower triangle:
do i=2,n
	do j=1, i-1
	b(i,j) =b(j,i)
	enddo
enddo

if(iunit ==2) close(2)

ENDIF kod_B_IF

if(kod_A >=2) then

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

read(1, '(a1)') cod_Eigen

pulsatii =cod_Eigen(1:1) =='p' .or. cod_Eigen(1:1) =='P'
frecv =cod_Eigen(1:1) =='f' .or. cod_Eigen(1:1) =='F'

read(1,*) kodPrintVP

! Store:
allocate(a0(n,n), b0(n,n))
a0 =a; b0 =b;

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, B:
do_PrintA =kodPrint ==1 .or. kodPrint >=3
do_PrintB =kodPrint ==2 .or. kodPrint >=3

A_Print: IF(do_PrintA) THEN

write(1,'(/,a)') 'Matricea A:'

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

ENDIF A_Print

B_Print: IF(do_PrintB) THEN

write(1,'(/,a)') 'Matricea B:'

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

ENDIF B_Print


! Subrutina IMSL:
CALL dGVCSP (N, A, N, B, N, EVAL, EVEC, N)

if(pulsatii .or. frecv) allocate (eval_pf(n))

IF(pulsatii) THEN
write(1, '(/,/, a)') ' Pulsatii:'
write(*, '(/, a)') ' Pulsatii:'

EVAL_pf =dsqrt(EVAL)

kod_eigen =1

ELSEIF(frecv) THEN
write(1, '(/,/, a)') ' Frecvente:'
write(*, '(/, a)') ' Frecvente:'

pi_2 =2*3.1415926535897932		!384626433832795
EVAL_pf =dsqrt(EVAL)/(pi_2)

kod_eigen =2

ELSE

write(1, '(/,/, a)') ' Valori proprii:'
write(*, '(/, a)') ' Valori proprii:'

kod_eigen =0

ENDIF



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

if(pulsatii .or. frecv) then
write(1, '(a, I<iw>, a, 1pg24.16)') '(',i,')', eval_pf(i)
write(*, '(a, I<iw>, a, 1pg24.16)') '(',i,')', eval_pf(i)

else

write(1, '(a, I<iw>, a, 1pg24.16)') '(',i,')', eval(i)
write(*, '(a, I<iw>, a, 1pg24.16)') '(',i,')', eval(i)
endif

enddo


! Write each vector separately:
IF(kodPrintVP /=0) THEN
write(1, '(/, a)') ' Vectori proprii:'

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

ido: do i =1,n
write(1, '(g24.16)') evec(i,j)
enddo ido

enddo jdo

ENDIF


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

if(pulsatii .or. frecv) then
call Gen_Eigen_Files(fname, n, eval_pf, evec, kod_eigen)

else

call Gen_Eigen_Files(fname, n, eval, evec, kod_eigen)
endif

ENDIF Gen_Eigen_Files_If


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


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

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

dif =0

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

y1 =matmul(a0, evec(:,ivp))
y2 =eval(ivp)*matmul(b0, evec(:,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


write(1, '(/, a)') ' Check result:'
write(1,'(a, 1pg10.3)')  ' Maximum absolute difference (A*x -lambda*B*x) =', 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

END

