!>
!> @file conmat2.tpl
!>
!> @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>
!>
!
!   In this version local matrices E and D are precalculted
!   (on all n1/n2 intervals and nterms weak-form terms
!
  INTERFACE
     SUBROUTINE coefeq(x, y, idt, idw, c)
       DOUBLE PRECISION, INTENT(in) :: x, y
       INTEGER, INTENT(out) :: idt(:,:), idw(:,:)
       DOUBLE PRECISION, INTENT(out) :: c(:)
     END SUBROUTINE coefeq
  END INTERFACE
  INTEGER, OPTIONAL :: maxder(2)   ! maximum oder of derivatives
!
  INTEGER :: n1, nidbas1, ndim1
  INTEGER :: n2, nidbas2, ndim2
  INTEGER :: ng1, ng2
  INTEGER :: i1, i2, ig1, ig2
  INTEGER :: igt1, igt2, igw1, igw2, irow, jcol
  INTEGER, ALLOCATABLE :: left1(:), left2(:)
!
  LOGICAL :: nlper1, nlper2
!
  INTEGER :: kterms         ! Number of terms in weak form
  INTEGER :: k, kmaxder, it1, iw1, it2, iw2
  INTEGER, ALLOCATABLE :: idert(:,:), iderw(:,:)  ! Derivative order
  DOUBLE PRECISION, ALLOCATABLE  :: coefs(:,:,:)  ! Terms in weak form
  DOUBLE PRECISION :: dummy(mat%nterms)
!
  DOUBLE PRECISION, POINTER :: xg1(:,:), xg2(:,:), wg1(:,:), wg2(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: fun1(:,:,:), fun2(:,:,:)
  DOUBLE PRECISION, ALLOCATABLE :: mata(:,:,:,:), matc(:,:), matcd(:,:,:)
  DOUBLE PRECISION, ALLOCATABLE :: matd(:,:,:,:,:), mate(:,:,:,:,:)
!===========================================================================
!              1.0 Prologue
!
!   Properties of spline space
!
  CALL get_dim(spl%sp1, ndim1, n1, nidbas1)
  CALL get_dim(spl%sp2, ndim2, n2, nidbas2)
  nlper1 = spl%sp1%period
  nlper2 = spl%sp2%period
!
!   Gauss points and weights on all intervals
!
  xg1 => spl%sp1%gausx  ! xg1(ng1,n1)
  wg1 => spl%sp1%gausw  ! wg1(ng1,n1)
  ng1 = SIZE(xg1,1)
  xg2 => spl%sp2%gausx
  wg2 => spl%sp2%gausw
  ng2 = SIZE(xg2,1)
!
!   Derivative orders in the weak form
!
  kterms = mat%nterms
  ALLOCATE(idert(kterms,2))
  ALLOCATE(iderw(kterms,2))
  CALL coefeq(xg1(1,1), xg2(1,1), idert, iderw, dummy)
!
!   Precalc matrix E  in dimension 1
!
  kmaxder = 1
  IF(PRESENT(maxder)) kmaxder = maxder(1)
  ALLOCATE(fun1(0:nidbas1,0:kmaxder,ng1))
  ALLOCATE(left1(ng1))
  ALLOCATE(mate(0:nidbas1,0:nidbas1,ng1,kterms,n1))
  DO i1=1,n1
     left1 = i1
     CALL basfun(xg1(:,i1), spl%sp1, fun1, left1)
     DO k=1,kterms
        DO ig1=1,ng1
           DO iw1=0,nidbas1
              DO it1=0,nidbas1
                 mate(it1,iw1,ig1,k,i1) = wg1(ig1,i1) * &
                      &  fun1(it1,idert(k,1),ig1) * &
                      &  fun1(iw1,iderw(k,1),ig1)
              END DO
           END DO
        END DO
     END DO
  END DO
  DEALLOCATE(fun1)
  DEALLOCATE(left1)
!
!   Precalc matrix D  in dimension 2
!
  kmaxder = 1
  IF(PRESENT(maxder)) kmaxder = maxder(2)
  ALLOCATE(fun2(0:nidbas2,0:kmaxder,ng2))
  ALLOCATE(left2(ng2))
  ALLOCATE(matd(0:nidbas2,0:nidbas2,ng2,kterms,n2))
  DO i2=1,n2
     left2 = i2
     CALL basfun(xg2(:,i2), spl%sp2, fun2, left2)
     DO k=1,kterms
        DO ig2=1,ng2
           DO iw2=0,nidbas2
              DO it2=0,nidbas2
                 matd(it2,iw2,ig2,k,i2) = wg2(ig2,i2) * &
                      &  fun2(it2,idert(k,2),ig2) * &
                      &  fun2(iw2,iderw(k,2),ig2)
              END DO
           END DO
        END DO
     END DO
  END DO
  DEALLOCATE(fun2)
  DEALLOCATE(left2)
!===========================================================================
!              2.0 Assembly loop
!
!   Physical coefficients in Weak form
!
  ALLOCATE(coefs(kterms,ng1,ng2))
  ALLOCATE(matc(ng1,ng2))
!
!   Allocate local matrix A
!
  ALLOCATE(mata(0:nidbas1,0:nidbas1,0:nidbas2,0:nidbas2))
  ALLOCATE(matcd(ng1,0:nidbas2,0:nidbas2))
!
  DO i1=1,n1
     DO i2=1,n2
!
!   Coefficients of the weak form
!
        DO ig1=1,ng1
           DO ig2=1,ng2
              CALL coefeq(xg1(ig1,i1), xg2(ig2,i2), &
                   &      idert, iderw, coefs(:,ig1,ig2))
           END DO
        END DO
!
!   Compute local matrix: A <- E*(C*D^T) + A
!
        mata = 0.0d0
        DO k=1,kterms
!
           matc(1:ng1,1:ng2) = coefs(k,1:ng1,1:ng2)
!
           CALL dgemm('N', 'T', ng1, (nidbas2+1)*(nidbas2+1), ng2, 1.0d0, &
                &      matc, ng1, matd(0,0,1,k,i2), (nidbas2+1)*(nidbas2+1), 0.0d0, &
                &      matcd, ng1)
           CALL dgemm('N', 'N', (nidbas1+1)*(nidbas1+1), (nidbas2+1)*(nidbas2+1), &
                &      ng1, 1.0d0, mate(0,0,1,k,i1), (nidbas1+1)*(nidbas1+1), matcd, ng1, 1.0d0, &
                &      mata, (nidbas1+1)*(nidbas1+1))
!
        END DO
!
!   Map local matrix A to global matrix
!
        DO it1=0,nidbas1
           igt1 = i1+it1; IF(nlper1) igt1 = MODULO(igt1-1,n1) + 1
           DO it2=0,nidbas2
              igt2 = i2+it2; IF(nlper2) igt2 = MODULO(igt2-1, n2) + 1
              irow = igt2 + (igt1-1)*n2
              DO iw1=0,nidbas1
                 igw1 = i1+iw1; IF(nlper1) igw1 = MODULO(igw1-1,n1) + 1
                 DO iw2=0,nidbas2
                    igw2 = i2+iw2; IF(nlper2) igw2 = MODULO(igw2-1, n2) + 1
                    jcol = igw2 + (igw1-1)*n2
                    CALL updtmat(mat, irow, jcol, mata(it1,iw1,it2,iw2))
                 END DO
              END DO
           END DO
        END DO
!
     END DO
  END DO
!===========================================================================
!              9.0  Epilogue
!
  DEALLOCATE(idert, iderw, coefs)
  DEALLOCATE(mata)
  DEALLOCATE(matc)
  DEALLOCATE(matd)
  DEALLOCATE(matcd)
  DEALLOCATE(mate)