[Fortran] 纯文本查看 复制代码
program polesearching
implicit none
!!! two variables
integer Npar,maxiter,inter
parameter(Npar=2)
real*8 par(Npar),del(Npar),error,fvalue,yourf2v
!! one variable to search zero
real*8 yourf1v,x,lower,upper,x0,rr,repsdzero
integer mxfdzero
external yourf1v
!! two variable for subroutine: nelmin
real*8 xmin(Npar),step(Npar),reqmin
integer konvge,kcount,icount,numres,ifault
!! two variable for subroutine: MINIMIZE_2D
real*8 par0(Npar)
integer info
external subyourf2v
maxiter=1500
inter=0
error= 1.d-8
!!
par(1)= 1011.0d0
par(2)= 0.0d0
del(1)= 0.001d0
del(2)= 0.001d0
print*, yourf2v(par)
call minsimplex2(Npar,par,del,error,maxiter,fvalue,yourf2v,inter)
write(*,'(A10,2F20.10,E12.2)') 'solution', par(1),par(2),fvalue
end program polesearching
function yourf2v(par)
implicit none
real*8 yourf2v
integer Npar
parameter(Npar=2)
real*8 par(Npar),x,y
complex*16 ui,s,E
real*8 mp,mq,mpi,mk,meta,metap,Fpi,mu
real*8 ct,st
real*8 asL10n1, asL10n2, asL10n3
integer nchan
parameter (nchan=3)
complex*16 gij10(nchan,nchan),gij102(nchan,nchan),
-tijgij10(nchan,nchan),ro10(nchan,nchan),
-i33tijgij10(nchan,nchan),ss(nchan,nchan),
-pwTij10TollerM33MLo(nchan,nchan),s01(nchan,nchan),s11(nchan,nchan)
-, EtauNijMlo10(nchan,nchan),DeltaLoguNijMlo10(nchan,nchan),
-EtauNijMlo1d2(nchan,nchan),DeltaLoguNijMlo1d2(nchan,nchan)
integer info !qiu inverse matrix
integer ipvt(nchan) !qiu inverse matrix
complex*16 work(nchan),detin(2),deterin !qiu inverse matrix
complex*16 pwTij10p2M33MAnl(nchan,nchan)
complex*16 i33umat(nchan,nchan)
real*8 pi
common /pi/ pi
real*8 asl10
common /lec10/ asl10
common /physicalmass/mpi,mk,meta,metap,fpi,mu
complex*16 GfdrI
external GfdrI
pi= dacos(-1.0d0)
mpi=137.3d0
mk=495.6d0
meta=547.9d0
metap=957.7d0
Fpi =92.4d0
mu=770.0d0
asL10 = -1.44d0
asL10n1 = asl10
asL10n2 = asL10
asL10n3 = asL10
st=-0.279d !theta is -16.2
ct=0.96d0 !theta is -16.2
ui=(0.d0,1.d0)
E=par(1)+ui*par(2)
s=E**2
gij10(1,1)=GfdrI(asL10n1,s,mpi**2,meta**2)
gij10(1,2)=(0.0d0,0.0d0)
gij10(1,3)=(0.0d0,0.0d0)
gij10(2,1)=(0.0d0,0.0d0)
gij10(2,2)=GfdrI(asL10n2,s,mk**2,mk**2)
gij10(2,3)=(0.0d0,0.0d0)
gij10(3,1)=(0.0d0,0.0d0)
gij10(3,2)=(0.0d0,0.0d0)
gij10(3,3)=GfdrI(asL10n3,s,mpi**2,metaP**2)
pwTij10p2M33MAnl(1,1)=(mpi**2*(ct**2 - 2.0*Sqrt(2.0)*ct*st
- + 2.0*st**2))/(3.0*Fpi**2)
pwTij10p2M33MAnl(1,2)= (ct*(3.0*meta**2 + 8.0*mk**2 + mpi**2
- - 9.0*s) + 2.0*Sqrt(2.0)*(2.0*mk**2 + mpi**2)*st)/
-(6.0*Sqrt(6.0)*Fpi**2)
pwTij10p2M33MAnl(1,3)=(mpi**2*(Sqrt(2.0)*ct**2 -
-ct*st - Sqrt(2.0)*st**2))/(3.0*Fpi**2)
pwTij10p2M33MAnl(2,1)=(ct*(3.0*meta**2 +
-8.0*mk**2 + mpi**2 - 9.0*s) +
- 2.0*Sqrt(2.0)*(2.0*mk**2 + mpi**2)*st)/(6.0*Sqrt(6.0)*Fpi**2)
pwTij10p2M33MAnl(2,2)=s/(4.0*Fpi**2)
pwTij10p2M33MAnl(2,3)= -(2.0*Sqrt(2.0)*ct*(2.0*mk**2 + mpi**2) -
- (3.0*metaP**2 + 8.0*mk**2 + mpi**2 - 9.0*s)*st)/
- (6.0*Sqrt(6.0)*Fpi**2)
pwTij10p2M33MAnl(3,1)=(mpi**2*(Sqrt(2.0)*ct**2 -
-ct*st - Sqrt(2.0)*st**2))/(3.0*Fpi**2)
pwTij10p2M33MAnl(3,2)= -(2.0*Sqrt(2.0)*ct*(2.0*mk**2 + mpi**2) -
- (3.0*metaP**2 + 8.0*mk**2 + mpi**2 - 9.0*s)*st)/
- (6.0*Sqrt(6.0)*Fpi**2)
pwTij10p2M33MAnl(3,3)=(mpi**2*(2.0*ct**2 + 2.0*Sqrt(2.0)*ct*st +
-st**2))/(3.0*Fpi**2)
i33umat(1,1)=(1.0d0,0.0d0)
i33umat(1,2)=(0.0d0,0.0d0)
i33umat(1,3)=(0.0d0,0.0d0)
i33umat(2,1)=(0.0d0,0.0d0)
i33umat(2,2)=(1.0d0,0.0d0)
i33umat(2,3)=(0.0d0,0.0d0)
i33umat(3,1)=(0.0d0,0.0d0)
i33umat(3,2)=(0.0d0,0.0d0)
i33umat(3,3)=(1.0d0,0.0d0)
tijgij10=MATMUL(pwTij10p2M33MAnl,gij10)
i33tijgij10=i33umat- tijgij10
call zgefa( i33tijgij10,nchan,nchan,ipvt,info) !! call both
call zgedi(i33tijgij10,nchan,nchan,ipvt,detin,work,11) !! call both
deterin=detin(1)*10.**detin(2) !! determinant
print*, deterin
pwTij10TollerM33MLo=MATMUL(i33tijgij10,pwTij10p2M33MAnl) !no problem
yourf2v=deterin
end
SUBROUTINE minsimplex2
1 (NVAR,PUNT,DELT,ERROR,MAXITER,VALOR,FUNCIO,INTER)
IMPLICIT REAL*8 (a-h,o-z)
PARAMETER(nmax=30)
DIMENSION PUNT(nvar),DELT(nvar)
DIMENSION Y(nmax+1), P(nmax,nmax+1)
common /iierror/ierror !(mio)
EXTERNAL FUNCIO
DO j=1,nvar+1
DO i=1, nvar
p(i,j)=PUNT(i)
ENDDO
IF(j .GT. 1) THEN
p(j-1,j)=punt(j-1)+delt(j-1)
ENDIF
ENDDO
10 CONTINUE
CALL SIMPLEX2
1 (P,Y,nvar,nmax,error,ierror,maxiter,funcio,millor,ipitjor)
IF (INTER .EQ. 1) THEN
WRITE(*,*)
WRITE(*,*)' ============================================='
WRITE(*,'(1X,3A20)') 'Par. Inic.','Par. Obtingu.', 'ERROR'
DO i=1,NVAR
WRITE(*,'(1X, 3F20.8)') PUNT(i),P(i,millor),
1 ABS(P(i,millor)-P(i,ipitjor))/2
ENDDO
WRITE(*,'(1X,A13,1P,E20.8,5x,a13,E20.8)')
1 ' FUNC. MILLOR', y(millor), ' FUNC. PITJOR',y(ipitjor)
IF (ierror .EQ. 0) THEN
WRITE(*,*) ' SIMPLEX arrib?a converg妌cia'
WRITE(*,*) ' Entreu 1 -> dividir error per 10, 0 -> finir'
READ(*,*) icas
IF (icas .EQ. 1) THEN
error=error/10
WRITE(*,*) ' Nou error=',error
GOTO 10
ENDIF
ELSE
WRITE(*,*) ' No s''ha aconseguit la converg妌cia'
WRITE(*,*) ' Entreu 1 per ITERAR, 0 per FINIR'
READ(*,*) icas
IF (icas .EQ. 1) GOTO 10
ENDIF
ENDIF
c Actualitzar DELTA i PUNT
valor=Y(millor)
DO i=1, nvar
PUNT(i)=P(i,millor)
vmin=P(i,1)
c print *,'kkk'
vmax=P(i,1)
DO j=2,nvar+1
vmin=MIN(vmin,p(i,j))
vmax=MAX(vmax,p(i,j))
ENDDO
DELT(I)=vmax-vmin
ENDDO
RETURN
END
SUBROUTINE simplex2
1 (P,Y,NVAR,NDIM,ERROR,IERROR,MAXITER,FMV,MILLOR,PITJOR)
IMPLICIT REAL*8 (a-h, o-z)
PARAMETER(nmax=30, unitat=1.d0)
DIMENSION p(ndim,ndim+1), y(ndim+1), iordre(nmax+1)
DIMENSION cdg(nmax+1), pn(nmax+1), pp(nmax+1)
INTEGER millor, pitjor , quasipitjor
common /espolo/iespolo !(mio)
IF (nvar .GT. nmax) THEN
WRITE(*,*) ' Mass variables in SIMPLEX'
STOP
ENDIF
IERROR=0
nvert=nvar+1
DO 1000 iter=1, maxiter
CALL ordenar2(y,iordre,nvert)
millor=iordre(1)
pitjor=iordre(nvert)
quasipitjor=iordre(nvert-1)
den = ABS(y(pitjor)+y(millor))/2.d0
c(mio) IF ((y(pitjor)-y(millor)) .LT. error*den) RETURN
aakk=abs(y(pitjor)/y(millor))-1.d0 !(mio)
errorkk=1.d-6 !(mio)
IF (aakk.LT.errorkk) then
iespolo=0
return
endif
c print *,'a',iter,y(millor),y(pitjor),y(pitjor)/y(millor)
c + ,(y(pitjor)-y(millor))/den
if (y(millor).LT.error) then
iespolo=1
return
endif
IF (iter .EQ. maxiter) GOTO 1000
c Moviments ... Primer busque CDG
DO i=1,nvar
cdg(i)=0
DO j=1, nvert
IF (j .NE. pitjor) cdg(i)=cdg(i)+P(i,j)
ENDDO
cdg(i)=cdg(i)/nvar
ENDDO
CALL combinacio2(-1.d0, P(1,pitjor), 2.d0, CDG, pn, nvar)
F1=fmv(pn)
F2=fmv(pp)
IF (f1 .GT. f2) THEN
y(pitjor)=f2
DO i=1,nvar
p(i,pitjor)=pp(i)
ENDDO
GOTO 1000
ENDIF
GOTO 2000
40 if (F1 .lt .Y(QUASIPITJOR)) goto 2000
IF (f1 .GT. y(pitjor)) GOTO 3000
4000 y(pitjor)=f1
DO i=1,nvar
p(i,pitjor)=pn(i)
ENDDO
ccc (m) hay problemas si pitjor es igual a cdg:
akk1=abs(P(1,pitjor)/cdg(1)-1.d0)
akk2=abs(P(2,pitjor)/cdg(2)-1.d0)
if ((akk1.LT.1.d-6).and.(akk2.LT.1.d-6)) then
iespolo=0
return
endif
c Contracci?del pitjor punt (interior o exterior)
3000 CALL combinacio2(0.5d0, P(1,pitjor), 0.5d0, CDG, pn, nvar)
F1=fmv(pn)
c Acceptable?
IF (f1 .LT. y(quasipitjor)) GOTO 2000
c No, per?si 俿 MILLOR que PITJOR l'accepte i contrac de nou
IF (F1 .LT. y(pitjor)) GOTO 4000
c CAS 4: Contracci?de tot el splex
DO j=1,nvert
IF (j .NE. millor) THEN
CALL combinacio2(0.5d0,P(1,j),0.5d0,P(1,millor),P(1,j),nvar)
Y(j)=fmv(pn)
ENDIF
ENDDO
goto 1000
c Acceptar pn i comen嘺r de nou
2000 CONTINUE
y(pitjor)=f1
DO i=1,nvar
p(i,pitjor)=pn(i)
ENDDO
1000 CONTINUE
c S'arriba ac?si no convergeix en MAXITER
c!(mio) WRITE(*,*) ' Massa iteracions en SIMPLEX'
ierror=1
RETURN
END
!!!!!!!!!!!!!!!!!
SUBROUTINE combinacio2(a, P, b, Q, R, n)
IMPLICIT REAL*8 (a-h, o-z)
DIMENSION p(n), q(n), r(n)
DO 10 i=1, n
10 r(i)=a*p(i)+b*q(i)
RETURN
END
!!!!!!!!!!!!!!!
SUBROUTINE ordenar2(f,iordre,n)
IMPLICIT REAL*8 (a-h,o-z)
DIMENSION F(n), IORDRE(n)
iordre(1)=1
DO 10 j=2,n
x=f(j)
DO 20 i=j-1,1,-1
IF (x .GT. f(iordre(i)) ) GOTO 10
20 iordre(i+1)=iordre(i)
i=0
10 iordre(i+1)=j
RETURN
END
F 币
贡献
权杖
窥视卡
发表于 2019-2-27 17:25:38
微信
收藏
点赞
点踩
提升卡
置顶卡
变色卡
千斤顶
显身卡
发表于 2019-2-27 18:18:55
楼主