!>
!> @file conmat.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/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
!
!   In this version s[lines are precalculted
!   (on all n1/n2 intervals
!
  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
  LOGICAL, OPTIONAL :: nat_order   ! Natural ordering for 2d-1d mapping
!
  INTEGER :: n1, nidbas1, ndim1, n1e
  INTEGER :: n2, nidbas2, ndim2, n2e
  INTEGER :: ng1, ng2
  INTEGER :: i1, i2, ig1, ig2
  INTEGER :: igt1, igt2, igw1, igw2, irow, jcol
  INTEGER, ALLOCATABLE :: left1(:), left2(:)
!
  LOGICAL :: nlper1, nlper2, nlnat
!
  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, POINTER :: xg1(:,:), xg2(:,:), wg1(:,:), wg2(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: fun1(:,:,:,:), fun2(:,:,:,:)
  DOUBLE PRECISION, ALLOCATABLE :: mata(:,:,:,:), matc(:,:)
  DOUBLE PRECISION, ALLOCATABLE :: matg(:,:,:), matf(:,:,:), matcg(:,:,:)
!===========================================================================
!              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
!
  n1e = n1+nidbas1  ! Number of elements in 1st coordinate
  n2e = n2+nidbas2  ! Number of elements in 2nd coordinate
  iF(nlper2) n2e = n2
!
!   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)
!
!   Splines on all intervals
!
  kmaxder = 1
  IF(PRESENT(maxder)) kmaxder = maxder(1)
  ALLOCATE(fun1(0:nidbas1,0:kmaxder,ng1,n1))
  ALLOCATE(left1(ng1))
  DO i1=1,n1
     left1 = i1
     CALL basfun(xg1(:,i1), spl%sp1, fun1(:,:,:,i1), left1)
  END DO
  DEALLOCATE(left1)
!
  kmaxder = 1
  IF(PRESENT(maxder)) kmaxder = maxder(2)
  ALLOCATE(fun2(0:nidbas2,0:kmaxder,ng2,n2))
  ALLOCATE(left2(ng2))
  DO i2=1,n2
     left2 = i2
     CALL basfun(xg2(:,i2), spl%sp2, fun2(:,:,:,i2), left2)
  END DO
  DEALLOCATE(left2)
!
!   Ordering in local to global matrix mapping
!
  nlnat = .FALSE.
  IF(PRESENT(nat_order)) nlnat = nat_order
!===========================================================================
!              2.0 Assembly loop
!
!   Weak form
!
  kterms = mat%nterms
  ALLOCATE(idert(kterms,2))
  ALLOCATE(iderw(kterms,2))
  ALLOCATE(coefs(kterms,ng1,ng2))
!
!   Allocate local matrices
!
  ALLOCATE(mata(0:nidbas1,0:nidbas1,0:nidbas2,0:nidbas2))
  ALLOCATE(matc(ng1,ng2))
  ALLOCATE(matg(0:nidbas2,0:nidbas2,ng2))
  ALLOCATE(matf(0:nidbas1,0:nidbas1,ng1))
  ALLOCATE(matcg(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)
!
           DO it1=0,nidbas1
              DO iw1=0,nidbas1
                 DO ig1=1,ng1
                    matf(it1,iw1,ig1) = wg1(ig1,i1) * &
                         &  fun1(it1,idert(k,1),ig1,i1) * &
                         &  fun1(iw1,iderw(k,1),ig1,i1)
                 END DO
              END DO
           END DO
!
           DO it2=0,nidbas2
              DO iw2=0,nidbas2
                 DO ig2=1,ng2
                    matg(it2,iw2,ig2) = wg2(ig2,i2) * &
                         &  fun2(it2,idert(k,2),ig2,i2) * &
                         &  fun2(iw2,iderw(k,2),ig2,i2)
                 END DO
              END DO
           END DO
!
           CALL dgemm('N', 'T', ng1, (nidbas2+1)*(nidbas2+1), ng2, 1.0d0, &
                &      matc, ng1, matg, (nidbas2+1)*(nidbas2+1), 0.0d0, &
                &      matcg, ng1)
           CALL dgemm('N', 'N', (nidbas1+1)*(nidbas1+1), (nidbas2+1)*(nidbas2+1), &
                &      ng1, 1.0d0, matf, (nidbas1+1)*(nidbas1+1), matcg, 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 = glmap(igt1, igt2, n1e, n2e)
              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 = glmap(igw1, igw2, n1e, n2e)
                    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(fun1)
  DEALLOCATE(fun2)
  DEALLOCATE(idert, iderw, coefs)
  DEALLOCATE(mata)
  DEALLOCATE(matc)
  DEALLOCATE(matg)
  DEALLOCATE(matcg)
  DEALLOCATE(matf)
!
CONTAINS
  INTEGER FUNCTION glmap(i,j,n1,n2)
    INTEGER, INTENT(in) :: i,j,n1,n2
    IF(nlnat) THEN
       glmap = (j-1)*n1 + i
    ELSE
       glmap = (i-1)*n2 + j
    END IF
  END FUNCTION glmap