!>
!> @file pde1d_eig_ge.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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!  Solving the following 1d differential  eqation using splines:
!
!    -d/dr[r d/dr] f = k^2 r^(k-1), with f(r=1) = 0
!    exact solution: f(r) = 1 - r^k
!
  USE bsplines
  USE matrix
  USE futils
  USE conmat_mod
  IMPLICIT NONE
  INTEGER :: nx, nidbas, ngauss, kdiff
  INTEGER :: i, nrank, kl, ku
  LOGICAL :: nlppform
  DOUBLE PRECISION :: coefs(5)
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: xgrid, rhs, sol
  TYPE(spline1d) :: splx
  TYPE(gemat) :: mat
!
  CHARACTER(len=128) :: file='pde1d.h5'
  INTEGER :: fid
  DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: arr
  DOUBLE PRECISION :: seconds, t0, tmat, teig, tarpack
!
  DOUBLE PRECISION, ALLOCATABLE :: mata(:,:), eigvecsr(:,:), eigvecsl(:,:), &
       &                           wr(:), wi(:), work(:)
  INTEGER :: j, info
!
  INTEGER, PARAMETER :: maxn=256, maxnev=maxn, maxncv=maxn
  DOUBLE PRECISION :: v(maxn,maxncv),workd(3*maxn), workev(3*maxncv), &
       &              d(maxncv,2), resid(maxn), w(maxncv,maxn), &
       &              zero=0.0d0, tol=0.0d0, sigmar, sigmai
  DOUBLE PRECISION, ALLOCATABLE :: workl(:)
  DOUBLE PRECISION, EXTERNAL :: dnrm2
  INTEGER :: nev=10, ncv=30, ishfts=1, maxitr=300, nconv, &
       &     mode1=1, ierr, lworkl
  INTEGER :: ido, ipntr(11), iparam(11)
  
  CHARACTER(len=1) :: bmat='I'
  CHARACTER(len=2) :: which='SA'
  LOGICAL :: rvec, select(maxncv)
  
!
  INTERFACE
     SUBROUTINE dismat(spl, mat)
       USE bsplines
       USE matrix
       TYPE(spline1d), INTENT(in) :: spl
       TYPE(gbmat), INTENT(inout) :: mat
     END SUBROUTINE dismat
     SUBROUTINE disrhs(kdiff, spl, rhs)
       USE bsplines
       INTEGER, INTENT(in) :: kdiff
       TYPE(spline1d), INTENT(in) :: spl
       DOUBLE PRECISION, INTENT(out) :: rhs(:)
     END SUBROUTINE disrhs
     SUBROUTINE meshdist(coefs, x, nx)
       DOUBLE PRECISION, INTENT(in) :: coefs(5)
       INTEGER, INTENT(iN) :: nx
       DOUBLE PRECISION, INTENT(inout) :: x(0:nx)
     END SUBROUTINE meshdist
  END INTERFACE
!
  NAMELIST /newrun/ nx, nidbas, ngauss, kdiff, nlppform, coefs, &
       &            nev, ncv, which
!===========================================================================
!              1.0 Prologue
!
!   Read in data specific to run
!
  nx = 8     ! Number oh intevals in x
  nidbas = 3 ! Degree of splines
  ngauss = 4 ! Number of Gauss points/interval
  kdiff = 2  ! Exponent of differential problem
  nlppform = .TRUE. ! Use PPFORM for gridval or not
  coefs(1:5) = (/0.2d0, 1.d0, 0.d0, 0.d0, 1.d0/) ! Mesh point distribution function
                                                 ! see function FDIST in MESHDIST
!
  READ(*,newrun)
  WRITE(*,newrun)
!
!   Define grid on x axis
!
  ALLOCATE(xgrid(0:nx))
  xgrid(0) = 0.0d0
  xgrid(nx) = 1.0d0
  CALL meshdist(coefs, xgrid, nx)
  WRITE(*,'(a/(10f8.3))') 'XGRID', xgrid(0:nx)
  WRITE(*,'(a/(10f8.3))') 'Diff XGRID', (xgrid(i)-xgrid(i-1), i=1,nx)
  WRITE(*,'(a/(10f8.3))') 'Diff XGRID**2', (xgrid(i)**2-xgrid(i-1)**2, i=1,nx)
!
!   Create hdf5 file
!
  CALL creatf(file, fid, 'PDE1D Result File')
  CALL attach(fid, '/', 'NX', nx)
  CALL attach(fid, '/', 'NIDBAS', nidbas)
  CALL attach(fid, '/', 'NGAUSS', ngauss)
  CALL attach(fid, '/', 'KDIFF', kdiff)
!===========================================================================
!              2.0 Discretize the PDE
!
!   Set up spline
!
  t0 = seconds()
  CALL set_spline(nidbas, ngauss, xgrid, splx, nlppform=nlppform)
  CALL get_dim(splx, nrank)       ! Rank of the FE matrix
  WRITE(*,'(a,l6)') 'splx%nlequid', splx%nlequid
!
!   FE matrix assembly
!
  kl = nidbas
  ku = kl
!!$  CALL init(kl, ku, nrank, 1, mat)
  WRITE(*,'(a,3i6)') 'kl, ku, nrank', kl, ku, nrank
  CALL init(nrank, 1, mat)
!!$  CALL dismat(splx, mat)
  CALL conmat(splx, mat, coefeq)
!
  ALLOCATE(arr(nrank))
!!$  WRITE(*,'(/a)') 'Matrice before BC'
!!$  DO i=1,nrank
!!$     CALL getrow(mat, i, arr)
!!$     WRITE(*,'(12f8.3)') arr, SUM(arr)
!!$  END DO
!
!   RHS assembly
!
  ALLOCATE(rhs(nrank), sol(nrank))
  CALL disrhs(kdiff, splx, rhs)
!!$  WRITE(*,'(/a,/(12(f8.3)))') 'RHS before BC', rhs
!
!   Set BC f(r=1) = 0 on matrix
!
  arr(1:nrank-1) = 0.0d0
  arr(nrank) = 1.0d0
  CALL putrow(mat, nrank, arr)
  CALL putcol(mat, nrank, arr)
  tmat = seconds() - t0
!
!===========================================================================
!              3.0  Eigevalue problem
!
!   Using Lapack dgeev
!
  t0 = seconds()
  ALLOCATE(mata(nrank,nrank))
  ALLOCATE(eigvecsr(nrank,nrank))
  ALLOCATE(eigvecsl(nrank,nrank))
  ALLOCATE(work(4*nrank))
  ALLOCATE(wr(nrank), wi(nrank))
  mata(:,:) = mat%val(:,:)
  CALL putarr(fid, '/MAT', mata, 'matrix A')
!!$  CALL dsyev('V', 'L', nrank, mata, nrank, eigvals, work, SIZE(work), info)
  CALL dgeev('N', 'V', nrank, mata, nrank, wr, wi, eigvecsl, SIZE(eigvecsl,1), &
       &     eigvecsr, SIZE(eigvecsr,1), work, SIZE(work), info)
  teig = seconds() - t0
  PRINT*,'Info from DGEEV', info, arr(1)
  IF(info.EQ.0) THEN
     CALL putarr(fid, '/REIGVS', wr, 'Real of eigenvalues of A')
     CALL putarr(fid, '/IEIGVS', wi, 'Imag of eigenvalues of A')
     CALL putarr(fid, '/EIGVECL', eigvecsl, 'left eigenvalues of A')
     CALL putarr(fid, '/EIGVECR', eigvecsr, 'right eigenvalues of A')
     WRITE(*,'(a/(10f10.4))') 'Real eigval', wr
     WRITE(*,'(a/(10f10.4))') 'Imag eigval', wi
  END IF
!
!   Using Arpack
!
  ido = 0
  iparam(1) = ishfts
  iparam(3) = maxitr
  iparam(7) = mode1
!
 lworkl  = 3*ncv**2+6*ncv
 ALLOCATE(workl(lworkl))
  
!
  t0 = seconds()
  DO   ! ARPACK reverse communication loop 
     CALL dnaupd(ido, bmat, nrank, which, nev, tol, resid, &
          &      ncv, v, maxn, iparam, ipntr, workd, workl,&
          &      lworkl, info)
     IF(ido.EQ.-1 .OR. ido.EQ.1) THEN
        CALL av(nrank,workd(ipntr(1)), workd(ipntr(2)))
        CYCLE
     END IF
     PRINT*, 'INFO =', info
     IF(info .LT. 0) THEN
        PRINT*, 'Error in dnaupd with info =', info
     ELSE
        rvec = .TRUE.
        CALL dneupd(rvec, 'A', select, d, d(1,2), v, size(v,1), &
             &      sigmar, sigmai, workev, bmat, nrank, which, nev, tol, &
             &      resid, ncv, v, size(v,1), iparam, ipntr, workd, workl, &
             &      lworkl, ierr )
        IF( ierr .NE. 0 ) THEN
           PRINT*,'Error in dneupd with ierr =', ierr
        ELSE
           nconv = iparam(5)
           PRINT*, '--- Real eigenvalues and comprison with Lapack results ---'
! eiegvalues and diff with Lapack results
           DO j=1,nconv
              WRITE(*,'(i3,3(1pe12.4))') j, d(j,1), wr(j), wr(j)-d(j,1)
           END DO
           PRINT*, '--- Imag eigenvalues and comprison with Lapack results ---'
           DO j=1,nconv
              WRITE(*,'(i3,3(1pe12.4))') j, d(j,2), wi(j), wi(j)-d(j,2)
           END DO
!!$              CALL av(nrank,v(1,j), w)   ! d(1,j) is the j^th eigenvalue
!!$              CALL daxpy(nrank, -d(j,1), v(1,j), 1, w, 1)
!!$              d(j,2) = dnrm2(nrank, w, 1)
!!$              d(j,2) = d(j,2)/abs(d(j,1))
           CALL putarr(fid, '/EIGVS', d(1:nconv,:), 'ARPACK eigenvalues of A')
!!$              CALL putarr(fid, '/EIGVECS', v, 'ARPACK eigenvectors of A')
!!$           END DO
        END IF
        EXIT
     END IF
  END DO   ! End of ARPACK reverse communication loop
  tarpack = seconds()-t0
  
!===========================================================================
!              9.0  Epilogue
!
  WRITE(*,'(a)') '---'
  WRITE(*,'(a,1pe12.3)') 'Matrice construction time (s) ', tmat
  WRITE(*,'(a,1pe12.3)') 'Lapack: Matrice eigvalue time (s)', teig
  WRITE(*,'(a,1pe12.3)') 'Arpack: Matrice eigvalue time (s)', tarpack
!
  DEALLOCATE(xgrid, rhs, sol)
  DEALLOCATE(arr)
  CALL destroy(mat)
  CALL destroy_sp(splx)
  CALL closef(fid)
CONTAINS
  SUBROUTINE av(n,v,w)
!
! Matrix vector product: w <- Av
    INTEGER, INTENT(in) :: n
    DOUBLE PRECISION, INTENT(in) :: v(*)
    DOUBLE PRECISION, INTENT(out) :: w(*)
    w(1:n) = vmx(mat,v(1:n))
  END SUBROUTINE av
  SUBROUTINE coefeq(x, idt, idw, c)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(out) :: idt(:), idw(:)
    DOUBLE PRECISION, INTENT(out) :: c(:)
!
! Weak form = Int(x*dw/dx*dt/dx)dx
!
    c(1) = x
    idt(1) = 1
    idw(1) = 1
  END SUBROUTINE coefeq
END PROGRAM main
!
!+++
SUBROUTINE meshdist(c, x, nx)
!
!   Construct an 1d  non-equidistant mesh given a
!   mesh distribution function defined in FDIST
!
  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
!+++
SUBROUTINE dismat(spl, mat)
!
!   Assembly FE matrix mat using spline spl
!
  USE bsplines
  USE matrix
  IMPLICIT NONE
  TYPE(spline1d), INTENT(in) :: spl
  TYPE(gbmat), INTENT(inout) :: mat
  INTEGER :: nrank, nx, nidbas, ngauss
  INTEGER :: i, igauss, iterm, iw, jt, irow, jcol
  DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
  DOUBLE PRECISION :: contrib
!
  INTEGER :: kterms                ! Number of terms in weak form
  INTEGER, ALLOCATABLE :: idert(:), iderw(:)  ! Derivative order
  DOUBLE PRECISION, ALLOCATABLE :: coefs(:)
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  CALL get_dim(spl, nrank, nx, nidbas)
  ALLOCATE(fun(0:nidbas,0:1)) ! Spline and its 1st derivative
!
!   Weak form
!
  kterms = mat%nterms
  ALLOCATE(idert(kterms), iderw(kterms), coefs(kterms))
!
!   Gauss quadature
!
  CALL get_gauss(spl, ngauss)
  ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
  DO i=1,nx
     CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
     DO igauss=1,ngauss
        CALL basfun(xgauss(igauss), spl, fun, i)
        CALL coefeq(xgauss(igauss), idert, iderw, coefs)
        DO iterm=1,kterms
           DO jt=0,nidbas
              DO iw=0,nidbas
                 contrib = fun(jt,idert(iterm)) * coefs(iterm) * &
                      &    fun(iw,iderw(iterm)) * wgauss(igauss)
                 irow=i+iw; jcol=i+jt
                 CALL updtmat(mat, irow, jcol, contrib)
              END DO
           END DO
        END DO
     END DO
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(fun)
  DEALLOCATE(xgauss, wgauss)
  DEALLOCATE(iderw, idert, coefs)
!
CONTAINS
  SUBROUTINE coefeq(x, idt, idw, c)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(out) :: idt(:), idw(:)
    DOUBLE PRECISION, INTENT(out) :: c(SIZE(idt))
!
! Weak form = Int(x*dw/dx*dt/dx)dx
!
    c(1) = x
    idt(1) = 1
    idw(1) = 1
  END SUBROUTINE coefeq
END SUBROUTINE dismat
!+++
SUBROUTINE disrhs(kdiff, spl, rhs)
!
!   Assenbly the RHS using spline spl
!
  USE bsplines
  IMPLICIT NONE
  INTEGER, INTENT(in) :: kdiff
  TYPE(spline1d), INTENT(in) :: spl
  DOUBLE PRECISION, INTENT(out) :: rhs(:)
  INTEGER :: nrank, nx, nidbas, ngauss
  INTEGER :: i, igauss, left
  DOUBLE PRECISION, ALLOCATABLE :: fun(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: xgauss(:), wgauss(:)
  DOUBLE PRECISION:: contrib
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  CALL get_dim(spl, nrank, nx, nidbas)
!!$  WRITE(*,'(/a, 5i3)') 'nrank, nx, nidbas =', nrank, nx, nidbas
!
  ALLOCATE(fun(nidbas+1,1)) ! needs only basis functions (no derivatives)
!
!   Gauss quadature
!
  CALL get_gauss(spl, ngauss)
!!$  WRITE(*,'(/a, i3)') 'Gauss points and weights, ngauss =', ngauss
  ALLOCATE(xgauss(ngauss), wgauss(ngauss))
!===========================================================================
!              2.0 Assembly loop
!
  rhs(1:nrank) = 0.0d0
!
  DO i=1,nx
     CALL get_gauss(spl, ngauss, i, xgauss, wgauss)
!!$     WRITE(*,'(a,i3,(8f10.4))') 'i=',i, xgauss, wgauss
     DO igauss=1,ngauss
        CALL basfun(xgauss(igauss), spl, fun, i)
        left = i-1
!!$        WRITE(*,'(a,5i4)') '>>>igauss, left', igauss, left
        contrib = wgauss(igauss)*rhseq(xgauss(igauss), kdiff)
        rhs(i:i+nidbas) = rhs(i:i+nidbas) + contrib*fun(:,1)
     END DO
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(fun)
  DEALLOCATE(xgauss, wgauss)
!
CONTAINS
  DOUBLE PRECISION FUNCTION rhseq(x,k)
    DOUBLE PRECISION, INTENT(in) :: x
    INTEGER, INTENT(in) :: k
    rhseq = k*k*x**(k-1)
  END FUNCTION rhseq
END SUBROUTINE disrhs