doxygen/v8.037/mnimpr_8f_source.html

 *
 * $Id: mnimpr.F,v 1.1.1.1 1996/03/07 14:31:30 mclareni Exp $
 *
 * $Log: mnimpr.F,v $
 * Revision 1.1.1.1  1996/03/07 14:31:30  mclareni
 * Minuit
 *
 *
       SUBROUTINE mnimpr(FCN,FUTIL)
 *
 * $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
 *
 * $Log: d506dp.inc,v $
 * Revision 1.1.1.1  1996/03/07 14:31:32  mclareni
 * Minuit
 *
 *
 *
 *
 * d506dp.inc
 *
 C ************ DOUBLE PRECISION VERSION *************
       IMPLICIT DOUBLE PRECISION (a-h,o-z)
 CC        Attempts to improve on a good local minimum by finding a
 CC        better one.   The quadratic part of FCN is removed by MNCALF
 CC        and this transformed function is minimized using the simplex
 CC        method from several random starting points.
 CC        ref. -- Goldstein and Price, Math.Comp. 25, 569 (1971)
 CC
 *
 * $Id: d506cm.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
 *
 * $Log: d506cm.inc,v $
 * Revision 1.1.1.1  1996/03/07 14:31:32  mclareni
 * Minuit
 *
 *
 *
 *
 * d506cm.inc
 *
       parameter(mne=100 , mni=50)
       parameter(mnihl=mni*(mni+1)/2)
       CHARACTER*10 CPNAM
       COMMON
      1/mn7nam/ cpnam(mne)
      2/mn7ext/ u(mne)     ,alim(mne)  ,blim(mne)
      3/mn7err/ erp(mni)   ,ern(mni)   ,werr(mni)  ,globcc(mni)
      4/mn7inx/ nvarl(mne) ,niofex(mne),nexofi(mni)
      5/mn7int/ x(mni)     ,xt(mni)    ,dirin(mni)
      6/mn7fx2/ xs(mni)    ,xts(mni)   ,dirins(mni)
      7/mn7der/ grd(mni)   ,g2(mni)    ,gstep(mni) ,gin(mne) ,dgrd(mni)
      8/mn7fx3/ grds(mni)  ,g2s(mni)   ,gsteps(mni)
      9/mn7fx1/ ipfix(mni) ,npfix
      a/mn7var/ vhmat(mnihl)
      b/mn7vat/ vthmat(mnihl)
      c/mn7sim/ p(mni,mni+1),pstar(mni),pstst(mni) ,pbar(mni),prho(mni)
 C
       parameter(maxdbg=10, maxstk=10, maxcwd=20, maxp=30, maxcpt=101)
       parameter(zero=0.0,  one=1.0,   half=0.5)
       COMMON
      d/mn7npr/ maxint ,npar   ,maxext ,nu
      e/mn7iou/ isysrd ,isyswr ,isyssa ,npagwd ,npagln ,newpag
      e/mn7io2/ istkrd(maxstk) ,nstkrd ,istkwr(maxstk) ,nstkwr
      f/mn7tit/ cfrom  ,cstatu ,ctitl  ,cword  ,cundef ,cvrsn ,covmes
      g/mn7flg/ isw(7) ,idbg(0:maxdbg) ,nblock ,icomnd
      h/mn7min/ amin   ,up     ,edm    ,fval3  ,epsi   ,apsi  ,dcovar
      i/mn7cnv/ nfcn   ,nfcnmx ,nfcnlc ,nfcnfr ,itaur,istrat,nwrmes(2)
      j/mn7arg/ word7(maxp)
      k/mn7log/ lwarn  ,lrepor ,limset ,lnolim ,lnewmn ,lphead
      l/mn7cns/ epsmac ,epsma2 ,vlimlo ,vlimhi ,undefi ,bigedm,updflt
      m/mn7rpt/ xpt(maxcpt)    ,ypt(maxcpt)
      n/mn7cpt/ chpt(maxcpt)
      o/mn7xcr/ xmidcr ,ymidcr ,xdircr ,ydircr ,ke1cr  ,ke2cr
       CHARACTER CTITL*50, CWORD*(maxcwd), CUNDEF*10, CFROM*8,
      +          cvrsn*6,  covmes(0:3)*22, cstatu*10, chpt*1
       LOGICAL   LWARN, LREPOR, LIMSET, LNOLIM, LNEWMN, LPHEAD
       EXTERNAL fcn,futil
       dimension dsav(mni), y(mni+1)
       parameter(alpha=1.,beta=0.5,gamma=2.0)
       DATA rnum/0./
       IF (npar .LE. 0)  RETURN
       IF (amin .EQ. undefi)  CALL mnamin(fcn,futil)
       cstatu = 'UNCHANGED '
       itaur = 1
       epsi = 0.1*up
       npfn=nfcn
       nloop = word7(2)
       IF (nloop .LE. 0)  nloop = npar + 4
       nparx = npar
       nparp1=npar+1
       wg = 1.0/float(npar)
       sigsav = edm
       apsi = amin
          DO 2 i= 1, npar
          xt(i) = x(i)
          dsav(i) = werr(i)
            DO 2 j = 1, i
            ndex = i*(i-1)/2 + j
            p(i,j) = vhmat(ndex)
     2      p(j,i) = p(i,j)
       CALL mnvert(p,maxint,maxint,npar,ifail)
       IF (ifail .GE. 1)  GO TO 280
 C               Save inverted matrix in VT
          DO 12 i= 1, npar
          ndex = i*(i-1)/2
            DO 12 j= 1, i
            ndex = ndex + 1
    12      vthmat(ndex) = p(i,j)
       loop = 0
 C
    20 CONTINUE
          DO 25 i= 1, npar
          dirin(i) = 2.0*dsav(i)
          CALL mnrn15(rnum,iseed)
    25    x(i) = xt(i) + 2.0*dirin(i)*(rnum-0.5)
       loop = loop + 1
       reg = 2.0
       IF (isw(5) .GE. 0)   WRITE (isyswr, 1040) loop
    30 CALL  mncalf(fcn,x,ycalf,futil)
       amin = ycalf
 C                                        . . . . set up  random simplex
       jl = nparp1
       jh = nparp1
       y(nparp1) = amin
       amax = amin
          DO 45 i= 1, npar
          xi = x(i)
          CALL mnrn15(rnum,iseed)
          x(i) = xi - dirin(i) *(rnum-0.5)
          CALL mncalf(fcn,x,ycalf,futil)
          y(i) = ycalf
          IF (y(i) .LT. amin)  THEN
             amin = y(i)
             jl = i
          ELSE IF (y(i) .GT. amax)  THEN
             amax = y(i)
             jh = i
          ENDIF
             DO 40 j= 1, npar
    40       p(j,i) = x(j)
          p(i,nparp1) = xi
          x(i) = xi
    45    CONTINUE
 C
       edm = amin
       sig2 = edm
 C                                        . . . . . . .  start main loop
    50 CONTINUE
       IF (amin .LT. zero)  GO TO 95
       IF (isw(2) .LE. 2)  GO TO 280
       ep = 0.1*amin
       IF (sig2 .LT. ep   .AND. edm.LT.ep  )     GO TO 100
       sig2 = edm
       IF ((nfcn-npfn) .GT. nfcnmx)  GO TO 300
 C         calculate new point * by reflection
       DO 60 i= 1, npar
       pb = 0.
       DO 59 j= 1, nparp1
    59 pb = pb + wg * p(i,j)
       pbar(i) = pb - wg * p(i,jh)
    60 pstar(i)=(1.+alpha)*pbar(i)-alpha*p(i,jh)
       CALL mncalf(fcn,pstar,ycalf,futil)
       ystar = ycalf
       IF(ystar.GE.amin) GO TO 70
 C         point * better than jl, calculate new point **
       DO 61 i=1,npar
    61 pstst(i)=gamma*pstar(i)+(1.-gamma)*pbar(i)
       CALL mncalf(fcn,pstst,ycalf,futil)
       ystst = ycalf
    66 IF (ystst .LT. y(jl))  GO TO 67
       CALL mnrazz(ystar,pstar,y,jh,jl)
       GO TO 50
    67 CALL mnrazz(ystst,pstst,y,jh,jl)
       GO TO 50
 C         point * is not as good as jl
    70 IF (ystar .GE. y(jh))  GO TO 73
       jhold = jh
       CALL mnrazz(ystar,pstar,y,jh,jl)
       IF (jhold .NE. jh)  GO TO 50
 C         calculate new point **
    73 DO 74 i=1,npar
    74 pstst(i)=beta*p(i,jh)+(1.-beta)*pbar(i)
       CALL mncalf(fcn,pstst,ycalf,futil)
       ystst = ycalf
       IF(ystst.GT.y(jh)) GO TO 30
 C     point ** is better than jh
       IF (ystst .LT. amin)  GO TO 67
       CALL mnrazz(ystst,pstst,y,jh,jl)
       GO TO 50
 C                                        . . . . . .  end main loop
    95 IF (isw(5) .GE. 0)  WRITE (isyswr,1000)
       reg = 0.1
 C                                        . . . . . ask if point is new
   100 CALL mninex(x)
       CALL fcn(nparx,gin,amin,u,4,futil)
       nfcn = nfcn + 1
       DO 120 i= 1, npar
       dirin(i) = reg*dsav(i)
       IF (abs(x(i)-xt(i)) .GT. dirin(i)) GO TO 150
   120 CONTINUE
       GO TO 230
   150 nfcnmx = nfcnmx + npfn - nfcn
       npfn = nfcn
       CALL mnsimp(fcn,futil)
       IF (amin .GE. apsi)  GO TO 325
       DO 220 i= 1, npar
       dirin(i) = 0.1 *dsav(i)
       IF (abs(x(i)-xt(i)) .GT. dirin(i)) GO TO 250
   220 CONTINUE
   230 IF (amin .LT. apsi)  GO TO 350
       GO TO 325
 C                                        . . . . . . truly new minimum
   250 lnewmn = .true.
       IF (isw(2) .GE. 1) THEN
           isw(2) = 1
           dcovar = max(dcovar,half)
       ELSE
           dcovar = 1.
       ENDIF
       itaur = 0
       nfcnmx = nfcnmx + npfn - nfcn
       cstatu = 'NEW MINIMU'
       IF (isw(5) .GE. 0)      WRITE (isyswr,1030)
       RETURN
 C                                        . . . return to previous region
   280 IF (isw(5) .GT. 0) WRITE (isyswr,1020)
       GO TO 325
   300 isw(1) = 1
   325 DO 330 i= 1, npar
       dirin(i) = 0.01*dsav(i)
   330 x(i) = xt(i)
       amin = apsi
       edm = sigsav
   350 CALL mninex(x)
       IF (isw(5) .GT. 0)    WRITE (isyswr,1010)
       cstatu= 'UNCHANGED '
       CALL mnrset(0)
       IF (isw(2) .LT. 2)  GO TO 380
       IF (loop .LT. nloop .AND. isw(1) .LT. 1)  GO TO 20
   380 CALL mnprin (5,amin)
       itaur = 0
       RETURN
  1000 FORMAT (54h an improvement on the previous minimum has been found)
  1010 FORMAT (51h improve has returned to region of original minimum)
  1020 FORMAT (/44h covariance matrix was not positive-definite)
  1030 FORMAT (/38h improve has found a truly new minimum/1h ,37(1h*)/)
  1040 FORMAT (/18h start attempt no.,i2,  20h to find new minimum)
       END
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst

z
nodes z
Definition: cstdatmos0.fNew2.h:93

up
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data h g *is for param c g data up(2, 1)/7.0d0/

e
dE dx *! Nuc Int sampling table e
Definition: cblkMuInt.h:130

i
nodes i
Definition: cstdatmos0.fNew2.h:93

mnrset
subroutine mnrset(IOPT)
Definition: mnrset.f:10

true
block data cblkElemag data *AnihiE ! Eposi< 1 TeV, anihilation considered *X0/365.667/, ! radiation length of air in kg/m2 *Ecrit/81.e-3/, ! critical energy of air in GeV *MaxComptonE/1./, ! compton is considered below 1 GeV *MaxPhotoE/1.e-3/, ! above this, PhotoElectric effect neg. *MinPhotoProdE/153.e-3/, ! below 153 MeV, no gp --> hadrons ! scattering const not MeV *Knockon true
Definition: cblkElemag.h:7

mnvert
subroutine mnvert(A, L, M, N, IFAIL)
Definition: mnvert.f:25

mnprin
subroutine mnprin(INKODE, FVAL)
Definition: mnprin.f:10

mnimpr
subroutine mnimpr(FCN, FUTIL)
Definition: mnimpr.f:10

xs
! constants thru Cosmos real ! if multiplied to deg radian Torad ! light velocity m sec ! infinty ! kg m2 *Togpcm2 g cm2 ! g cm2 *Tokgpm2 kg m2 ! cm *Tom m ! m *Tocm cm ! g cm3 *Tokgpm3 kg m3 ! kg m3 *Togpcm3 g cm3 ! sec *Tonsec nsec ! Tesla m ! Avogadro *A2deninv ! mfp *n * xs
Definition: Zglobalc.h:18

mnamin
subroutine mnamin(FCN, FUTIL)
Definition: mnamin.f:10

p
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7

j
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36

to
integer maxbin nregion c minsave drx2 ! drawing region real maxdep integer maxpos integer to
Definition: Zfit.h:15

mncalf
subroutine mncalf(FCN, PVEC, YCALF, FUTIL)
Definition: mncalf.f:10

o
struct ob o[NpMax]
Definition: Zprivate.h:34

mninex
subroutine mninex(PINT)
Definition: mninex.f:10

y
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25

d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130

mnsimp
subroutine mnsimp(FCN, FUTIL)
Definition: mnsimp.f:25

a
nodes a
Definition: cstdatmos0.fNew2.h:93

b
dE dx *! Nuc Int sampling table b
Definition: cblkMuInt.h:130

m
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data m
Definition: ZlatfitBD.h:35

h
dE dx *! Nuc Int sampling table h
Definition: cblkMuInt.h:130

g
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130

mnrazz
subroutine mnrazz(YNEW, PNEW, Y, JH, JL)
Definition: mnrazz.f:10

maxp
!onst int maxp
Definition: Zprivate.h:3

mnrn15
subroutine mnrn15(VAL, INSEED)
Definition: mnrn15.f:10

n
integer n
Definition: Zcinippxc.h:1

softenpik::half
integer, parameter half
Definition: csoftenPiK.f:108

x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21

f
dE dx *! Nuc Int sampling table f
Definition: cblkMuInt.h:130

c
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130