
module MatlabUtils
 implicit none

 contains
  function kron(X,Y) result(Z)
  implicit none
   complex,intent(in) :: X(:,:)
   complex,intent(in) :: Y(:,:)
   complex,allocatable ::Z(:,:)
   !local varibales
   integer :: m,n
   integer :: p,q
   integer :: mp,nq
   integer :: i,j
   integer :: mpmin,mpmax
   integer :: nqmin,nqmax

   m=size(X,dim=1)
   n=size(X,dim=2)

   p=size(Y,dim=1)
   q=size(Y,dim=2)

   mp=m*p
   nq=n*q
   allocate(Z(mp,nq))
   !forall(i=1:m,j=1:n)
   do j=1,n
     do i=1,m
       mpmin = 1+p*(i-1)
       mpmax = p+p*(i-1)
       nqmin = 1+q*(j-1)
       nqmax = q+q*(j-1)
       Z(mpmin:mpmax,nqmin:nqmax) = X(i,j)*Y
     end do
   end do
  !end forall
  end function kron

  function linspace(ymin,ymax,numY) result(y)
   implicit none
   real,intent(in) :: ymin
   real,intent(in) :: ymax
   integer,intent(in) :: numY
   real,allocatable :: y(:)
   real :: dy
   integer :: ic

   if(numY .gt. 0) then
    allocate(y(numY))
    if(numY .gt. 1) then
     dy = (ymax-ymin)/(numY-1)
    else
     dy = 0
    endif
    y = (/((ymin+(ic-1)*dy),ic=1,numY)/)
   else
    print*,'ErrMessage:numY must be a positive integer'
    stop
   endif
  end function linspace

  function getMat(m,n,val,MatType) result(z)
   ! val = 1.0 for Ones
   ! val = 0.0 for Zeros
   implicit none
   integer,intent(in) :: m
   integer,intent(in),optional :: n
   real,intent(in),optional :: val
   character(len=1),intent(in),optional :: MatType
   integer :: nloc,ic,jc
   real, allocatable  :: z(:,:)
   if (present(n)) then
    allocate(z(m,n))
    nloc = n
   else
    allocate(z(m,m))
    nloc=m
   endif

     if(present(val)) then
      if(present(MatType)) then
       select case(MatType)
        case('I')
         if(m==nloc) then
          do ic=1,m
           z(ic,ic) = val
          end do
         else
          do jc=1,nloc
           do ic=1,m
            if(ic==jc) z(ic,jc)=val
           end do
          end do
         endif
       end select
      else
       z = val
      endif
     else
      z = 0.0
      if(present(MatType)) then
       select case(MatType)
        case('I')
         if(m==nloc) then
          do ic=1,m
           z(ic,ic) = 1.0
          end do
         else
          do jc=1,nloc
           do ic=1,m
            if(ic==jc) z(ic,jc)=1.0
           end do
          end do
         endif
       end select
      endif
     endif

  end function getMat

  function getVec(m,val) result(z)
   ! val = 1.0 for Ones
   ! val = 0.0 for Zeros
   implicit none
   integer,intent(in) :: m
   real,intent(in),optional :: val
   real, allocatable  :: z(:)

    allocate(Z(m))
     if(present(val)) then
       Z = val
     else
       Z = 0.0
     endif
  end function getVec

  function diag(v,k) result(Z)
  implicit none

  real,intent(in) :: v(:)
  integer,intent(in) :: k

  real,allocatable :: Z(:,:)

  integer :: lv,absk
  integer :: m,n
  integer :: ic,jc,icmin,icmax
  lv = size(v,dim=1)
  absk = abs(k)
  m = lv + absk
  n = m
  Z = getMat(m)
    icmin = 1 + (absk-k)/2
    icmax = m - (absk+k)/2
   if(k .ge. 0) then

    !(icmax-icmin = m-1 -abs(k) +1)
    do ic=icmin,icmax
     Z(ic,ic+k) = v(ic)
    end do
   else
    do ic=icmin,icmax
     Z(ic,ic+k) = v(ic-absk)
    end do

   endif
  end function diag
   function trace(Matrix) result(tra)
  complex,dimension(:,:),intent(in) :: Matrix
  complex :: tra
  integer :: n,i
  N = size(Matrix,1)

   tra=(0.0,0.0)
   do i=1,N
    tra = tra + Matrix(i,i)
   end do
 end function
end module MatlabUtils

module qMatsolve

contains
subroutine solveFixedPoint(Epieta,Imat,Alp,beta,numItr,gs)
 use lapack95
 use blas95
 implicit none
 complex,intent(in) :: Epieta ! E + i*eta; i(complex)
 complex,intent(in) :: Imat(:,:) ! Complex Identity matrix
 complex,intent(in) :: Alp(:,:)
 complex, intent(in) :: beta(:,:)
 integer, intent(in) :: numItr
 complex,intent(inout) :: gs(:,:)

 ! local variables
 integer,allocatable :: ipiv(:)
 complex,allocatable :: betaTran_gs_beta(:,:)
 complex,allocatable :: gs_beta(:,:),gsold(:,:)
 real,allocatable :: lwork(:)
 real,external :: clange
 real :: FnormDiffX,FnormResidue
 integer :: N,ic
 integer :: info

 N = size(gs,dim=1)
 allocate(ipiv(N))
 allocate(lwork(N))
 allocate(gs_beta(N,N),gsold(N,N))
 ! initial guess
 gs = Epieta*Imat - Alp
 call getrf(gs,ipiv,info) ! LU factorize
 call getri(gs,ipiv,info) ! invert and overwrite
 gsold = gs
  ! call gemm3m(a, b, c [,transa][,transb] [,alpha][,beta])
 ! C := alpha*op(A)*op(B) + beta*C,
 !
  print*,'iterno:','	  FnormResidue	 ','FnormDiffX	','info'
  do ic=1,numItr
  !fixedpoint=.true.
  !do while(fixedpoint)
  gs_beta=0.0
   call gemm3m(gs,beta,gs_beta,'N','N',(1.0,0.0),(0.0,0.0))
   gs = Epieta*Imat - Alp
   call gemm3m(beta,gs_beta,gs,'C','N',(-1.0,0.0),(1.0,0.0))
   call getrf(gs,ipiv,info) ! LU factorize
   call getri(gs,ipiv,info) ! invert and overwrite
   FnormDiffX=clange('F',N,N,gs-gsold,N,lwork)
   gs = (0.75*gs+0.25*gsold)
   gsold = gs
   gs_beta = matmul(Epieta*Imat - Alp,gs) &
             - matmul(matmul(matmul(transpose(beta),gs),beta),gs)-Imat
   FnormResidue = clange('F',N,N,gs_beta,N,lwork)
   lwork=0.0
   print*,ic,FnormResidue,FnormDiffX,info
  end do

end subroutine solveFixedPoint
subroutine solveQuadMatEqnv2(A,B,C,Xt,sigma,Xinit)
 use lapack95
 use blas95
 use f95_PRECISION, only: WP=>SP
 implicit none
 ! Solves a Quaratic MAtrix Equation
 ! A*X*B*X + C*X+I=0 ! I = Identitiy MAtrix
 ! input variable
 complex(WP),intent(in) :: A(:,:)
 complex(WP),intent(in) :: B(:,:)
 complex(WP),intent(in) :: C(:,:)
 complex(WP),intent(in),optional :: Xinit(:,:)
 ! output variables
 complex(WP),intent(inout) :: Xt(:,:) ! Solution
 complex(WP),intent(inout) :: sigma(:,:) ! Solution
 ! local variables
 real(WP),dimension(:),allocatable :: lwork
 complex(WP),dimension(:,:),allocatable:: lFs,lCs,lDs,lEs,lAs,lBs,lXs
 complex(WP),dimension(:,:),allocatable :: QX
 complex(WP),dimension(:,:),allocatable :: lW1,lZ1
 complex(WP),dimension(:,:),allocatable :: lW2,lZ2
 complex(WP),dimension(:,:),allocatable :: QXZ2,CsZ2,ldXs
 complex(WP),dimension(:),allocatable :: lalpha,lbeta
 complex(WP),dimension(:,:),allocatable :: Xold,X(:,:)
 complex(WP),dimension(:,:),allocatable :: Imat
 complex(WP) :: cOne,czero
 integer(WP) :: info
 integer(WP) :: ic
 integer(WP) :: lmr
 integer(WP) :: ijob=1
 integer(WP),allocatable :: ipiv(:)
 real(WP) :: dif,scaleFac
 real(WP) :: fnormA,fnormB,fnormC,fnormX,fnormQ,fnormX0
 real(WP),external :: clange,second
 real(WP) ::errMaxNewtonExLS,errMaxFixedPoint,errMaxQX,errMaxDiffX
 real :: t,err,fnormdiffX,time1,time2
 real :: errMQ
 integer :: MaxITs
 logical :: fixedPoint,NewtonExLS
 cOne = (1.0,0.0)
 czero = (0.0,0.0)
  lmr = size(A,1)
 !---------- STEP 1 -----------------------
  ! Initialize solver
  allocate(lAs(lmr,lmr),lBs(lmr,lmr),lCs(lmr,lmr))
  allocate(lDs(lmr,lmr),lEs(lmr,lmr),lFs(lmr,lmr))
  allocate(lW1(lmr,lmr),lZ1(lmr,lmr))
  allocate(lW2(lmr,lmr),lZ2(lmr,lmr))

  allocate(Imat(lmr,lmr))
  allocate(QX(lmr,lmr))
  allocate(ipiv(lmr))
  allocate(QXZ2(lmr,lmr),CsZ2(lmr,lmr),ldXs(lmr,lmr))
  allocate(Xold(lmr,lmr),X(lmr,lmr))
  allocate(lwork(lmr))

  Imat=czero
  do ic=1,lmr
   Imat(ic,ic)=(1.0,0.0)
  end do

  errMAxFixedPoint = 1.0E-5
  errMaxDiffX =  1E-4

  ! Intial Guess
  if(present(Xinit)) then
   X = Xinit
  else
   X = -Imat
   lAs = C
   call gesv( lAs,X)
  end if
  Xold = X
 ! print*,X
 ! stop
 ! intialize fixed point
   call gemm3m(A,X,lAs,'N','N',cOne,czero)
   lBs = C
   call gemm3m(lAs,B,lBs,'N','N',cOne,cOne)
 !print*,lBs
 !stop

 print*,'iterno:','	  FnormResidue	 ','FnormDiffX	','info'
 !---------- STEP 2 -----------------------
 ! Fixed point iteration
  fixedpoint = .true.
  ic=0
  do while(fixedpoint)
   ic=ic+1
  ! X = -Imat
   lCs = lBs
  ! call gesv(lCs,X,info=info)
  ! call zgesv( lmr, lmr, lCs, lmr, ipiv, X, lmr, info )
   X = -Imat
    call getrf(lCs,ipiv,info)
    call getrs(lCs, ipiv, X ,info=info )
!   call getri(X,ipiv,info)
   X = 0.75*X + 0.25*Xold
   fnormDiffX = clange('F', lmr,lmr, X-Xold,lmr, lwork )
   Xold = X
   !call gemm3m(A,X,lAs,'N','N',cOne,czero)
   lAs = matmul(A,X)
   !lBs = C
   !call gemm3m(lAs,B,lBs,'N','N',cOne,cOne)
   lBs = matmul(lAs,B) + C
   !QX = Imat
   !call gemm3m(lBs,X,QX,'N','N',cOne,cOne)
   QX = matmul(lBs,X) + Imat
   fnormQ = clange('F', lmr,lmr,QX,lmr, lwork )

   !/clange('F', lmr,lmr, X+Xold,lmr, lwork )
   fixedpoint = (fnormQ > errMaxFixedPoint).and.(ic < 350)
   !.and.(fnormDiffX > errMaxdiffX)
   print*,ic,fnormQ,fnormDiffX,info
  end do
  Xt = X
  sigma = lBs -C

  deallocate(lAs,lBs,lCs)
  deallocate(lDs,lEs,lFs)
  deallocate(lW1,lZ1)
  deallocate(lW2,lZ2)
  deallocate(Imat)
  deallocate(QX)
  deallocate(QXZ2,CsZ2,ldXs)
  deallocate(Xold)
  deallocate(lwork)

end subroutine solveQuadMatEqnv2

end module qMatsolve


module GreensFunctionVar
  use F95_PRECISION, only: WP=>SP

  implicit none
  integer,parameter :: Nu=4
  complex(WP),dimension(:,:),allocatable :: alpU,beta1,beta2,betaW,betaL
  complex(WP),dimension(:,:),allocatable :: betaWTrans,betaLTrans
  complex,dimension(:,:),allocatable :: alpha,beta,betaTrans
  complex,allocatable :: sigmaL(:,:),sigmaR(:,:)
  complex,allocatable :: Aleft(:,:),Cleft(:,:)
  complex,allocatable :: Aright(:,:),Cright(:,:)
  complex,allocatable :: Bleft_right(:,:),XL(:,:),XR(:,:),X(:,:)
  complex,allocatable :: XLOld(:,:),XROld(:,:)
  complex,allocatable :: ImatW(:,:),ImatWu(:,:),ImatWl(:,:) ! Idendity matrices of size NW
  complex,allocatable :: Imat(:,:) ! Identity matrix of size Nw*Nu

  complex :: t
  integer :: Nw,NL
  integer :: Nc
  real :: eta
  integer :: GNRtypeID
  character(len=5) :: GNRARM='GNRARM'
  character(len=5) :: GNRZig='GNRZIG'
  complex,parameter :: cOne=(1.0,0.0)
  complex,parameter :: cZero=(0.0,0.0)

	contains
  subroutine initGreensFunctionVar()
   use MatlabUtils
   implicit none
   integer :: inf
   integer :: i,j

   allocate(alpU(Nu,Nu),betaW(Nu,Nu),betaL(Nu,Nu),betaWTrans(Nu,Nu),betaLTrans(Nu,Nu))
   allocate(beta1(Nu,Nu),beta2(Nu,Nu))
   t=(-1.0,0.0)
   !---------------------------------
   alpU = cZero
   alpU(1,2)= t
   alpU(2,3)= t
   alpU(3,4)= t
   alpU(2,1)= t
   alpU(3,2)= t
   alpU(4,3)= t
   !----------------------------------
   beta1=cZero
   beta1(2,1) = t
   beta1(3,4) = t
   !----------------------------------
   beta2 = cZero
   beta2(4,1) = t
   !----------------------------------
   print*,'Type the Number of Unit Cells, Nw, along Width:'
   print*,' an integer value such that 3 < NW < 10'
   read*,Nw
   !print*,'Type the Number of Unit Cells, NL, along Length:'
   !read*,NL
   NL=1
   print*,'select the configuration of device you want to simulate Armchair/Zigzag\n'
   print*,'[1] for Armchair\n'
   print*,'[2] for Zigzag\n'
   read*,GNRtypeID
   select case(GNRtypeID)
   case(1)
   betaW = beta1
   betaL = beta2
   case(2)
   betaW = beta2
   betaL = beta1
   case default
   betaW = beta1
   betaL = beta2
   end select
   betaWTrans = transpose(betaW)
   betaLTrans = transpose(betaL)

   IMatW = diag(getVec(Nw,1.0),0)*cOne
   IMatWu = diag(getVec(Nw-1,1.0),1)*cOne
   IMatWl = diag(getVec(Nw-1,1.0),-1)*cOne

   alpha = kron(ImatW,alpu) + kron(ImatWu,betaW) + kron(ImatWl,betaWTrans) ! diagonal
   beta  = kron(ImatW,betaL) ! off diagonal
   betaTrans = transpose(conjg(beta))
   !
   eta = 1E-4
   print*,'Enter the Energy eta for the contact selfenergy\n'
   print*,'A real Value such that 0 < eta < 1.0E-1'
   read*,eta

   Nc = Nw*Nu
   Imat = diag(getVec(Nc,1.0),0)*cOne

   Aleft = betaTrans
   Cleft = beta

   Aright = beta
   Cright = betaTrans

   allocate(sigmaL(Nc,Nc),sigmaR(Nc,Nc))
   allocate(XLold(Nc,Nc),XRold(Nc,Nc),X(Nc,Nc))

   print*,'Intialization SucessFull'
  end subroutine initGreensFunctionVar


 elemental function fermi(x) result(fe)
   real,intent(in) :: x
   real :: fe
   fe = 1.0/(1.0 + exp(x))
  end function fermi

 subroutine getContactSelfEnergyv2(Eo,LorR,initGuess,sigmaLR)

 use lapack95
 use blas95
 use qMatsolve
 implicit none
 ! input variables
 real,intent(in) :: Eo
 character,intent(in) :: LorR
 logical,intent(in) :: initGuess
 ! in/output variables
 complex,intent(inout)  :: sigmaLR(:,:)
 ! Local Variables
 complex :: energy
   energy = Eo*(1.0,0.0) + eta*(0.0,1.0)
   Bleft_right = -(energy*Imat - alpha)
   X = 0.0*Imat
   !sigmaLR = 0.0*Imat
   !*Imat
   select case(LorR)
   case('L')
   if(initGuess) then
   call solveQuadMatEqnv2(betaTrans, beta,Bleft_right,X,sigmaLR,XLOld)
   XLOld = X
   else
   call solveQuadMatEqnv2(betaTrans, beta,Bleft_right,X,sigmaLR)
   end if

   case('R')
   if(initGuess) then
   call solveQuadMatEqnv2(beta, betaTrans,Bleft_right,X,sigmaLR,XROld)
   XROld = X
   else
   call solveQuadMatEqnv2(beta, betaTrans,Bleft_right,X,sigmaLR)
   endif

   end select

end subroutine getContactSelfEnergyv2




end module GreensFunctionVar




program TestGF
 ! mkl modules
  use blas95
  use lapack95
  use f95_precision
 ! local modules
  use qMatSolve
  use GreensFunctionVar
  use MatlabUtils
  implicit none

  integer :: numE
  real :: Eo
  real :: Emin,Emax
 ! numE=1

 ! Emin = 0.0! eV
 ! Emax = 0.0 ! eV
 ! Eo = linspace(Emin,Emax,numE)

  call initGreensFunctionVar()
  print*,'Enter the Energy for contact selfenergy -1.0 < E < 1.0 :'
  read*,Eo
  ! call this subroutine to run an alternate implimentation as explined in m file'
  call getContactSelfEnergyv2(Eo,'L',.false.,sigmaL)
  ! Call the following to run an implimentation like MATLAB
  call solveFixedPoint(Eo*Cone + eta*(0.0,1.0),Imat,alpha,beta,350,sigmaL)

end program TestGF
