[求助] 寻找极点

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 splex
        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