mncont.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	mncont (FCN, KE1, KE2, NPTU, XPTU, YPTU, IERRF, FUTIL)

Function/Subroutine Documentation

mncont()

subroutine mncont	(	external	FCN,
			KE1,
			KE2,
			NPTU,
		dimension(nptu)	XPTU,
		dimension(nptu)	YPTU,
			IERRF,
		external	FUTIL
	)

Definition at line 10 of file mncont.f.

References a, b, c, d, e, f, g, h, softenpik::half, i, j, m, maxp, mncros(), mncuve(), mnfixp(), mnfree(), mninex(), mnmnot(), mnplot(), mnwarn(), n, o, p, parameter(), up(), x, xs, and z.

Referenced by mnexcm().

*
* $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       Find NPTU points along a contour where the function
CC             FMIN (X(KE1),X(KE2)) =  AMIN+UP
CC       where FMIN is the minimum of FCN with respect to all
CC       the other NPAR-2 variable parameters (if any).
CC   IERRF on return will be equal to the number of points found:
CC     NPTU if normal termination with NPTU points found
CC     -1   if errors in the calling sequence (KE1, KE2 not variable)
CC      0   if less than four points can be found (using MNMNOT)
CC     n>3  if only n points can be found (n < NPTU)
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
      dimension xptu(nptu), yptu(nptu), w(mni),gcc(mni)
      CHARACTER chere*10
      parameter(chere='MNContour ')
      LOGICAL ldebug
      EXTERNAL fcn,futil
C                 input arguments: parx, pary, devs, ngrid
      ldebug = (idbg(6) .GE. 1)
      IF (ke1.LE.0 .OR. ke2.LE.0)  GO TO 1350
      IF (ke1.GT.nu .OR. ke2.GT.nu)  GO TO 1350
      ki1 = niofex(ke1)
      ki2 = niofex(ke2)
      IF (ki1.LE.0 .OR. ki2.LE.0)  GO TO 1350
      IF (ki1 .EQ. ki2)  GO TO 1350
      IF (nptu .LT. 4)  GO TO 1400
C
      nfcnco = nfcn
      nfcnmx = 100*(nptu+5)*(npar+1)
C           The minimum
      CALL mncuve(fcn,futil)
      u1min = u(ke1)
      u2min = u(ke2)
      ierrf = 0
      cfrom = chere
      nfcnfr = nfcnco
      IF (isw(5) .GE. 0)  THEN
         WRITE (isyswr,'(1X,A,I4,A)')
     +   'START MNCONTOUR CALCULATION OF',nptu,' POINTS ON CONTOUR.'
         IF (npar .GT. 2) THEN
            IF (npar .EQ. 3) THEN
              ki3 = 6 - ki1 - ki2
              ke3 = nexofi(ki3)
              WRITE (isyswr,'(1X,A,I3,2X,A)')
     +        'EACH POINT IS A MINIMUM WITH RESPECT TO PARAMETER ',
     +        ke3, cpnam(ke3)
            ELSE
              WRITE (isyswr,'(1X,A,I3,A)')
     +        'EACH POINT IS A MINIMUM WITH RESPECT TO THE OTHER',
     +        npar-2, ' VARIABLE PARAMETERS.'
            ENDIF
         ENDIF
      ENDIF
C
C           Find the first four points using MNMNOT
C              ........................ first two points
      CALL mnmnot(fcn,ke1,ke2,val2pl,val2mi,futil)
      IF (ern(ki1) .EQ. undefi)  THEN
         xptu(1) = alim(ke1)
         CALL mnwarn('W',chere,'Contour squeezed by parameter limits.')
      ELSE
         IF (ern(ki1) .GE. zero)  GO TO 1500
         xptu(1) = u1min+ern(ki1)
      ENDIF
      yptu(1) = val2mi
C
      IF (erp(ki1) .EQ. undefi)  THEN
         xptu(3) = blim(ke1)
         CALL mnwarn('W',chere,'Contour squeezed by parameter limits.')
      ELSE
         IF (erp(ki1) .LE. zero)  GO TO 1500
         xptu(3) = u1min+erp(ki1)
      ENDIF
      yptu(3) = val2pl
      scalx = 1.0/(xptu(3) - xptu(1))
C              ........................... next two points
      CALL mnmnot(fcn,ke2,ke1,val2pl,val2mi,futil)
      IF (ern(ki2) .EQ. undefi)  THEN
         yptu(2) = alim(ke2)
         CALL mnwarn('W',chere,'Contour squeezed by parameter limits.')
      ELSE
         IF (ern(ki2) .GE. zero)  GO TO 1500
         yptu(2) = u2min+ern(ki2)
      ENDIF
      xptu(2) = val2mi
      IF (erp(ki2) .EQ. undefi)  THEN
         yptu(4) = blim(ke2)
         CALL mnwarn('W',chere,'Contour squeezed by parameter limits.')
      ELSE
         IF (erp(ki2) .LE. zero)  GO TO 1500
         yptu(4) = u2min+erp(ki2)
      ENDIF
      xptu(4) = val2pl
      scaly = 1.0/(yptu(4) - yptu(2))
      nowpts = 4
      next = 5
      IF (ldebug) THEN
         WRITE (isyswr,'(A)') ' Plot of four points found by MINOS'
         xpt(1) = u1min
         ypt(1) = u2min
         chpt(1) = ' '
         nall = min(nowpts+1,maxcpt)
         DO 85 i= 2, nall
           xpt(i) = xptu(i-1)
           ypt(i) = yptu(i-1)
   85    CONTINUE
           chpt(2)= 'A'
           chpt(3)= 'B'
           chpt(4)= 'C'
           chpt(5)= 'D'
         CALL mnplot(xpt,ypt,chpt,nall,isyswr,npagwd,npagln)
      ENDIF
C
C               ..................... save some values before fixing
      isw2 = isw(2)
      isw4 = isw(4)
      sigsav = edm
      istrav = istrat
      dc = dcovar
      apsi  = epsi*0.5
      abest=amin
      mpar=npar
      nfmxin = nfcnmx
      DO 125 i= 1, mpar
  125 xt(i) = x(i)
      DO 130 j= 1, mpar*(mpar+1)/2
  130 vthmat(j) = vhmat(j)
      DO 135 i= 1, mpar
      gcc(i) = globcc(i)
  135 w(i) = werr(i)
C                           fix the two parameters in question
      kints = niofex(ke1)
      CALL mnfixp (kints,ierr)
      kints = niofex(ke2)
      CALL mnfixp (kints,ierr)
C               ......................Fill in the rest of the points
      DO 900 inew= next, nptu
C            find the two neighbouring points with largest separation
      bigdis = 0.
         DO 200  iold = 1, inew-1
         i2 = iold + 1
         IF (i2 .EQ. inew) i2 = 1
         dist = (scalx*(xptu(iold)-xptu(i2)))**2 +
     +          (scaly*(yptu(iold)-yptu(i2)))**2
         IF (dist .GT. bigdis) THEN
            bigdis = dist
            idist = iold
         ENDIF
  200    CONTINUE
      i1 = idist
      i2 = i1 + 1
      IF (i2 .EQ. inew) i2 = 1
C                   next point goes between I1 and I2
      a1 = half
      a2 = half
  300 xmidcr = a1*xptu(i1) + a2*xptu(i2)
      ymidcr = a1*yptu(i1) + a2*yptu(i2)
      xdir = yptu(i2) - yptu(i1)
      ydir = xptu(i1) - xptu(i2)
      sclfac = max(abs(xdir*scalx), abs(ydir*scaly))
      xdircr = xdir/sclfac
      ydircr = ydir/sclfac
      ke1cr = ke1
      ke2cr = ke2
C                Find the contour crossing point along DIR
      amin = abest
      CALL mncros(fcn,aopt,iercr,futil)
      IF (iercr .GT. 1)  THEN
C              If cannot find mid-point, try closer to point 1
         IF (a1 .GT. half) THEN
            IF (isw(5) .GE. 0)
     +      WRITE (isyswr,'(A,A,I3,A)') ' MNCONT CANNOT FIND NEXT',
     +           ' POINT ON CONTOUR.  ONLY ',nowpts,' POINTS FOUND.'
            GO TO 950
         ENDIF
         CALL mnwarn('W',chere,'Cannot find midpoint, try closer.')
         a1 = 0.75
         a2 = 0.25
         GO TO 300
      ENDIF
C                Contour has been located, insert new point in list
         DO 830 move= nowpts,i1+1,-1
         xptu(move+1) = xptu(move)
         yptu(move+1) = yptu(move)
  830    CONTINUE
      nowpts = nowpts + 1
      xptu(i1+1) = xmidcr + xdircr*aopt
      yptu(i1+1) = ymidcr + ydircr*aopt
  900 CONTINUE
  950 CONTINUE
C
      ierrf = nowpts
      cstatu = 'SUCCESSFUL'
      IF (nowpts .LT. nptu)  cstatu = 'INCOMPLETE'
C                make a lineprinter plot of the contour
      IF (isw(5) .GE. 0) THEN
         xpt(1) = u1min
         ypt(1) = u2min
         chpt(1) = ' '
         nall = min(nowpts+1,maxcpt)
         DO 1000 i= 2, nall
           xpt(i) = xptu(i-1)
           ypt(i) = yptu(i-1)
           chpt(i)= 'X'
 1000    CONTINUE
         WRITE (isyswr,'(A,I3,2X,A)') ' Y-AXIS: PARAMETER ',ke2,
     +        cpnam(ke2)
         CALL mnplot(xpt,ypt,chpt,nall,isyswr,npagwd,npagln)
         WRITE (isyswr,'(25X,A,I3,2X,A)') 'X-AXIS: PARAMETER ',
     +         ke1,cpnam(ke1)
      ENDIF
C                 print out the coordinates around the contour
      IF (isw(5) .GE. 1)  THEN
         npcol = (nowpts+1)/2
         nfcol = nowpts/2
         WRITE (isyswr,'(/I5,A,G13.5,A,G11.3)') nowpts,
     +    ' POINTS ON CONTOUR.   FMIN=',abest,'   ERRDEF=',up
         WRITE (isyswr,'(9X,A,3X,A,18X,A,3X,A)')
     +         cpnam(ke1),cpnam(ke2),cpnam(ke1),cpnam(ke2)
         DO 1050 line = 1, nfcol
           lr = line + npcol
           WRITE (isyswr,'(1X,I5,2G13.5,10X,I5,2G13.5)')
     +     line,xptu(line),yptu(line),lr,xptu(lr),yptu(lr)
 1050    CONTINUE
         IF (nfcol .LT. npcol) WRITE (isyswr,'(1X,I5,2G13.5)')
     +                         npcol,xptu(npcol),yptu(npcol)
      ENDIF
C                                    . . contour finished. reset v
      itaur = 1
      CALL mnfree(1)
      CALL mnfree(1)
      DO 1100 j= 1, mpar*(mpar+1)/2
 1100 vhmat(j) = vthmat(j)
      DO 1120 i= 1, mpar
      globcc(i) = gcc(i)
      werr(i) = w(i)
 1120 x(i) = xt(i)
      CALL mninex (x)
      edm = sigsav
      amin = abest
      isw(2) = isw2
      isw(4) = isw4
      dcovar = dc
      itaur = 0
      nfcnmx = nfmxin
      istrat = istrav
      u(ke1) = u1min
      u(ke2) = u2min
      GO TO 2000
C                                     Error returns
 1350 WRITE (isyswr,'(A)') ' INVALID PARAMETER NUMBERS.'
      GO TO 1450
 1400 WRITE (isyswr,'(A)') ' LESS THAN FOUR POINTS REQUESTED.'
 1450 ierrf = -1
      cstatu = 'USER ERROR'
      GO TO 2000
 1500 WRITE (isyswr,'(A)') ' MNCONT UNABLE TO FIND FOUR POINTS.'
      u(ke1) = u1min
      u(ke2) = u2min
      ierrf = 0
      cstatu = 'FAILED'
 2000 CONTINUE
      cfrom = chere
      nfcnfr = nfcnco
      RETURN

Here is the call graph for this function:

Here is the caller graph for this function: