PROGRAM NewOp_ngl
! "ngl" ecuatii (ord. 2)
! Dubla Precizie
! Cu sau fara mediere pt. U3 - pondere "niu"
! Coeficientii limita: p,teta, p1,teta1, p2,teta2
! Eroarea + Eroarea medie patratica - in Deplasare si Viteza.
!
! Ecuatii:
! - Neliniara, de forma:
! MU'' + G(U,U') +F(U) =P(t)
! Functiile G, F, P; jacobienii JF, JG; derivata DP: in FGP_J
!
! - Liniara:
! MU'' +CU' +KU =P(t)
! G = C; F =K;
!
! Febr. 2006
!
! Fisierul de intrare - continut:
!===============================
! 1) Titlul problemei
! 2) ngl
! 3) kod_M, kod_K	! cod pt. matricile M si K
!                   <0: nu se utilizeaza M (M=I), sau K; 
!                   =0: generata (subrutina);
!                   =1: citita din fisier exterior - matrice nxn
!                   >=2: citita triunghiul sup.

! 4) n_param: numar de parametri de care depinde sist. de ecuatii
!    4-1) valori parametri: numai pt. n_param >0
! 5) t0, h, tEnd	! 
! 6) u0, u10	! 
! 7) eps, lnit
! 8) kodm			! cod mediere: NU: =0; DA: #0
!	8-1) kodm, niu	! Daca kodm #0; niu : factor de mediere.
! 9) k_rez			! cod rezultate: 0/#0: pentru scriere fisier deplasari & viteze.
!	9-1) tStart_H, tEnd_H, ipas_H ! pentru k_rez #0
!   9-2) igl_Start, igl_End   ! gradele de libertate pt. care se scrie istoria
! 10) true_solution	! logic
! ========================================

use dfwin
use work

use DFLIB

implicit real(8) (a-h,o-z)
character(80) title
real(8), parameter :: PI=3.1415926535897932_8	!384626433832795

character(max_path) fname1
character(max_path), allocatable :: fname(:)

character(10) igl_char

real(4) s

!real(8) niu

common /coef/	beta_h3, gama_h2, delta_h, &
				teta_beta_h3, teta1_gama_h2, teta2_delta_h

data teta , beta  /.25_8, 6._8/       ! beta = 1/beta
data teta1, gama  / 3._8,  .5_8/        ! teta1= 1/teta1
data teta2, delta /.5_8,  1._8/

logical PrintIt
logical true_solution

real(8) norma_2

real(8), allocatable :: b(:)

!real(8) C10,C20,C1,C2,C11,C21,da,ust
!common/coef/C10,C20,C1,C2,C11,C21,da,ust

integer width
integer status
character(max_path) lpstrTitle

character(max_path) fname_M, fname_K
character(80) dummy

real(8), allocatable :: uError(:), uPrimeError(:)
real(8), allocatable :: Rel_uError(:), Rel_uPrimeError(:)


! coeficientii limita:
beta = 1d0/beta
teta1= 1d0/teta1
! delta = 1d0

lpstrTitle ='Select Input File'C
call GetFile(fname1, lpstrTitle, status)

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

! Read data; title.
read(1,'(a)') title

read(1,*) ngl


! Alocari:
allocate( u0(ngl), u10(ngl), u20(ngl), u30(ngl) )
u0=0; u10 =0; u20 =0; u30 =0;

allocate( uEnd(ngl), u1End(ngl), u2End(ngl), u3End(ngl) )
uEnd =0; u1End =0; u2End =0; u3End =0;

!allocate( m(ngl,ngl) )
allocate( f(ngl), g(ngl) )
allocate( jf(ngl, ngl), jg(ngl, ngl), jg1(ngl, ngl) )

allocate( p(ngl), dp(ngl) )
p =0; dp =0;

allocate( b(ngl) )
b =0;
allocate( pp(ngl) )
pp =0;

!read(1,*) e
read(1,*) kod_M, kod_K

    allocate( m(ngl,ngl) )
    m =0.
    do igl =1, ngl
    m(igl, igl) =1.
    enddo

if(kod_K >=0) then
    allocate( kk(ngl, ngl) )
    kk =0.
endif

read(1,*) n_param
if(n_param >=1) then

    allocate(param(n_param))
    read(1,*) param(:)
endif


read(1,*) t0, h, tEnd	!
read(1,*) u0, u10

read(1,*) eps, lnit		!, PrintIt

read(1,*)  kodm
	if(kodm /=0) read(1,*) niu

read(1,*) k_rez	! 0, #0 - for writing time-history: X, Xprime
	if(k_rez /=0) then
	read(1,*) tStart_H, tEnd_H, iPas_H
	read(1,*) igl_Start, igl_End
	endif

read(1,*) true_solution

!*
if(ngl ==1) then
C10 =0.; C20 =0.; C1 =0.; C2 =0.; C11 =0.; C21 =0.; da =0.; ust =0.;
call coef_lin()
endif
!*

! Coeficienti:
h2 = h*h; h3 = h2*h
beta_teta =beta*teta;	gama_teta1 =gama*teta1;	delta_teta2 =delta*teta2
beta_h3 =beta*h3; gama_h2= gama*h2; delta_h =delta*h
teta_beta_h3 = teta*beta_h3; teta1_gama_h2 =teta1*gama_h2; teta2_delta_h =teta2*delta_h


! Inital conditions:
! u0(ngl), u1(ngl)
!
! Calculeaza U20, U30:
! MU'' +G(U,U') +F(U) =P(t):  
!
! MU''' +JG(U,U')U' +JG1(U,U')U'' +JF(U)U' =DP(t)

!call Gen_M()
IF(kod_M ==0) THEN
call Gen_M()

ELSEIF(kod_M >=1) THEN    ! kod_M >=1
lpstrTitle ='Select Mass matrix "M" - Data File'C
call GetFile(fname_M, lpstrTitle, status)

open(unit =2, file =fname_M, status ='old', form ='formatted')
read(2,*) dummy

    if(kod_M ==1) then
    do i=1,ngl
    read(2,*) m(i,:)
    enddo
    
    else
    do i=1,ngl
    read(2,*) m(i,i:)
    enddo
   ! Umple simetric:
    do i =2,ngl
    do j =1,i-1
    m(i,j) =m(j,i)
    enddo
    enddo
    endif

close(2)
ENDIF

IF(kod_K ==0) THEN
call Gen_K()

ELSEIF(kod_K >=1) THEN     ! kod_K >=1
lpstrTitle ='Select Stiffness matrix "K" - Data File'C
call GetFile(fname_K, lpstrTitle, status)

open(unit =3, file =fname_K, status ='old', form ='formatted')
read(3,*) dummy

    if(kod_K ==1) then
    do i=1,ngl
    read(3,*) kk(i,:)
    enddo
        
    else
    do i=1,ngl
    read(3,*) kk(i,i:)
    enddo
   ! Umple simetric:
    do i =2,ngl
    do j =1,i-1
    kk(i,j) =kk(j,i)
    enddo
    enddo
    endif

close(3)
ENDIF

                                             
call FGP_J(t0, u0, u10)	! Cu t0, u0, u10, etc.

! M*U20 = P - G - F
b = P - G - F
kod =0; kodPrint =0;
call ELIM (ngl, M, b, prag, kod, kodPrint)	!
U20 =b

!*
!write(1,*); write(1,*) 'U20:', u20;
if(maxval(dabs(U20)) ==0) then
write(1,*); write(1,*) 'U20 =0.'

else
write(1,*); write(1,*) 'U20:', u20;
endif


! M*U30 = dP0 - jG*U10 - jG1*U20 - jF*U10
DO l=1,NGL
b1 =0d0
do j=1,ngl
b1 = b1 + jG(l,j)*U10(j) + jG1(l,j)*U20(j) + jF(l,j)*U10(j)
enddo

b(l) =dP(l) -b1
ENDDO

kod =0; kodPrint =0;
call Gen_M()

call ELIM (ngl, M, b, prag, kod, kodPrint)	!
U30 =b

!*
!write(1,*) 'U30:', u30
if(maxval(dabs(U30)) ==0) then
write(1,*); write(1,*) 'U30 =0.'

else
write(1,*); write(1,*) 'U30:', u30;
endif


No_Intervals =(tEnd-t0)/h +0.5

allocate (u(0:No_Intervals, ngl), uPrime(0:No_Intervals, ngl))

! Initialize:
u(0,:) =u0; uPrime(0,:) =u10
u(1:No_Intervals,:) =0.0; uPrime(1:No_Intervals,:) =0.0


! Open for time-history:
IF(k_rez /=0) THEN
n_igl =igl_End -igl_Start +1
allocate(fname(1:n_igl))

iPoint =index(fname1,'.',back=.TRUE.)

!iunit =6
do i =1, n_igl		! ngl <=9
igl =igl_Start +(i-1)
ic =inttochar(igl, igl_char, 10)
l_i =lstrlen(igl_char)

fname(i) =fname1(1:iPoint-1)//'-'//igl_char(:l_i)//'.dat'

iunit =6+i
open(iunit, file=fname(i), status='unknown')

enddo


ENDIF


print *, Title
print *, 'ngl:', ngl
print *, 't0, Step:', t0,h
print *, 'Initial cond.(system):',u0, u10
print *, 'No.of Intervals:', No_Intervals

print '(/, a)', ' * Working ...'


PrintIt =.false.

! Start Time:
s = secnds(0.0)
! Time-history integration:
!call gen_M()	! pentru Newop -> sub_ab

NewOp_DO: DO i =1, No_Intervals
t =t0+ i*h

kod_elim =0
call NewOp(t, PrintIt, kod_elim)

if(kod_elim <0) exit NewOp_DO

u(i,:) =uEnd; uPrime(i,:) =u1End
u0 =uEnd; u10 =u1End; u20 =u2End; u30 =u3End

ENDDO NewOp_DO

s = secnds(s)

if(kod_elim <0) then

write(1,'(/, a)') ' Pivot < prag. Solutia nu se poate obtine.'
goto 500
endif

iw =width(No_Intervals)
write(1, '(/, a, I<iw>)') ' No.of Intervals: ', No_Intervals

if(kodm.ne.0) write(1, '(/, a, F6.2)') '  *** Mediere: niu =', niu

print '(/, a, f10.2, a)',' Integration Time =', s,' sec.'

write(1,'(/, a,f10.2,a)') ' Integration Time =',s,' sec.'

write(1, '(/, a, g12.5)') ' tEnd:', t;

n_u =min(ngl, 6)
!write(1,'(/,a, <ngl>(1pG21.11))') ' End: U : ', u(No_intervals, :)
write(1,'(/,a)') ' End: U : '
write(1,'(<n_u>(1pG21.11))') u(No_intervals, :)

!write(1,'(a, <ngl>(1pG21.11))')   '      U'': ', uPrime(No_intervals, :)
write(1,'(a)')   '      U'': '
write(1,'(<n_u>(1pG21.11))')  uPrime(No_intervals, :)



IF(k_rez /=0) THEN
iStart =(tStart_H -t0)/h
iEnd =(tEnd_H -t0)/h

print '(/,a)',' Writing history ...'

IGL_DO: DO i =1,n_igl

igl =igl_Start +(i-1)
iunit =6+i

His_Do: DO it =iStart, iEnd, iPas_H
write(iunit, *) it*h, u(it, igl)

enddo His_Do

close(iunit)

ENDDO IGL_DO

!close(2); close(3); close(4); ;close(5)

ENDIF


deallocate( u0, u10, u20, u30 )
!deallocate( m )
deallocate( f, g )
deallocate( jf, jg, jg1 )
deallocate( p, dp )
deallocate( b )

! Now, calculate the error at tEnd and Maximumm error:

True_If: IF(true_solution) THEN
allocate( uExact(ngl), uError(ngl), uPrimeExact(ngl), uPrimeError(ngl) )
allocate(Rel_uError(ngl), Rel_uPrimeError(ngl))
uError =0.; uPrimeError =0.; Rel_uError =0; Rel_uPrimeError =0.;

allocate(om(ngl), aa(ngl), fi(ngl,ngl))

! Error:
call exact(tEnd, 0)		! uExact, uPrimeExact

!*
if(ngl <=6) then
write(1,'(/, a, <ngl>(1pG21.11))') ' U_exact :', uExact(:)
write(1,'(a, <ngl>(1pG21.11))') ' U''_exact:', uPrimeExact(:)
endif
!*

uError =uExact -uEnd
uPrimeError =uPrimeExact -u1End

do igl =1, ngl
if(uExact(igl) /=0) then
Rel_uError(igl) =uError(igl)/uExact(igl)
endif

if(uPrimeExact(igl) /=0) then
Rel_uPrimeError(igl) =uPrimeError(igl)/uPrimeExact(igl)
endif
enddo

write(1, '(/, a)') ' Relative Error at tEnd:'
!write(1,'(a, <ngl>(2x,1pD10.3))') ' U :', uError
write(1,'(a)') ' U :'
write(1,'(<n_u>(2x,1pD10.3))') Rel_uError

!write(1,'(a, <ngl>(2x,1pD10.3))') ' U'':', uPrimeError
write(1,'(a)') ' U'':'
write(1,'(<n_u>(2x,1pD10.3))') Rel_uPrimeError


! Maximum Error:
uError_max =0D0
uPrimeError_max =0D0

I_DO: DO i =1, No_Intervals
t =t0 +i*h
call exact(t, i)		! uExact, uPrimeExact

! Error:
uError =uExact -u(i,:)
uPrimeError =uPrimeExact -uPrime(i,:)

uNorm =norma_2(uError)

if(uNorm > uError_max) then
uError_max =uNorm
t_max =t
endif

uPrimeNorm =norma_2(uPrimeError)

if(uPrimeNorm > uPrimeError_max) then
uPrimeError_max =uPrimeNorm
t1_max =t
endif

END DO I_DO

write(1, '(/, a)') ' Maximum Error Norm in U and U'''
write(1,'(a, f10.6, a, 2x,1pD10.3)') ' U:', t_max,  ':', uError_max
write(1,'(a, f10.6, a, 2x,1pD10.3)') ' U''',t1_max, ':', uPrimeError_max

ENDIF True_If

500 continue
close (1)
deallocate( m )
deallocate(u, uPrime)
deallocate( uEnd, u1End, u2End, u3End )
deallocate( pp )

if(true_solution) then
deallocate( uExact, uError, uPrimeExact, uPrimeError, Rel_uError, Rel_uPrimeError )

deallocate(om, aa)
endif

print '(/, a)', ' * Finished.'

end program NewOp_ngl
