!include 'mkl_pardiso.f90'
module d1isoModel
    use kinddef
    !USE mkl_pardiso
    implicit none
    type :: d1isoGrid
        private
        real(rk), dimension(:), allocatable :: z
        real(rk), dimension(:), allocatable :: sigma
        integer(ik):: n
    contains
        procedure, public :: init => initd1isoModel
!        procedure, public :: printOnScreen => printd1isoModelOnScreen
        procedure, public :: computeEx
        procedure, public :: deallocate => clean_d1isoGrid
    end type d1isoGrid
contains
 subroutine initd1isoModel(this,n,z,sigma)
        class(d1isoGrid) :: this
        integer(ik) :: n
        real(rk), dimension(n) :: z
        real(rk), dimension(n-1) :: sigma
        allocate(this%z(n))
        allocate(this%sigma(n-1))

        this%z = z
        this%sigma = sigma
        this%n = n
    end subroutine initd1isoModel

    subroutine clean_d1isoGrid(this)
        class(d1isoGrid) :: this
        if (allocated(this%sigma))then
            deallocate(this%sigma)
        end if

        if (allocated(this%z))then
            deallocate(this%z)
        end if
        print*, 'memory used for d1isoGrid relased ...'
    end subroutine clean_d1isoGrid

 function computeEx(this,freq,Ex) result(info)

      use physicsparameters
        class(d1isoGrid), intent(in) :: this
        real(rk), intent(in) :: freq
        complex(ck), dimension(:), allocatable :: Ex
        integer :: info
        ! pardiso header file, parameters and variables, based on the intel example
        !TYPE(MKL_PARDISO_HANDLE), dimension(64)  :: pt
        integer(8), dimension(64) :: pt
        INTEGER :: iparm(64)
        INTEGER :: i, idum
        complex(rk) :: ddum!, waltime1, waltime2
        INTEGER :: maxfct, mnum, mtype, phase, nPardiso, nrhs, error, msglvl


        ! local arguments!
        real(rk) :: omega
        real(rk) :: znodes(this%n-1), C(this%n-3)
        integer  :: j
!        integer (ik) ::  icol, irow
        integer  :: n
        integer  :: elementIndex
        real(rk) :: zGhostNode
        complex (ck), dimension ((this%n -3)*3 + 2 + 1) :: values
        integer , dimension ((this%n -3)*3 + 2 + 1) :: columns
        integer , dimension (this%n) :: rowIndex
        complex (ck), dimension(this%n-1) :: b
        integer, dimension(this%n -1) :: perm
        n = this%n
        nPardiso = n-1

        ! compute znodes coordiantes
        do i = 2,n
            znodes(i-1) = (this%z(i) + this%z(i-1))*0.5_rk
        end do
        omega = 2*pi*freq
        allocate (Ex(this%n -1)) ! columns store respectively Ex, Ey and Ez.

        b = cmplx(0._rk,0._rk,kind=rk)

        b(1) = omega*mu0*i1

        C = computeC(n,this%z,znodes)
        !print*, znodes
        !print*, C

        zGhostNode = znodes(1) - (znodes(2) - znodes(1))
        !        print*, zGhostNode, znodes

        print*, 'forming the finite-differences sparse matrix and r.h.s. ...'
        values(1) = -((znodes(1) - zGhostNode)**(-1))
        columns(1) = 1
        rowIndex(1) = 1
        values(2) = ((znodes(1) - zGhostNode)**(-1))
        columns(2) = 2
        elementIndex = 2
        !print*, 'n = ', n

        do i = 2,n-2
            do j = -1,1
                elementIndex = elementIndex+1
                select case (j)
                    case(-1)
                        values(elementIndex) = C(i-1)*(znodes(i+1)-znodes(i))
                        !print*, c(i), values(elementIndex)
                        columns(elementIndex) = i+j
                        rowIndex(i) = elementIndex
                    case(0)
                        values(elementIndex) = -C(i-1)*(znodes(i+1)-znodes(i-1))
                        columns(elementIndex) = i+j
                    case(1)
                        values(elementIndex) = C(i-1)*(znodes(i)-znodes(i-1))
                        columns(elementIndex) = i+j
                end select
            end do
        end do
        elementIndex = elementIndex+1
        values(elementIndex) = 1._rk
        columns(elementIndex) = n-1
        rowIndex(n-1) = elementIndex
        rowIndex(n) = size(values)+1


        print*, values
        print*, columns
        print*, rowIndex

        ! calling the sparse solver
        !C..
        !C.. Set up PARDISO control parameter
        !C..
        do i = 1, 64
            iparm(i) = 0
        end do
!!$        iparm(1) = 0 ! no solver default
!!$        iparm(2) = 2 ! uses openMP
!!$        iparm(3) = 0 ! numbers of processors
!!$        iparm(4) = 62 ! no iterative-direct algorithm
!!$        iparm(5) = 0 ! no user fill-in reducing permutation
!!$        iparm(6) = 0 ! =0 solution on the first n compoments of x
!!$        iparm(7) = 0 ! not in use
!!$        iparm(8) = 0 ! numbers of iterative refinement steps
!!$        iparm(9) = 0 ! not in use
!!$        iparm(10) = 13 ! perturbe the pivot elements with 1E-13
!!$        iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
!!$        iparm(12) = 0 ! not in use
!!$        iparm(13) = 1 ! not in use
!!$        iparm(14) = 0 ! Output: number of perturbed pivots
!!$        iparm(15) = 0 ! not in use
!!$        iparm(16) = 0 ! not in use
!!$        iparm(17) = 0 ! not in use
!!$        iparm(18) = -1 ! Output: number of nonzeros in the factor LU
!!$        iparm(19) = -1 ! Output: Mflops for LU factorization
!!$        iparm(20) = 0 ! Output: Numbers of CG Iterations
        iparm(27) = 1
        error = 0 ! initialize error flag
        msglvl = 1 ! print statistical information
        mtype = 13 ! complex unsymmetric
        !C.. Initiliaze the internal solver memory pointer. This is only
        !C necessary for the FIRST call of the PARDISO solver.
        do i = 1, 64
!           pt( i )%DUMMY =  0
           pt(i) = 0
        end do

        maxfct = 1
        mnum = 1
        nrhs = 1
        !C.. Reordering and Symbolic Factorization, This step also allocates
        !C all memory that is necessary for the factorization
        phase = 11 ! only reordering and symbolic factorization
        print*, ' calling sparse solver'

        CALL pardiso (pt, maxfct, mnum, mtype, phase, nPardiso, values, rowIndex, columns, &
                    idum, nrhs, iparm, msglvl, ddum, ddum, error)
        WRITE(*,*) 'Reordering completed ... '
        IF (error .NE. 0) THEN
            WRITE(*,*) 'The following ERROR was detected: ', error
            STOP
        END IF
        WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
        WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
        !C.. Factorization.
        phase = 22 ! only factorization
        CALL pardiso (pt, maxfct, mnum, mtype, phase, nPardiso, values, rowIndex, columns, &
                   idum, nrhs, iparm, msglvl, ddum, ddum, error)
        WRITE(*,*) 'Factorization completed ... '
        IF (error .NE. 0) THEN
            WRITE(*,*) 'The following ERROR was detected: ', error
            STOP
        ENDIF
        !C.. Back substitution and iterative refinement
        iparm(8) = 2 ! max numbers of iterative refinement steps
        phase = 33 ! only factorization

        CALL pardiso (pt, maxfct, mnum, mtype, phase, nPardiso, values, rowIndex, columns, &
                   idum, nrhs, iparm, msglvl, b, Ex, error)
        WRITE(*,*) 'Solve completed ... '
        WRITE(*,*) 'The solution of the system is '
        DO i = 1, nPardiso
            WRITE(*,*) ' Ex(',i,') = ', Ex(i)
        END DO
        !C.. Termination and release of memory
        phase = -1 ! release internal memory
        CALL pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
                   idum, nrhs, iparm, msglvl, ddum, ddum, error)
        info = error
        return
    end function computeEx

    function computeC(n,z,znodes) result(C)
        integer(ik), intent(in) :: n
        real(rk), dimension(n) :: z
        real(rk), dimension(n-1) :: znodes
        real(rk), dimension(2:n-2) :: C
        integer(ik) :: i
        !C = 0._rk
        !print*, 'c ='
        !print*, C

        do i = 2,n-2
            C(i) = (z(i+1)-z(i))**(-1) * ((znodes(i+1)-znodes(i))*(znodes(i)-znodes(i-1)))**(-1)
        end do
        !print*, 'c ='
        !print*, C
        return
    end function computeC

end module d1isoModel
