!
! This file is part of SACAMOS, State of the Art CAble MOdels for Spice. 
! It was developed by the University of Nottingham and the Netherlands Aerospace 
! Centre (NLR) for ESA under contract number 4000112765/14/NL/HK.
! 
! Copyright (C) 2016-2018 University of Nottingham
! 
! SACAMOS is free software: you can redistribute it and/or modify it under the 
! terms of the GNU General Public License as published by the Free Software 
! Foundation, either version 3 of the License, or (at your option) any later 
! version.
! 
! SACAMOS 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.
! 
! A copy of the GNU General Public License version 3 can be found in the 
! file GNU_GPL_v3 in the root or at <http://www.gnu.org/licenses/>.
! 
! SACAMOS uses the EISPACK library (in /SRC/EISPACK). EISPACK is subject to 
! the GNU Lesser General Public License. A copy of the GNU Lesser General Public 
! License version can be found in the file GNU_LGPL in the root of EISPACK 
! (/SRC/EISPACK ) or at <http://www.gnu.org/licenses/>.
! 
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
!
! SUBROUTINE solve_real_symm(n, b, x,tol,itmax)
! SUBROUTINE dot ( n, x, y ,res)
! SUBROUTINE solve_complex_symm(n, b, x,tol,itmax)
! SUBROUTINE zdot ( n, x, y ,res)
!
! NAME
!    solve_real_symm
!
! DESCRIPTION
!
! Solve the system Ax=b using a conjugate gradient method. 
! The iterative solution halts when the norm of the residual is less than tol
!
! HISTORY
!
!     started 19/6/2018 CJS
!
! COMMENTS
! 
SUBROUTINE solve_real_symm(n, b, x,tol,itmax)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

integer,intent(IN) :: n

real(dp),intent(IN) :: b(n)
real(dp),intent(INOUT) :: x(n)   ! x contains the initial guess for x
real(dp),intent(IN) :: tol
integer,intent(IN) :: itmax

! local variables

real(dp) :: ak
real(dp) :: rk(n),rkb(n)
real(dp) :: pk(n),pkb(n)
real(dp) :: bk

real(dp) :: t1(n),t2(n)
real(dp) :: res(n),res_norm
real(dp) :: num,den

integer :: i,ii,itloop

logical :: minres

! START

write(*,*)'CALLED solve_real_symm with n=',n

! choice between normal conjugate gradient and minimum residual forms.

minres=.TRUE.

! set starting values
    
! calculate the residual res=b-A.xk
  CALL Aprod( n, x, t1 )    
  res(:)=b(:)-t1(:)

! calculate the size of the residual
  res_norm=0d0
  do i=1,n
    res_norm=res_norm+res(i)*res(i)
  end do
  res_norm=sqrt(res_norm)
    
  write(*,*)'iteration ',0,' norm=',res_norm
  
! initial r1=b-A.x
  CALL Aprod( n, x, t1 )    
  rk(:)=b(:)-t1(:) 
  
  if (minres) then
! For the minimum residual algorithm r1b=A.r1
    CALL Aprod( n, rk, rkb )    
  else
! For the normal algorithm  r1b=r1
    rkb(:)=rk(:)    
  end if

! p1=r1
  pk=rk
  
! p1b=r1b
  pkb=rkb

! iteration loop
  do itloop=1,itmax
  
! calculate ak=(rkb.rk)/(pkb.A.pk)
    
    CALL dot(n,rkb,rk,num)

    CALL Aprod( n, pk, t1 )    
    CALL dot(n,pkb,t1,den)
    
    ak=num/den    
       
! calculate the next value of the solution vector, xk+1=xk+ak*pk
    x(:)=x(:)+ak*pk(:)
    
! calculate the residual res=b-A.xk
    CALL Aprod( n, x, t2 )    
    res(:)=b(:)-t2(:)

! calculate the size of the residual
    res_norm=0d0
    do i=1,n
      res_norm=res_norm+res(i)*res(i)
    end do
    res_norm=sqrt(res_norm)
    
    if (res_norm.LT.tol) then
      write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak
      RETURN
    end if
    
! calculate rk+1=rk-ak*A.pk
    rk(:) =rk(:) -ak*t1(:)
   
! calculate rk+1b=rkb-ak*AT.pkb, note for the symmetric matrices here AT=A
    CALL ATprod( n, pkb, t2 )    
    rkb(:)=rkb(:)-ak*t2(:)
    
! calculate bk=(rk+1b.rk+1)/(rkb.rk)
    den=num
    CALL dot(n,rkb,rk,num)
    bk=num/den
    
    write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak,' bk=',bk
        
! calculate pk+1=rk+1-bk*pk
    pk(:)=rk(:)+bk*pk(:)
    
! calculate pk+1b=rk+1b-bk*pkb
    pkb(:)=rkb(:)+bk*pkb(:)

  end do ! next iteration
  
! We only get here if the iterative solution has not converged after itmax iterations

  run_status='ERROR in solve_real_symm: maximum number of iterations exceeded'
  CALL write_program_status()
  STOP 1

  RETURN


END SUBROUTINE solve_real_symm
!
! _________________________________________________________________________
!
!

SUBROUTINE dot ( n, x, y ,res)

! form the dot product of x and y

USE type_specifications

IMPLICIT NONE

integer n
real(dp) :: x(n), y(n), res

integer :: i

! START

res=0d0

do i=1,n

  res=res+x(i)*y(i)
  
end do

RETURN

END SUBROUTINE dot
!
! NAME
!    solve_complex_symm
!
! DESCRIPTION
!
! Solve the system Ax=b using a conjugate gradient method. 
! The iterative solution halts when the norm of the residual is less than tol
!
! HISTORY
!
!     started 19/6/2018 CJS
!
! COMMENTS
! 
 SUBROUTINE solve_complex_symm(n, b, x,tol,itmax)

USE type_specifications

IMPLICIT NONE

! variables passed to subroutine

integer,intent(IN) :: n

complex(dp),intent(IN) :: b(n)
complex(dp),intent(INOUT) :: x(n)   ! x contains the initial guess for x
real(dp),intent(IN) :: tol
integer,intent(IN) :: itmax

! local variables

complex(dp) :: ak
complex(dp) :: rk(n),rkb(n)
complex(dp) :: pk(n),pkb(n)
complex(dp) :: bk

complex(dp) :: t1(n),t2(n)
complex(dp) :: res(n)
real(dp)    :: res_norm
complex(dp) :: num,den

integer :: i,ii,itloop

logical :: minres

! START

write(*,*)'CALLED solve_real_symm with n=',n

! choice between normal conjugate gradient and minimum residual forms.

minres=.TRUE.

! set starting values
    
! calculate the residual res=b-A.xk
  CALL zAprod( n, x, t1 )    
  res(:)=b(:)-t1(:)

! calculate the size of the residual
  res_norm=0d0
  do i=1,n
    res_norm=res_norm+res(i)*res(i)
  end do
  res_norm=sqrt(res_norm)
    
  write(*,*)'iteration ',0,' norm=',res_norm
  
! initial r1=b-A.x
  CALL zAprod( n, x, t1 )    
  rk(:)=b(:)-t1(:) 
  
  if (minres) then
! For the minimum residual algorithm r1b=A.r1
    CALL zAprod( n, rk, rkb )    
  else
! For the normal algorithm  r1b=r1
    rkb(:)=rk(:)    
  end if

! p1=r1
  pk=rk
  
! p1b=r1b
  pkb=rkb

! iteration loop
  do itloop=1,itmax
  
! calculate ak=(rkb.rk)/(pkb.A.pk)
    
    CALL zdot(n,rkb,rk,num)

    CALL zAprod( n, pk, t1 )    
    CALL zdot(n,pkb,t1,den)
    
    ak=num/den    
       
! calculate the next value of the solution vector, xk+1=xk+ak*pk
    x(:)=x(:)+ak*pk(:)
    
! calculate the residual res=b-A.xk
    CALL zAprod( n, x, t2 )    
    res(:)=b(:)-t2(:)

! calculate the size of the residual
    res_norm=0d0
    do i=1,n
      res_norm=res_norm+abs(res(i))**2
    end do
    res_norm=sqrt(res_norm)
    
    if (res_norm.LT.tol) then
      write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak
      RETURN
    end if
    
! calculate rk+1=rk-ak*A.pk
    rk(:) =rk(:) -ak*t1(:)
   
! calculate rk+1b=rkb-ak*AT.pkb, note for the symmetric matrices here AT=A
    CALL zATprod( n, pkb, t2 )    
    rkb(:)=rkb(:)-ak*t2(:)
    
! calculate bk=(rk+1b.rk+1)/(rkb.rk)
    den=num
    CALL zdot(n,rkb,rk,num)
    bk=num/den
    
    write(*,*)'iteration ',itloop,' norm=',res_norm,' ak=',ak,' bk=',bk
        
! calculate pk+1=rk+1-bk*pk
    pk(:)=rk(:)+bk*pk(:)
    
! calculate pk+1b=rk+1b-bk*pkb
    pkb(:)=rkb(:)+bk*pkb(:)

  end do ! next iteration
  
! We only get here if the iterative solution has not converged after itmax iterations
 
  run_status='ERROR in solve_complex_symm: maximum number of iterations exceeded'
  CALL write_program_status()
  STOP 1

  RETURN


END SUBROUTINE solve_complex_symm
!
! _________________________________________________________________________
!
!

SUBROUTINE zdot ( n, x, y ,res)

! form the dot product of x and y

USE type_specifications

IMPLICIT NONE

integer n
complex(dp) :: x(n), y(n), res

integer :: i

! START

res=0d0

do i=1,n

  res=res+x(i)*y(i)     ! note no complex conjugate in the inner product here. This gives the biconjugate gradient method.
  
end do

RETURN

END SUBROUTINE zdot