! Exact solution:
subroutine exact(t)
use work

implicit real(8) (a-h, o-z)
real(8) y(ngl), dy(ngl)
real(8) m_fi(ngl), fi_y(ngl), fi_dy(ngl)		! _ = *; fi(ngl, ngl)
real(8) mm, d, y_st
real(8) pp0(ngl)



! Eigenvalues
! 1   917.0400956800341    
! 2   207.5974745328022    
! 3   1815.362429787164    

! Eigenvectors (of R)
!          1                       2                       3
!   1.000000000000000       1.000000000000000       1.000000000000000    
! -0.7472706718779475      0.9249016223752760      -2.864636719006213    
! -0.7668418199290550      0.4027533908560142       4.095576068071734    

! Defineste Constante:
Const_If: IF(t ==tEnd) THEN

! Omega:
open (2, file ='Ex1-OMEGA.dat', status ='old', form ='formatted')
do it =1,3
read(2,'(a)') titlu
enddo

do igl =1, ngl

read(2,*) i, om(i)
enddo

close(2)

! Vect. proprii:
open (2, file ='Ex1-VEC_X.dat', status ='old', form ='formatted')

do it =1,3
read(2,'(a)') titlu
enddo

do igl =1, ngl
read(2,*) j		! vector #j

	do i =1, ngl
	read(2,*) fi(i,j)
	enddo
enddo

close(2)


!pp(1) =25_8
pp0 =matmul(transpose(fi), pp)

aa =0.; m_fi =0.;

do igl =1,ngl
d =1./( 1-(omp/om(igl))**2 )

m_fi =matmul(M, fi(:, igl))
mm =dot_product(fi(:,igl), m_fi)

y_st =pp0(igl)/( mm*om(igl)**2)

aa(igl) =y_st*d
enddo


!*
write(11,*) 'Omega:', om

write(11, '(/,a)') ' FI:'
do i =1,ngl
write(11,*) fi(i,:)
enddo

write(11,'(/,a)') 'PP0:'
write(11,*) pp0

write(11,'(/,a)') ' A:'
write(11,*) aa
!*

ENDIF Const_If


! Solutia:
do igl =1,ngl
y(igl)  =aa(igl)*dsin(omp*t) -omp*aa(igl)/om(igl)*dsin(om(igl)*t)
dy(igl) =omp*aa(igl)*( dcos(omp*t) -dcos(om(igl)*t) )
enddo

uExact = matmul(fi, y)

uPrimeExact = matmul(fi, dy)


end subroutine exact
