! Exact solution:
subroutine exact(t, icall)
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)

character(260) fname_Val, fname_Vec
character(80) dummy
character(260) lpstrTitle
integer status


! Defineste Constante:
Const_If: IF(icall==0) THEN

! Read eigenvalues:
lpstrTitle ='Select Eigen_Values File'C
call GetFile(fname_Val, lpstrTitle, status)

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

! Omega:
do id =1,2
read(2,*) dummy
enddo

do ii =1,ngl
read(2,*) i, om(i)
enddo

close(2)

! Read eigenvectors:
lpstrTitle ='Select Eigen_Vectors File'C
call GetFile(fname_Vec, lpstrTitle, status)

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

! Eigenvectors:
do id =1,2
read(2,*) dummy
enddo

do jj =1,ngl
read(2,*) j

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

enddo

close(2)


!pp(1) =25_8
pp0 =0D0;
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
