! 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)

real(8) S_1(ngl,ngl)

! Eigenvalues
! 1  0.9454183848254148    
! 2   3.478563691072702    
! 3  7.6017924101882628E-02

! Eigenvectors
!          1                       2                       3
!  0.7548686059843572     -0.5224038898372542      0.3965697461760532    
!  0.4735760529129831E-01  0.6464764854471262      0.7614626786553659    
! -0.6541640810380120     -0.5560236772358149      0.5127446005886904    

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

om(1) = 7.6017924101882628D-02
om(2) = 0.9454183848254148_8    
om(3) = 3.478563691072702_8    

om =dsqrt(om)


fi(1,:) =(/ 0.3965697461760532_8,    0.7548686059843572_8,     -0.5224038898372542_8 /)
fi(2,:) =(/ 0.7614626786553659_8,	 0.4735760529129831D-01,    0.6464764854471262_8 /)
fi(3,:) =(/ 0.5127446005886904_8,	-0.6541640810380120_8,     -0.5560236772358149_8 /)

! S_1:
!  0.707106781186547      -0.316227766016838       0.158113883008419     
!						   0.632455532033676      -0.316227766016838     
!												   0.790569415042095     


s_1(1,:) =(/  0.707106781186547_8,     -0.316227766016838_8,     0.158113883008419_8 /)
s_1(2,:) =(/	0._8,					0.632455532033676_8,    -0.316227766016838_8 /)    
s_1(3,:)   =(/	0._8,					0._8,				     0.790569415042095_8 /)


! Eigenvectors: FI =S_1*X
!fi =matmul(S_1, fi)
do i=1,ngl
fi(:,i) =matmul(S_1, fi(:,i))
enddo


!p(1) =25_8
!pp0 =matmul(transpose(fi), pp)
do igl =1, ngl
pp0(igl)=dot_product(fi(:,igl), pp)
enddo

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(10,*) 'Omega:', om

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

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

write(10,'(/,a)') ' A:'
write(10,*) 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
