!>
!> @file test_gvec1d.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>
!>
!
!    Test implementation of 1D vectors with arbirary
!    vector bounds and ghost cell width.
!
!  T.M. Tran (09/2013)
!
MODULE gvector
  USE iso_fortran_env, ONLY : rkind => real64
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: rkind, gvector_1d, disp, norm2, &
       &    OPERATOR(+), OPERATOR(-), OPERATOR(*), &
       &    ASSIGNMENT(=)

  TYPE gvector_1d
     INTEGER :: s, e, g
     REAL(rkind), ALLOCATABLE :: val(:)
  END TYPE gvector_1d

  INTERFACE gvector_1d
     MODULE PROCEDURE constructor
  END INTERFACE gvector_1d
  INTERFACE OPERATOR(+)
     MODULE PROCEDURE add_scal
     MODULE PROCEDURE add_vec
  END INTERFACE OPERATOR(+)
  INTERFACE OPERATOR(-)
     MODULE PROCEDURE minus_vec
     MODULE PROCEDURE substract_vec
  END INTERFACE OPERATOR(-)
  INTERFACE OPERATOR(*)
     MODULE PROCEDURE scale_left
     MODULE PROCEDURE scale_right
  END INTERFACE OPERATOR(*)
  INTERFACE ASSIGNMENT(=)
     MODULE PROCEDURE from_scal
     MODULE PROCEDURE from_vec
  END INTERFACE ASSIGNMENT(=)
  INTERFACE norm2
     module procedure norm2_gvector_1d
  END INTERFACE norm2

CONTAINS
  FUNCTION constructor(s, e, g) RESULT(res)
    INTEGER, INTENT(in)           :: s, e
    INTEGER, OPTIONAL, INTENT(in) :: g
    TYPE(gvector_1d)              :: res
    INTEGER :: lb, ub
    res%g= 0
    IF(PRESENT(g)) res%g=g
    res%s=s
    res%e=e
    lb = res%s-res%g
    ub = res%e+res%g
    ALLOCATE(res%val(lb:ub))
    res%val = -9999.0
  END FUNCTION constructor

  FUNCTION add_vec(lhs, rhs) RESULT(res)
    TYPE(gvector_1d), INTENT(in) :: lhs, rhs
    TYPE(gvector_1d)             :: res
    res = gvector_1d(lhs%s, lhs%e, lhs%g)
    res%val(res%s:res%e) = lhs%val(res%s:res%e) + rhs%val(res%s:res%e)
  END FUNCTION add_vec

  FUNCTION add_scal(lhs, rhs) RESULT(res)
    TYPE(gvector_1d), INTENT(in) :: lhs
    REAL(rkind), INTENT(in)      :: rhs
    TYPE(gvector_1d)             :: res
    res = gvector_1d(lhs%s, lhs%e, lhs%g)
    res%val(res%s:res%e) = lhs%val(res%s:res%e) + rhs
  END FUNCTION add_scal

  FUNCTION minus_vec(this) RESULT(res)
    TYPE(gvector_1d), INTENT(in) :: this
    TYPE(gvector_1d)             :: res
    res = gvector_1d(this%s, this%e, this%g)
    res%val(res%s:res%e) = -this%val(res%s:res%e)
  END FUNCTION minus_vec

  FUNCTION substract_vec(lhs, rhs) RESULT(res)
    TYPE(gvector_1d), INTENT(in) :: lhs, rhs
    TYPE(gvector_1d)             :: res
    res = gvector_1d(lhs%s, lhs%e, lhs%g)
    res = lhs + (-rhs)
  END FUNCTION substract_vec

  FUNCTION scale_left(lhs, rhs) RESULT(res)
    REAL(rkind), INTENT(in)      :: lhs
    TYPE(gvector_1d), INTENT(in) :: rhs
    TYPE(gvector_1d)             :: res
    res = gvector_1d(rhs%s, rhs%e, rhs%g)
    res%val(res%s:res%e) = lhs * rhs%val(res%s:res%e)
  END FUNCTION scale_left

  FUNCTION scale_right(lhs, rhs) RESULT(res)
    TYPE(gvector_1d), INTENT(in) :: lhs
    REAL(rkind), INTENT(in)      :: rhs
    TYPE(gvector_1d)             :: res
    res = gvector_1d(lhs%s, lhs%e, lhs%g)
    res%val(res%s:res%e) = rhs * lhs%val(res%s:res%e)
  END FUNCTION scale_right

  SUBROUTINE from_vec(lhs, rhs)
    TYPE(gvector_1d), INTENT(inout) :: lhs
    REAL(rkind), INTENT(in)         :: rhs(:)
    INTEGER :: n
    n = lhs%e - lhs%s + 1
    IF(SIZE(rhs) .NE. n) THEN
       PRINT*, 'from_vec: sizes of rhs and lhs not equal!'
       STOP
    END IF
    lhs%val(lhs%s:lhs%e) = rhs(1:n)
  END SUBROUTINE from_vec

  SUBROUTINE from_scal(lhs, rhs)
    TYPE(gvector_1d), INTENT(inout) :: lhs
    REAL(rkind), INTENT(in)         :: rhs
    lhs%val(lhs%s:lhs%e) = rhs
  END SUBROUTINE from_scal

  SUBROUTINE disp(str,this)
    CHARACTER(len=*), INTENT(in) :: str
    TYPE(gvector_1d), INTENT(in) :: this
    WRITE(*,'(/a,3i6)') str//': s, e, g =', this%s, this%e, this%g
    WRITE(*,'(10(1pe12.3))') this%val
  END SUBROUTINE disp

  FUNCTION norm2_gvector_1d(this) RESULT(res)
    TYPE(gvector_1d), INTENT(in) :: this
    REAL(rkind)                  :: res
    res = NORM2(this%val(this%s:this%e))
  END FUNCTION norm2_gvector_1d
END MODULE gvector

PROGRAM main
  USE gvector
  IMPLICIT NONE
  INTEGER :: s=0, e=5, g=1
  INTEGER :: i, lb, ub
  REAL(rkind) :: a=0.1
  TYPE(gvector_1d) :: v1, v2, v3
!
  lb = s-g
  ub = e+g
  v1 = gvector_1d(s, e, g)
  v1%val(s:e) = [ (i, i=s,e) ]
  CALL disp('v1', v1)
!
  v2 = v1 + a*v1
  CALL disp('v1+a*v1', v2)
!
  v3 = v1 - v1*a
  CALL disp('v1-v1*a', v3)
!
  WRITE(*,'(a,1pe12.3)') 'norm of v1      =', NORM2(v1)
  WRITE(*,'(a,1pe12.3)') 'norm of v1-a*v1 =', NORM2(v1-a*v1)
!
  v1 = 0.0d0
  CALL disp('Should be all zero', v1)
  v2 = [ 1.d0, 2.d0, 3.d0, 4.d0, 5.d0, 6.d0 ]
  CALL disp('Should be (1. 2. 3. 4. 5. 6.)', v2)

END PROGRAM main