scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
quantum-espresso/GWW/gww/fit_multipole.f90

974 lines
22 KiB
Fortran
Raw Blame History

!
! Copyright (C) 2001-2013 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!this contains routines to fit a function on the positive part
!of the imaginary axes with a multipole expansion
MODULE global_minpack
!this module conatins global variables(sigh!) for using old FORTRAN77
! minpack routine
USE kinds, ONLY : DP
IMPLICIT NONE
SAVE
INTEGER, PARAMETER :: maxm=400!max number of samples
INTEGER, PARAMETER :: maxpole=30
INTEGER :: n_poles
COMPLEX(kind=DP) :: c_target(maxm)
REAL(kind=DP) :: freq(maxm)
END MODULE
SUBROUTINE fit_multipole(n,m,z,s,a_0,a,b,delta,thres,maxiter)
!fits with the function f(z)=a_0+\sum_{i=1,m} a_i/(z-b_i)
!the values z_j,s_j
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
implicit none
INTEGER, INTENT(in) :: n!numer of sampled values
INTEGER, INTENT(in) :: m!number of parameters a
COMPLEX(kind=DP), INTENT(in) :: z(n)!where
COMPLEX(kind=DP), INTENT(in) :: s(n)!values s(z_j) to be fitted
COMPLEX(kind=DP), INTENT(inout) :: a_0
COMPLEX(kind=DP), INTENT(inout) :: a(m)
COMPLEX(kind=DP), INTENT(inout) :: b(m)
REAL(kind=DP), INTENT(in) :: delta!parameter for steepest descend
REAL(kind=DP), INTENT(in) :: thres!threshold for convergence
INTEGER, INTENT(in) :: maxiter!maximum number of iterations
REAL(kind=DP) :: chi0, chi1
INTEGER :: i,j,it
COMPLEX(kind=DP) :: cc, grad
COMPLEX(kind=DP) :: new_a_0, new_a(m), new_b(m), old_b(m), old_a(m), old_a_0
REAL(kind=DP) :: rr,rc
REAL(kind=DP) :: dd,ddb
LOGICAL :: random
INTEGER :: ip, im
ddb= 0.01d0
dd = delta
random=.true.
ip=1
im=1
!calculates initial chi
chi0=0.d0
do i=1,n
cc = func(z(i))-s(i)
!cc = (func(z(i))-s(i))/s(i)
chi0=chi0+cc*conjg(cc)
enddo
write(stdout,*) 'a_0', a_0!ATTENZIONE
write(stdout,*) 'a', a
write(stdout,*) 'b', b
write(stdout,*) 'z,s' , z(1),s(1),func(z(1))
write(stdout,*) 'z,s' , z(n),s(n),func(z(n))
do it=1,maxiter
!updates a_0
grad=(0.d0,0.d0)
do i=1,n
grad=grad+(func(z(i))-s(i))
enddo
new_a_0=a_0-dd*grad
if(it==1) write(stdout,*) 'Grad a_0', grad!ATTENZIONE
!updates a(:)
do j=1,m
grad=(0.d0,0.d0)
do i=1,n
grad=grad+(func(z(i))-s(i))/(conjg(z(i))-conjg(b(j)))
enddo
new_a(j)=a(j)-grad*dd
if(it==1) write(stdout,*) 'Grad a', grad!ATTENZIONE
enddo
!updates b(:)
if(.not.random) then
do j=1,m
grad=(0.d0,0.d0)
do i=1,n
grad=grad+(func(z(i))-s(i))*conjg(a(j))/((conjg(z(i))-conjg(b(j)))*(conjg(z(i))-conjg(b(j))))
enddo
new_b(j)=b(j)-grad*dd
if(it==1) write(stdout,*) 'Grad b', grad!ATTENZIONE
enddo
endif
!calculates new chi
old_a_0=a_0
a_0=new_a_0
old_a(:)=a(:)
a(:)=new_a(:)
if(.not. random) b(:)=new_b(:)
chi1=0.d0
do i=1,n
cc = func(z(i))-s(i)
!cc = (func(z(i))-s(i))/s(i)
chi1=chi1+cc*conjg(cc)
enddo
if(chi1 > chi0) then
a_0=old_a_0
a(:)=old_a(:)
write(stdout,*) 'Routine fit_multipole: chi1 > chi0 '
!return
dd=dd*0.1d0
endif
! if((chi0-chi1) <= thres) then
! write(stdout,*) 'Reached threshold', chi0
! return
! endif
chi0=chi1
if(random) then!minimize b in a random fashion
! if(mod(it,50000) == 1) ddb=ddb*0.1d0
if(mod(ip,10)==0) then
ddb=ddb*0.1d0
ip=1
write(stdout,*) 'Random plus'
endif
if(mod(im,100)==0) then
ddb=ddb*10.d0
im=1
write(stdout,*) 'Random minus'
endif
do j=1,m
old_b(:)=b(:)
call random_number(rr)
call random_number(rc)
b(j)=b(j)+ddb*cmplx(rr,rc)
if(aimag(b(j)) >= 0.d0) b(j)=cmplx(real(b(j)),-aimag(b(j)))
chi1=0.d0
do i=1,n
cc = func(z(i))-s(i)
!cc = (func(z(i))-s(i))/s(i)
chi1=chi1+cc*conjg(cc)
enddo
if(chi1<chi0) then
! write(*,*) 'PASSED'
chi0=chi1
ip=ip+1
im=1
else
b(:)=old_b(:)
ip=1
im=im+1
endif
enddo
endif
!write(stdout,*) 'chi0',chi0!ATTENZIONE
enddo
write(stdout,*) 'Routine fit_multipole: maxiter reached ', chi0
return
CONTAINS
FUNCTION func(zz)
COMPLEX(kind=DP) :: func
COMPLEX(kind=DP) :: zz
INTEGER :: ii
func=a_0
do ii=1,m
func=func+a(ii)/(zz-b(ii))
enddo
return
END FUNCTION func
END SUBROUTINE fit_multipole
SUBROUTINE fit_multipole_verlet(n,m,z,s,a_0,a,b,thres,maxiter,ma_0,ma,mb,dt,frice)
!fits with the function f(z)=a_0+\sum_{i=1,m} a_i/(z-b_i)
!the values z_j,s_j
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
implicit none
INTEGER, INTENT(in) :: n!numer of sampled values
INTEGER, INTENT(in) :: m!number of parameters a
COMPLEX(kind=DP), INTENT(in) :: z(n)!where
COMPLEX(kind=DP), INTENT(in) :: s(n)!values s(z_j) to be fitted
COMPLEX(kind=DP), INTENT(inout) :: a_0
COMPLEX(kind=DP), INTENT(inout) :: a(m)
COMPLEX(kind=DP), INTENT(inout) :: b(m)
REAL(kind=DP), INTENT(in) :: thres!threshold for convergence
INTEGER, INTENT(in) :: maxiter!maximum number of iterations
REAL(kind=DP) :: ma_0!fictitious mass relative to a_0
REAL(kind=DP) :: ma!fictitious mass relative to a's
REAL(kind=DP) :: mb!fictitious mass relative to b's
REAL(kind=DP) :: dt!
REAL(kind=DP) :: frice
REAL(kind=DP) :: chi0, chi1
INTEGER :: i,j,it
COMPLEX(kind=DP) :: cc
COMPLEX(kind=DP) :: a_0p,a_0m,ap(m),am(m),bp(m),bm(m)
COMPLEX(kind=DP) :: va_0,va(m),vb(m),fa_0,fa(m),fb(m),va_0m,vam(m),vbm(m),va_0p,vap(m),vbp(m)
!calculates initial chi
chi0=0.d0
do i=1,n
cc = func(z(i))-s(i)
chi0=chi0+cc*conjg(cc)
enddo
write(*,*) 'a_0', a_0!ATTENZIONE
write(*,*) 'a', a
write(*,*) 'b', b
write(*,*) 'z,s' , z(1),s(1),func(z(1))
write(*,*) 'z,s' , z(n),s(n),func(z(n))
write(*,*) 'M a_0:',ma_0
write(*,*) 'M a:',ma
write(*,*) 'M b:',mb
write(*,*) 'Dt:', dt
write(*,*) 'frice:', frice
!first SD step
!updates a_0
fa_0=(0.d0,0.d0)
do i=1,n
fa_0=fa_0-(func(z(i))-s(i))
enddo
!updates a(:)
do j=1,m
fa(j)=(0.d0,0.d0)
do i=1,n
fa(j)=fa(j)-(func(z(i))-s(i))/(conjg(z(i))-conjg(b(j)))
enddo
enddo
!updates b(:)
do j=1,m
fb(j)=(0.d0,0.d0)
do i=1,n
fb(j)=fb(j)-(func(z(i))-s(i))*conjg(a(j))/((conjg(z(i))-conjg(b(j)))*(conjg(z(i))-conjg(b(j))))
enddo
enddo
!step
a_0p=a_0+0.5d0*fa_0*dt*dt/ma_0
ap(:)=a(:)+0.5d0*fa(:)*(dt)**2.d0/ma
bp(:)=b(:)+0.5d0*fb(:)*(dt)**2.d0/mb
write(*,*) 'fa_0', fa_0
write(*,*) 'fa', fa
write(*,*) 'fb', fb
a_0m=a_0
a_0=a_0p
am(:)=a(:)
a=ap(:)
bm(:)=b(:)
b=bp(:)
write(*,*) 'a_0', a_0,ma_0
write(*,*) 'a', a
write(*,*) 'b', b
!set intial velocities
va_0m=(0.d0,0.d0)
vam(:)=(0.d0,0.d0)
vbm(:)=(0.d0,0.d0)
va_0=fa_0*dt/ma_0
va(:)=fa(:)*dt/ma
vb(:)=fb(:)*dt/mb
do it=1,maxiter
!calculates forces
!updates a_0
fa_0=(0.d0,0.d0)
do i=1,n
fa_0=fa_0-(func(z(i))-s(i))
enddo
fa_0=(0.d0,0.d0)!ATTENZIONE
!updates a(:)
do j=1,m
fa(j)=(0.d0,0.d0)
do i=1,n
fa(j)=fa(j)-(func(z(i))-s(i))/(conjg(z(i))-conjg(b(j)))
enddo
enddo
!updates b(:)
do j=1,m
fb(j)=(0.d0,0.d0)
do i=1,n
fb(j)=fb(j)-(func(z(i))-s(i))*conjg(a(j))/((conjg(z(i))-conjg(b(j)))*(conjg(z(i))-conjg(b(j))))
enddo
enddo
!add of frice
!write(*,*) 'fa_0', fa_0
! write(*,*) 'fa', fa
! write(*,*) 'fb', fb
fa_0=fa_0-frice*va_0
fa(:)=fa(:)-frice*va(:)
fb(:)=fb(:)-frice*vb(:)
! write(*,*) 'fa_0', fa_0
! write(*,*) 'fa', fa
! write(*,*) 'fb', fb
!update positions
a_0p=2.d0*a_0-a_0m+(fa_0/ma_0)*(dt)**2.d0
ap(:)=2.d0*a(:)-am(:)+(fa(:)/ma)*(dt)**2.d0
bp(:)=2.d0*b(:)-bm(:)+(fb(:)/mb)*(dt)**2.d0
a_0m=a_0
am(:)=a(:)
bm(:)=b(:)
a_0=a_0p
a(:)=ap(:)
b(:)=bp(:)
!calculates new chi
chi1=0.d0
do i=1,n
cc = func(z(i))-s(i)
chi1=chi1+cc*conjg(cc)
enddo
!calculates velocities but at t-dt/2 ...
va_0=(a_0-a_0m)/dt
va(:)=(a(:)-am(:))/dt
vb(:)=(b(:)-bm(:))/dt
! if((chi0-chi1) <= thres) return
chi0=chi1
if(mod(it,100)==1) write(*,*) chi0!ATTENZIONE
enddo
write(stdout,*) 'Routine fit_multipole: maxiter reached '
return
CONTAINS
FUNCTION func(zz)
COMPLEX(kind=DP) :: func
COMPLEX(kind=DP) :: zz
INTEGER :: ii
func=a_0
do ii=1,m
func=func+a(ii)/(zz-b(ii))
enddo
return
END FUNCTION func
END SUBROUTINE fit_multipole_verlet
!the following routines are for MINPACK
SUBROUTINE fcn_set(m, npoles, omegas,ctarget)
!added internal parmeters to set up function
!the parameters are in the order re(a_0),im(a_0),re(a_1)...re(a_n),im(a_1)..im(a_n),re(b_1)..re(b_n), im(b_1)..im(b_n)
use kinds, ONLY : DP
use io_global, ONLY : stdout
use global_minpack
implicit none
INTEGER :: m !number of variables
INTEGER :: npoles
REAL(kind=DP) :: omegas(m)
COMPLEX(kind=DP) :: ctarget(m)
n_poles=npoles
freq(1:m)=omegas(1:m)
c_target(1:m)=ctarget(1:m)
return
END SUBROUTINE
SUBROUTINE fcn(m,n,x,fvec,iflag)
!added internal parmeters to set up function
!the parameters are in the order re(a_0),im(a_0),re(a_1)...re(a_n),im(a_1)..im(a_n),re(b_1)..re(b_n), im(b_1)..im(b_n)
use kinds, ONLY : DP
use io_global, ONLY : stdout
use global_minpack
implicit none
INTEGER :: m !number of variables
INTEGER :: n!total number of parameters
REAL(kind=DP) :: x(n)!parameters
REAL(kind=DP) :: fvec(m)!evaluated error function
INTEGER :: iflag!not used
COMPLEX(kind=DP) a_0, a(maxpole), b(maxpole)
INTEGER i,j
COMPLEX(kind=DP) :: func,zz
if(m>maxm) then
write(stdout,*) 'FCN: MAXN TOO SMALL'
stop
endif
!check number of parameters
if(n /= 4*n_poles+2) then
write(stdout,*) 'FCN: WRONG NUMBER OF PARAMETERS',n,n_poles
stop
endif
if(n_poles>maxpole) then
write(stdout,*) 'FCN: MAXPOLE TOO SMALL'
stop
endif
!set up parameters
a_0=cmplx(x(1),x(2))
do i=1,n_poles
a(i)=cmplx(x(i*2+1),x(i*2+2))
enddo
do i=1,n_poles
!b(i)=cmplx(x((i+n_poles)*2+1),-(x((i+n_poles)*2+2))**2.d0)
b(i)=cmplx(x((i+n_poles)*2+1),x((i+n_poles)*2+2))!ATTENZIONE
enddo
!perform calculaation
do i=1,m
fvec(i)=0.d0
func=a_0
zz=cmplx(0.d0,freq(i))
do j=1,n_poles
func=func+a(j)/(zz-b(j))
enddo
func=func-c_target(i)
!fvec(i)=sqrt(real(func*conjg(func)))
fvec(i)=dble(func*conjg(func))
! do j=1,n_poles
! func=func+a(j)/(zz-b(j))
! enddo
! func=(func-c_target(i))/c_target(i)
! fvec(i)=real(func*conjg(func))
enddo
! write(*,*) 'fcn', fvec(1)
return
END SUBROUTINE fcn
SUBROUTINE fit_multipole_verlet2(n,m,z,s,a_0,a,b,thres,n_max_iterations, chi, dt, frice)
!fits with the function f(z)=a_0+\sum_{i=1,m} a_i/(z-b_i)
!the values z_j,s_j
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
implicit none
INTEGER, INTENT(in) :: n!numer of sampled values
INTEGER, INTENT(in) :: m!number of parameters a
COMPLEX(kind=DP), INTENT(in) :: z(n)!where
COMPLEX(kind=DP), INTENT(in) :: s(n)!values s(z_j) to be fitted
COMPLEX(kind=DP), INTENT(inout) :: a_0
COMPLEX(kind=DP), INTENT(inout) :: a(m)
COMPLEX(kind=DP), INTENT(inout) :: b(m)
REAL(kind=DP), INTENT(in) :: thres!threshold for convergence
INTEGER, INTENT(in) :: n_max_iterations!maximum number of search iterations
REAL(kind=DP), INTENT(out) :: chi!the final chi
REAL(kind=DP), INTENT(in) :: dt!time step
REAL(kind=DP), INTENT(in) :: frice!frice
REAL(kind=DP) :: chi0, chi1, dtt
INTEGER :: i,j,it
COMPLEX(kind=DP) :: cc
REAL(kind=DP), ALLOCATABLE :: omegas(:)
REAL(kind=DP), ALLOCATABLE :: x0(:),x(:),v(:),f(:),ma(:),x1(:)
INTEGER :: iflag
INTEGER :: np
np=2+4*m!number of parameters
allocate(omegas(n))
allocate(x0(np),x(np),v(np),f(np),ma(np),x1(np))
omegas(1:n)=aimag(z(1:n))
!calculates initial chi
chi0=0.d0
do i=1,n
cc = func(z(i))-s(i)
chi0=chi0+cc*conjg(cc)
enddo
write(stdout,*) 'Chi0 initial:', chi0
!set in variables
x(1)=real(a_0)
x(2)=aimag(a_0)
do i=1,m
x(i*2+1)=real(a(i))
x(i*2+2)=aimag(a(i))
x((i+m)*2+1)=real(b(i))
x((i+m)*2+2)=aimag(b(i))
enddo
!set up fcn function paramters
call fcn_set(n, m, omegas,s)
!intial values
ma(:)=1.d0
x0(:)=x(:)
call fcn_point(n,np,x,chi1,f)
write(stdout,*) 'VERLET2', chi1
x(:)=x0(:)+f(:)*dt/ma(:)
v(:)=0.d0
chi0=chi1
dtt=dt
do it=1,n_max_iterations
call fcn_point(n,np,x,chi1,f)
if(chi1 >= chi0) dtt=dtt/10.d0
chi0=chi1
if(mod(it,1000)==1) write(stdout,*) 'VERLET2', it, chi1
!x1(:)=2.d0*x(:)-x0(:)+f(:)*(dt**2.d0)/ma(:)-frice*v(:)*(dt**2.d0)
x1(:)=x(:)+f(:)*dtt/ma(:)
v(:)=(x1(:)-x0(:))/(2.d0*dt)
x0(:)=x(:)
x(:)=x1(:)
enddo
!set back parameters
a_0=cmplx(x(1),x(2))
do i=1,m
a(i)=cmplx(x(i*2+1),x(i*2+2))
enddo
do i=1,m
!b(i)=cmplx(x((i+m)*2+1),(-1.d0)*((x((i+m)*2+2))**2.d0))
b(i)=cmplx(x((i+m)*2+1),x((i+m)*2+2))
enddo
chi=chi1
write(stdout,*) 'FINAL CHI', chi!ATTENZIONE
deallocate(x,v,f,x0,ma,x1)
return
CONTAINS
FUNCTION func(zz)
COMPLEX(kind=DP) :: func
COMPLEX(kind=DP) :: zz
INTEGER :: ii
func=a_0
do ii=1,m
func=func+a(ii)/(zz-b(ii))
enddo
return
END FUNCTION func
END SUBROUTINE fit_multipole_verlet2
SUBROUTINE fit_multipole_minpack(n,m,z,s,a_0,a,b,thres,n_max_iterations, chi)
!fits with the function f(z)=a_0+\sum_{i=1,m} a_i/(z-b_i)
!the values z_j,s_j
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
implicit none
INTEGER, INTENT(in) :: n!numer of sampled values
INTEGER, INTENT(in) :: m!number of parameters a
COMPLEX(kind=DP), INTENT(in) :: z(n)!where
COMPLEX(kind=DP), INTENT(in) :: s(n)!values s(z_j) to be fitted
COMPLEX(kind=DP), INTENT(inout) :: a_0
COMPLEX(kind=DP), INTENT(inout) :: a(m)
COMPLEX(kind=DP), INTENT(inout) :: b(m)
REAL(kind=DP), INTENT(in) :: thres!threshold for convergence
INTEGER, INTENT(in) :: n_max_iterations!maximum number of search iterations
REAL(kind=DP), INTENT(out) :: chi!the final chi
REAL(kind=DP) :: chi0, chi1
INTEGER :: i,j,it
COMPLEX(kind=DP) :: cc
REAL(kind=DP), ALLOCATABLE :: omegas(:)
REAL(kind=DP), ALLOCATABLE :: variables(:)
INTEGER :: iflag
INTEGER :: info
INTEGER :: lwa
INTEGER :: np
INTEGER, ALLOCATABLE :: iwa(:)
INTEGER, ALLOCATABLE :: ipvt(:)
REAL(kind=DP), ALLOCATABLE :: wa(:)
REAL(kind=DP), ALLOCATABLE :: fvec(:),fjac(:,:)
EXTERNAL fcn,fcn2,fcnj
INTEGER :: ldfjac
ldfjac=n
allocate(omegas(n))
allocate(variables(2+4*m))
allocate(iwa(2+4*m))
lwa=n*(2+4*m)+5*(2+4*m)+n
allocate(wa(lwa))
allocate(fvec(n))
np=2+4*m
allocate(ipvt(np))
allocate(fjac(n,np))
write(stdout,*) 'Allocated'
FLUSH(stdout)
omegas(1:n)=aimag(z(1:n))
!calculates initial chi
chi0=0.d0
do i=1,n
cc = func(z(i))-s(i)
chi0=chi0+cc*conjg(cc)
enddo
! write(*,*) 'a_0', a_0!ATTENZIONE
! write(*,*) 'a', a
! write(*,*) 'b', b
! write(*,*) 'z,s' , z(1),s(1),func(z(1))
! write(*,*) 'z,s' , z(n),s(n),func(z(n))
write(stdout,*) 'Chi0 initial:', chi0
FLUSH(stdout)
!set in variables
variables(1)=real(a_0)
variables(2)=aimag(a_0)
do i=1,m
variables(i*2+1)=real(a(i))
variables(i*2+2)=aimag(a(i))
variables((i+m)*2+1)=real(b(i))
!variables((i+m)*2+2)=sqrt(abs(aimag(b(i))))
variables((i+m)*2+2)=aimag(b(i))!ATTENZIONE
enddo
!set up fcn function paramters
call fcn_set(n, m, omegas,s)
!call minpack driver routine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
info=1
call fcn(n,np,variables,fvec,info)
!write(*,*) 'FVEC', fvec(1:10)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!call lmdif1(fcn,n,np,n_max_iterations,variables,fvec,thres,info,iwa,wa,lwa)
call lmder1(fcnj,n,np,variables,fvec,fjac,ldfjac,thres,info,ipvt,wa,lwa,n_max_iterations)
write(stdout,*) ' INFO :', info,thres!ATTENZIONE
!set back parameters
a_0=dcmplx(variables(1),variables(2))
do i=1,m
a(i)=dcmplx(variables(i*2+1),variables(i*2+2))
enddo
do i=1,m
!b(i)=dcmplx(variables((i+m)*2+1),(-1.d0)*((variables((i+m)*2+2))**2.d0))
b(i)=dcmplx(variables((i+m)*2+1),variables((i+m)*2+2))
enddo
!recalculate chi and write result on stdout
!calculates initial chi
chi0=0.d0
do i=1,n
cc = (func(z(i))-s(i))
chi0=chi0+cc*conjg(cc)
enddo
!write(stdout,*) 'a_0', a_0!ATTENZIONE
!write(stdout,*) 'a', a
!write(stdout,*) 'b', b
!write(stdout,*) 'z,s' , z(1),s(1),func(z(1))
!write(stdout,*) 'z,s' , z(n),s(n),func(z(n))
write(stdout,*) 'Minpack fit chi0 :', chi0
chi=chi0
deallocate(omegas)
deallocate(variables)
deallocate(iwa)
deallocate(wa)
deallocate(fvec)
deallocate(ipvt)
deallocate(fjac)
return
CONTAINS
FUNCTION func(zz)
COMPLEX(kind=DP) :: func
COMPLEX(kind=DP) :: zz
INTEGER :: ii
func=a_0
do ii=1,m
func=func+a(ii)/(zz-b(ii))
enddo
return
END FUNCTION func
END SUBROUTINE fit_multipole_minpack
SUBROUTINE fcn2(m,n,x,fvec,iflag)
!added internal parmeters to set up function
!the parameters are in the order re(a_0),im(a_0),re(a_1)...re(a_n),im(a_1)..im(a_n),re(b_1)..re(b_n), im(b_1)..im(b_n)
use kinds, ONLY : DP
use io_global, ONLY : stdout
implicit none
INTEGER :: m !number of variables
INTEGER :: n!total number of parameters
REAL(kind=DP) :: x(n)!parameters
REAL(kind=DP) :: fvec(m)!evaluated error function
INTEGER :: iflag!not used
fvec(1:m)=0.d0
END SUBROUTINE
SUBROUTINE fcnj(m,n,x,fvec,fjac,ldfjac,iflag)
!this version calculates also the jacobian
!the parameters are in the order re(a_0),im(a_0),re(a_1)...re(a_n),im(a_1)..im(a_n),re(b_1)..re(b_n), im(b_1)..im(b_n)
use kinds, ONLY : DP
use io_global, ONLY : stdout
use global_minpack
implicit none
INTEGER :: m !number of variables
INTEGER :: n!total number of parameters
REAL(kind=DP) :: x(n)!parameters
REAL(kind=DP) :: fvec(m)!evaluated error function
REAL(kind=DP) :: fjac(ldfjac,n)
INTEGER :: ldfjac!leading dimension of fjac
INTEGER :: iflag! =1 calculate fvec, =2 calculate fjac
COMPLEX(kind=DP) a_0, a(maxpole), b(maxpole),g,h
INTEGER i,j
COMPLEX(kind=DP) :: func,zz
if(m>maxm) then
write(stdout,*) 'FCN: MAXN TOO SMALL'
stop
endif
!set up parameters
a_0=dcmplx(x(1),x(2))
do i=1,n_poles
a(i)=dcmplx(x(i*2+1),x(i*2+2))
enddo
do i=1,n_poles
!b(i)=dcmplx(x((i+n_poles)*2+1),-(x((i+n_poles)*2+2))**2.d0)
b(i)=dcmplx(x((i+n_poles)*2+1),x((i+n_poles)*2+2))
enddo
if(iflag==1) then
!perform calculaation
do i=1,m
fvec(i)=0.d0
func=a_0
zz=dcmplx(0.d0,freq(i))
do j=1,n_poles
func=func+a(j)/(zz-b(j))
enddo
func=func-c_target(i)
fvec(i)=dble(func*conjg(func))
enddo
else if(iflag==2) then
do j=1,m
fjac(j,:)=0.d0
!calculate g_j
g=a_0
zz=cmplx(0.d0,freq(j))
do i=1,n_poles
g=g+a(i)/(zz-b(i))
enddo
g=g-c_target(j)
!now term a_0
fjac(j,1)=2.d0*real(g)
fjac(j,2)=2.d0*aimag(g)
!now terms a_i
do i=1, n_poles
h=(1.d0,0.d0)/(zz-b(i))
fjac(j,i*2+1)=2.d0*real(h*conjg(g))
fjac(j,i*2+2)=-2.d0*aimag(h*conjg(g))
enddo
!now terms b_i
do i=1, n_poles
h=a(i)/((zz-b(i))**2.d0)
fjac(j,(i+n_poles)*2+1)=2.d0*real(h*conjg(g))
!fjac(j,(i+n_poles)*2+2)=-2.d0*aimag(h*conjg(g))*(-2.d0*x((i+n_poles)*2+2))
fjac(j,(i+n_poles)*2+2)=-2.d0*aimag(h*conjg(g))
enddo
enddo
endif
return
END SUBROUTINE fcnj
SUBROUTINE fcn_point(m,n,x,value,nabla)
!this version calculates also the jacobian
!the parameters are in the order re(a_0),im(a_0),re(a_1)...re(a_n),im(a_1)..im(a_n),re(b_1)..re(b_n), im(b_1)..im(b_n)
use kinds, ONLY : DP
use io_global, ONLY : stdout
use global_minpack
implicit none
INTEGER :: m !number of variables
INTEGER :: n!total number of parameters
REAL(kind=DP) :: x(n)!parameters
REAL(kind=DP) :: value!total error
REAL(kind=DP) :: nabla(n)!derivatives of total error respect to parameters
COMPLEX(kind=DP) a_0, a(maxpole), b(maxpole),g,h
INTEGER i,j
COMPLEX(kind=DP) :: func,zz
if(m>maxm) then
write(stdout,*) 'FCN: MAXN TOO SMALL'
stop
endif
!set up parameters
a_0=cmplx(x(1),x(2))
do i=1,n_poles
a(i)=cmplx(x(i*2+1),x(i*2+2))
enddo
do i=1,n_poles
! b(i)=cmplx(x((i+n_poles)*2+1),-(x((i+n_poles)*2+2))**2.d0)
b(i)=cmplx(x((i+n_poles)*2+1),x((i+n_poles)*2+2))
enddo
!perform calculation of value
value=0.d0
do i=1,m
func=a_0
zz=cmplx(0.d0,freq(i))
do j=1,n_poles
func=func+a(j)/(zz-b(j))
enddo
func=func-c_target(i)
value=value + dble(func*conjg(func))
enddo
nabla(:)=0.d0
do j=1,m
!calculate g_j
g=a_0
zz=cmplx(0.d0,freq(j))
do i=1,n_poles
g=g+a(i)/(zz-b(i))
enddo
g=g-c_target(j)
!now term a_0
nabla(1)=nabla(1)+2.d0*real(g)
nabla(2)=nabla(2)+2.d0*aimag(g)
!now terms a_i
do i=1, n_poles
h=(1.d0,0.d0)/(zz-b(i))
nabla(i*2+1)=nabla(i*2+1)+2.d0*real(h*conjg(g))
nabla(i*2+2)=nabla(i*2+1)-2.d0*aimag(h*conjg(g))
enddo
!now terms b_i
do i=1, n_poles
h=a(i)/((zz-b(i))**2.d0)
nabla((i+n_poles)*2+1)=nabla((i+n_poles)*2+1)+2.d0*real(h*conjg(g))
nabla((i+n_poles)*2+2)=nabla((i+n_poles)*2+2)-2.d0*aimag(h*conjg(g))
enddo
enddo
nabla(:)=-nabla(:)
return
END SUBROUTINE fcn_point
Reference in New Issue View Git Blame Copy Permalink