!>
!> @file driv4.f90
!>
!> @brief
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> SPClibs is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> SPClibs is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @author
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
! Computation of croos mass matrix between two splines sp1 & sp2 
! sp1 and sp2 can be splines of any type (i.e. either set up with set_spline or
! set_splcoef) and of any order.

!
  USE bsplines
  USE matrix
!
  IMPLICIT NONE
  TYPE(gbmat) :: matm
  INTEGER :: mrows, ncols, kl, ku
  DOUBLE PRECISION, ALLOCATABLE :: avec(:,:), bvec(:,:), matfull(:,:)
!
  TYPE(zgbmat) :: zmatm
  DOUBLE COMPLEX, ALLOCATABLE :: zavec(:,:), zbvec(:,:)
  DOUBLE PRECISION :: dznrm2
!
  TYPE(spline1d) :: sp1, sp2
  INTEGER :: nx, nidbas1, nidbas2, ngauss
  INTEGER :: i, j
  DOUBLE PRECISION :: a, b, coefx(5)
  DOUBLE PRECISION, ALLOCATABLE :: xgrid(:)
  DOUBLE PRECISION, DIMENSION(:, :), POINTER :: MassMat
  LOGICAL :: periodic1, periodic2
  NAMELIST /newrun/ nx, a, b, coefx, nidbas1, nidbas2, periodic1, periodic2
!===========================================================================
!                  1.0 Set up grids
!
!   Read in data specific to run
!
  nx = 8    ! Number of intevals in x
  a = 0.0d0  ! Left boundary of interval
  b = 1.0d0  ! Right boundary of interval
  coefx(1:5) = (/1.0d0, 0.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
  periodic1 = .FALSE.
  periodic2 = .FALSE.
  nidbas1 = 3
  nidbas2 = 2
  READ(*,newrun)
  WRITE(*,newrun)
!
!   Define grid/knots
!
  ALLOCATE(xgrid(0:nx))
  xgrid(0 ) = a
  xgrid(nx) = b
  CALL meshdist(coefx, xgrid, nx)
  WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
!===========================================================================
!                  2.0 Set up splines
!
  ngauss = 1 ! Gauss points initialized with set_spline are in fact not used
             ! for computing cross mass matrix
  ! First spline set up as for solving a PDE with FEMs
  CALL set_spline(nidbas1, ngauss, xgrid, sp1, periodic1)

  ! Second spline set up as for interpolation
!!$  CALL set_splcoef(nidbas2, xgrid, sp2, periodic2)
  CALL set_spline(nidbas2, ngauss, xgrid, sp2, periodic2)

  WRITE(*,'(/a,i6,a,i6/(10(f8.3)))') 'KNOTS of spline sp1', LBOUND(sp1%knots), &
       &                       ':',UBOUND(sp1%knots), sp1%knots
  WRITE(*,'(/a,i6,a,i6/(10(f8.3)))') 'KNOTS of spline sp2', LBOUND(sp2%knots), &
       &                       ':',UBOUND(sp2%knots), sp2%knots
  WRITE(*,'(3(a,i5, 2x))') 'NX =', nx, 'DIM sp1 =', sp1%dim, 'DIM sp2 =', sp2%dim
!===========================================================================
!                  3.0 Compute cross mass matrix
!
  CALL CompMassMatrix(sp1, sp2, a, b, MassMat)

  WRITE(*, "(a)") "Cross-mass matrix between splines sp1 & sp2:"
  DO i = 1, SIZE(MassMat, 1)
     WRITE(*, "(15f10.5)") (MassMat(i, j), j = 1, MIN(SIZE(MassMat, 2), 15))
  END DO
!
!  Should equal to 1 for splines i "not close to  the boundaries":
!                  p1 .LT. i .LE. N
!
  WRITE(*,'(/a/(15f8.5))') 'Sum of cols * NX', SUM(MassMat,dim=2)*REAL(nx,8)
!===========================================================================
!                  3.0  Use DGB matrice
!
!!$  mrows = nx+nidbas1
!!$  ncols = nx+nidbas2
  CALL get_dim(sp1, mrows)
  CALL get_dim(sp2, ncols)
  kl = nidbas1
  ku = nidbas2
  CALL init(kl, ku, ncols, 1, matm, mrows=mrows)
  WRITE(*,'(/a, 2i3)') 'Band matrix:, kl, ku =', kl, ku
!
  CALL CompMassMatrix(sp1, sp2, a, b, matm)
!
  DO i=1,SIZE(matm%val,1)
     WRITE(*,'(15f10.5)') matm%val(i,:)
  END DO
!
  WRITE(*,'(/a)') 'Full matrix'  
  ALLOCATE(matfull(mrows,ncols))
  matfull = 0.0d0
  DO i=1,mrows
     CALL getrow(matm, i, matfull(i,:))
     WRITE(*,'(15f10.5)') matfull(i,:)
  END DO
  WRITE(*,'(a,1pe12.4)') 'error =', MAXVAL(ABS(matfull-MassMat))
!
!   Check VMX
  ALLOCATE(avec(ncols,2))
  ALLOCATE(bvec(mrows,2))
  avec = 1.0d0
  bvec = vmx(matm,avec)*REAL(nx,8)
  WRITE(*,'(a)') 'M*a, with a=1'
  DO j=1,2
     WRITE(*,'(15f8.5)') bvec(:,j)
  END DO
!===========================================================================
!                  4.0 Test complex version
!
  CALL init(kl, ku, ncols, 1, zmatm, mrows=mrows)
  CALL CompMassMatrix(sp1, sp2, a, b, zmatm)
  ALLOCATE(zavec(ncols,2))
  ALLOCATE(zbvec(mrows,2))
  zavec = (1.0d0,0.0d0)
  zbvec = vmx(zmatm,zavec)*REAL(nx,8)
  zbvec = zbvec-bvec
  WRITE(*,'(/a)') 'Check complex version'
  WRITE(*,'(a,2(1pe12.4))') 'Norm of errors =', &
       &              (dznrm2(mrows, zbvec(1,j), 1), j=1,2)
  zmatm%val = zmatm%val-matm%val
  WRITE(*,'(a,1pe12.4)') 'Diff of matrix elements =', MAXVAL(ABS(zmatm%val))
!===========================================================================
!                  9.0 Epilogue
!
  DEALLOCATE(MassMat)
  DEALLOCATE(xgrid)
  CALL destroy_sp(sp1)
  CALL destroy_sp(sp2)
  call destroy(matm)

END PROGRAM main
!+++
SUBROUTINE meshdist(c, x, nx)
!
!   Construct an 1d  non-equidistant mesh given a
!   mesh distribution function.
!
  IMPLICIT NONE
  DOUBLE PRECISION, INTENT(in) :: c(5)
  INTEGER, INTENT(iN) :: nx
  DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
  INTEGER :: nintg
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xint, fint
  DOUBLE PRECISION :: a, b, dx, f0, f1, scal
  INTEGER :: i, k
!
  a=x(0)
  b=x(nx)
  nintg = 10*nx
  ALLOCATE(xint(0:nintg), fint(0:nintg))
!
!  Mesh distribution
!
  dx = (b-a)/REAL(nintg)
  xint(0) = a
  fint(0) = 0.0d0
  f1 = fdist(xint(0))
  DO i=1,nintg
     f0 = f1
     xint(i) = xint(i-1) + dx
     f1 = fdist(xint(i))
     fint(i) = fint(i-1) + 0.5*(f0+f1)
  END DO
!
!  Normalization
!
  scal = REAL(nx) / fint(nintg)
  fint(0:nintg) = fint(0:nintg) * scal
!!$  WRITE(*,'(a/(10f8.3))') 'FINT', fint
!
!  Obtain mesh point by (inverse) interpolation
!
  k = 1
  DO i=1,nintg-1
     IF( fint(i) .GE. REAL(k) ) THEN
        x(k) = xint(i) + (xint(i+1)-xint(i))/(fint(i+1)-fint(i)) * &
             &   (k-fint(i))
        k = k+1
     END IF
  END DO
!
  DEALLOCATE(xint, fint)
CONTAINS
  DOUBLE PRECISION FUNCTION fdist(x)
    DOUBLE PRECISION, INTENT(in) :: x
    fdist = c(1) + c(2)*x + c(3)*EXP(-((x-c(4))/c(5))**2)
  END FUNCTION fdist
END SUBROUTINE meshdist
!+++