! Problema Generalizata
! Calculul pulsatiilor proprii & a vectorilor proprii - cu matricea R
! Valorile proprii (lambda) & vectorii proprii (y): Fisiere generate de QR
! Matricea S^(-1) & Factorul de scalare pt. R: Fisiere generate de General_R
! omega =sqrt(lambda*scal_R)
! x =S^(-1)*y
!
! Optional:
! Cu Proba: K*x -lambda*M*x ~ 0
!
! Folderul de lucru va contine 1 fisier de intrare (v.mai jos) & urmatoarele 4 fisiere:
! 1) nume_exemplu-VAL.dat   (lambda); 
! 2) nume_exemplu-VEC.dat   (vectorii proprii y - ai matricii R);
! 3) nume_exemplu-Scal_R.dat  (factorul de scalare);
! 4) nume_exemplu-S^-1.dat  (matricea S^(-1).
! - Fisierul #1 se selecteaza.
! - Daca nu se gasesc fisierele 2-4 (cu numele de mai sus): se cere selectarea fisierelor 2-4.
!
! Fisier de date (intrare):
! ------------------------
! 1) Titlu: <=100 caractere
! 2) n: Ordinul matricilor
! 3) Pulsatii/Frecvente: T/F
! 4) Cod normare:
!  <0: nu se normeaza
!   0: norma-infinit
!	1: norma-1
!	2: norma-2
! >=3: norma-infinit
! [5) Proba: Text>= 1 char. Primul character e semnificativ: 
!     P, C, F/char = Da/Nu; F =Full_Proba.
! - Daca se cere proba: Matricile M si K se citesc dintr-un fisier exterior (se cere selectia).
!   Fisierul are structura (linii):
!   - Antet_M (<=80 caractere)
!	- matricea M: triunghiul superior
!	- m_factor (real)
!   - Antet_K (<=80 caractere; in linia imediat urmatoare lui "m_factor")
!	- matricea K: triunghiul superior
!	- k_factor (real)
!
! - Data (5) se poate introduce sau nu.
!
!  6) kod_print_MK = Cod tiparire matrici M & K:
!     <=0: NU; >=1: DA
!   Numai pentru Proba sau Full_Proba.
!   Daca lipseste: se ia kod_print_MK =0 ]
!
! Fisiere de iesire:
! 1) nume_exemplu-OMEGA.dat; sau: nume_exemplu-FRECV.dat
! 2) nume_exemplu-VEC_X.dat
! Pulsatiile/Frecventele sunt sortate crescator.
!
Program Main_Retrieve
use dfwin
use dflib

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

character($maxpath) Input_File
character($maxpath) fname_Val, fname_Vec_y, fname_Scal, fname_S_1, fname_MK
character($maxpath) fname_Omega, fname_Vec_x

character(80) dummy, Titlu*100

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

real(8), allocatable :: lambda(:), omega(:)
real(8), allocatable :: Vec_y(:,:), Vec_x(:,:), S_1(:,:)	! S^(-1)

integer width

logical do_pulsatii
character(max_path) title; integer status;

real(8), allocatable :: omega_Sort(:)
integer, allocatable :: index_omega_Sort(:)
integer i_cur(1)

character(260) pathbuf;
! Proba:
logical proba, full_proba;
character(1) w1; integer kod_print_MK; logical print_MK;
real(8), allocatable :: K(:,:), M(:,:), x1(:), y1(:), y2(:)
real(8) m_factor, k_factor
real(8) lambda1
real(8) dif, dify

pi2 =6.283185307179586476925286766559


! Select Input File:
title ='Select Input File'C
call GetFile(Input_File, status, title)

! Retrieve "VAL" File:
title ='Select Eigen_Values ("-VAL") File in Working Directory'C
call GetFile(fname_Val, status, title)

result = SPLITPATHQQ (fname_Val, drive, dir, name, ext)
ln =len_trim(name)


!fname_Val   =name(:ln-4)//'-VAL.dat'
fname_Vec_y =name(:ln-4)//'-VEC.dat'C

fname_Scal =name(:ln-4)//'-Scal_R.dat'C
fname_S_1 =name(:ln-4)//'-S^-1.dat'C

! Selected:
open(1, file =fname_Val, status ='old', form ='formatted')

! Other 3:
!#2
len2 =findfileqq(fname_Vec_y, " ", pathbuf)
if(len2 ==0) then
! Retrieve File:
title ='Select Eigen_Vector ("-VEC") File in Working Directory'C
call GetFile(fname_Vec_y, status, title)
endif

open(2, file =fname_Vec_y, status ='old', form ='formatted')

!#3
len3 =findfileqq(fname_Scal, " ", pathbuf)
if(len3 ==0) then
! Retrieve File:
title ='Select Scale ("-Scal_R") File in Working Directory'C
call GetFile(fname_Scal, status, title)
endif

open(3, file =fname_Scal, status ='old', form ='formatted')

!#4
len4 =findfileqq(fname_S_1, " ", pathbuf)
if(len4 ==0) then
! Retrieve File:
title ='Select S^-1 ("-S^-1") File in Working Directory'C
call GetFile(fname_S_1, status, title)
endif

open(4, file =fname_S_1, status ='old', form ='formatted')

! Read data:
open(7, file =Input_File, status ='old', form ='formatted')
read(7, '(a)') Titlu
read(7,*) n
read(7,*) do_pulsatii
read(7,*) kod_norm

! Read "Proba":
proba =.FALSE.
read(7,'(a1)', err =110, end =110) w1
proba = w1 =='P' .or. w1 =='p' .or. w1 =='C' .or. w1 =='c' .OR. w1 =='F' .or. w1 =='f'
full_proba = w1 =='F' .or. w1 =='f'

110 if(.not. proba) close(7)

kod_print_MK =0
if(proba) then
read(7,*, err =120, end =120) kod_print_MK
endif

120 print_MK = kod_print_MK >=1

allocate(lambda(n)); allocate(omega(n));
allocate(Vec_y(n,n)); allocate(Vec_x(n,n));

lambda =0.;
omega=0.;
Vec_y =0.;
Vec_x =0.;

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

! Read data:
! Lambda:
read(1,*) dummy
do ii =1,n
read(1,*) i, lambda(i)
enddo

close(1)

! Vec_y:
read(2,*) dummy
do jj =1,n
read(2,*) j

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

enddo

close(2)

! R_Scale:
read(3,*) dummy
read(3,*) R_factor

close(3)


! S^-1 matrix:
!S_1 =0.
read(4,*) dummy
do i=1,n
read(4,*) S_1(i,i:)
enddo
!! Fill symetric:
do i =2,n
    do j=1,i-1
    S_1(i,j) =S_1(j,i)
    enddo
enddo

! Compute Eigen_Values "omega":
VAL_DO: &
do iv =1,n
!omega(iv) =dsqrt(lambda(iv)*R_factor)
lambda(iv) =lambda(iv)*R_factor
omega(iv) =dsqrt(lambda(iv))

enddo VAL_DO

! Compute Eigen_Vectors "x":
VEC_DO: &
do ivp =1,n
Vec_x(:,ivp) =matmul(S_1, Vec_y(:,ivp))
enddo VEC_DO

! Normalize:
IF(kod_norm >=0) THEN

call Normalize(n, Vec_x, kod_norm)
ENDIF


! SORTEAZA "OMEGA" crescator:
allocate(omega_Sort(n)); allocate (index_omega_Sort(n));
omega_Sort =0; index_omega_Sort =0;

i_om_Sort =0
omega_max =maxval(omega)

do i_om =1,n
omega_curr =minval(omega)
i_cur =minloc(omega)

i_om_Sort =i_om_Sort+1
omega_Sort(i_om_Sort) =omega_curr
index_omega_Sort(i_om_Sort) =i_cur(1)

omega(i_cur) =omega_max*10_8
enddo


!*
!do i =1,n
!write(10,*) i, index_omega_Sort(i)
!enddo
!*

! Write Output Files:
if(do_pulsatii) then
fname_Omega =name(:ln-4)//'-OMEGA.dat'
else
fname_Omega =name(:ln-4)//'-FRECV.dat'
endif

fname_Vec_X =name(:ln-4)//'-VEC_X.dat'

open(1, file =fname_Omega, status ='unknown', form ='formatted')
open(2, file =fname_Vec_x, status ='unknown', form ='formatted')

write(1,*) Titlu
write(2,*) Titlu

if(do_pulsatii) then
write(1, '(a,/)') 'PULSATII PROPRII (omega):'
else
write(1, '(a,/)') 'FRECVENTE PROPRII (omega/(2*PI)):'
endif

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

if(do_pulsatii) then
write(1, '(1x, I<iwi>, 1pg24.16)') iv, omega_Sort(iv)    !omega(iv)

else
write(1, '(1x, I<iwi>, 1pg24.16)') iv, omega_Sort(iv)/pi2    !omega(iv)/pi2
endif

enddo

close(1)


write(2, '(a, \)') 'VECTORII PROPRII   '

select case(kod_norm)
case (0)
write(2, '(a)') '(Norma-infinit)'

case(1)
write(2, '(a)') '(Norma-1)'

case(2)
write(2, '(a)') '(Norma-2)'

case default    ! (3:)
write(2, '(a)') '(Norma-infinit)'

end select

! Acorda vectorii proprii cu sortarea:
do ivp =1,n
write(2, '(/, I8)') ivp

ivp_Sort =index_omega_Sort(ivp)
do l =1,n
write(2, '(g24.16)') VEC_x(l, ivp_Sort)  !VEC_x(l, ivp)
enddo

enddo

close(2)

! Proba - Complex:
Proba_IF: IF(proba) THEN
! Citeste K si M - Triunghiul superior:
title ='Select "M & K" File'C
call GetFile(fname_MK, status, title)
open (3, file =fname_MK, status ='old', form ='formatted')

allocate(K(n,n), M(n,n))
K =0D0; M=0D0;
allocate(x1(n), y1(n), y2(n))
x1 =0; y1 =0.; y2 =0.;

! M:
read(3,'(a80)') dummy   ! Antet_M
do i =1,n
read(3,*) m(i, i:n)
enddo

read(3,*) m_factor
M =M*m_factor

! M-simetric:
do i=2,n
    do j=1,i-1
    m(i,j) =m(j,i)
    enddo
enddo

! K:
read(3,'(a80)') dummy
do i =1,n
read(3,*) k(i, i:n)
enddo

read(3,*) k_factor
K =K*k_factor

! K-simetric:
do i=2,n
    do j=1,i-1
    k(i,j) =k(j,i)
    enddo
enddo

close(3)

!open(1, file =Input_File, status ='old', form ='formatted', access ='Append')

!* Print:
Print_MK_If: &
IF(print_MK) THEN
write(7, '(/,a)') 'M'
do i=1,n
write(7, '(5(1pg23.16))') m(i,i:)
enddo

write(7, '(/,a)') 'K'
do i=1,n
write(7, '(5(1pg23.16))') k(i,i:)
enddo

ENDIF print_MK_If
!*

!write(7, '(//,a)') 'Lambda'
!do iv =1,n
!iwi =width(iv)
!iv_Sort =index_omega_Sort(iv)
!write(7, '(1x, I<iwi>, 1pg24.16)') iv, lambda(iv_Sort)    !omega(iv)
!enddo
!*

write(7,'(//,a)') ' Proba:'

write(7,'(a,13x,a,\)') ' K*x -lambda*M*x:'

! Proba: K*x -lambda*M*x ~0
dif =0

DO ivp =1,n
ivp_Sort =index_omega_Sort(ivp)
x1 =VEC_x(:,ivp_Sort)
lambda1 =lambda(ivp_Sort)

!x1 =VEC_x(:,ivp)
!lambda1 =lambda(ivp)

IF(full_proba) write(7, '(/, I3)') ivp  !ivp_Sort

y1 =matmul(K, x1)
y2 =lambda1*matmul(M, x1)

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

!ivmax =ivp
ivmax =ivp_Sort
icoord =i
endif
enddo

ENDDO

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

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

close(7)
ENDIF Proba_IF

END Program Main_Retrieve
