!
! This file is part of SACAMOS, State of the Art CAble MOdels in 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-2017 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 .
!
! 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 .
!
! The University of Nottingham can be contacted at: ggiemr@nottingham.ac.uk
!
! SUBROUTINE solve_real_symm(n, b, x,tol)
! SUBROUTINE dot ( n, x, y ,res)
! SUBROUTINE solve_complex_symm(n, b, x,tol)
! 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)
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
! 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,2*n
! 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
! finish
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)
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
! 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,2*n
! 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
! finish
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