mnpsdf.f

Go to the documentation of this file.

*
* $Id: mnpsdf.F,v 1.2 1996/03/15 18:02:50 james Exp $
*
* $Log: mnpsdf.F,v $
* Revision 1.2  1996/03/15 18:02:50  james
*     Modified Files:
* mnderi.F eliminate possible division by zero
* mnexcm.F suppress print on STOP when print flag=-1
*          set FVAL3 to flag if FCN already called with IFLAG=3
* mninit.F set version 96.03
* mnlims.F remove arguments, not needed
* mnmigr.F VLEN -> LENV in debug print statement
* mnparm.F move call to MNRSET to after NPAR redefined, to zero all
* mnpsdf.F eliminate possible division by zero
* mnscan.F suppress printout when print flag =-1
* mnset.F  remove arguments in call to MNLIMS
* mnsimp.F fix CSTATU so status is PROGRESS only if new minimum
* mnvert.F eliminate possible division by zero
*
* Revision 1.1.1.1  1996/03/07 14:31:31  mclareni
* Minuit
*
*
      SUBROUTINE mnpsdf
*
* $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        calculates the eigenvalues of v to see if positive-def.
CC        if not, adds constant along diagonal to make positive.
*
* $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
      CHARACTER CHBUFF*12
      dimension s(mni)
      epsmin = 1.0e-6
      epspdf = max(epsmin, epsma2)
      dgmin = vhmat(1)
C                        Check if negative or zero on diagonal
      DO 200 i= 1, npar
      ndex = i*(i+1)/2
      IF (vhmat(ndex) .LE. zero) THEN
          WRITE (chbuff(1:3),'(I3)') i
          CALL mnwarn('W',cfrom,
     +'Negative diagonal element'//chbuff(1:3)//' in Error Matrix')
      ENDIF
      IF (vhmat(ndex) .LT. dgmin)  dgmin = vhmat(ndex)
  200 CONTINUE
      IF (dgmin .LE. zero) THEN
         dg = (one+epspdf) - dgmin
         WRITE (chbuff,'(E12.2)') dg
         CALL mnwarn('W',cfrom,
     +     chbuff//' added to diagonal of error matrix')
      ELSE
         dg = zero
      ENDIF
C                    Store VHMAT in P, make sure diagonal pos.
      DO 213 i= 1, npar
      ndex = i*(i-1)/2
      ndexd = ndex + i
      vhmat(ndexd) = vhmat(ndexd) + dg
      IF (vhmat(ndexd) .LE. zero)   vhmat(ndexd) = 1.0
      s(i) = 1.0/sqrt(vhmat(ndexd))
      DO 213 j= 1, i
      ndex =  ndex + 1
  213 p(i,j) = vhmat(ndex) * s(i)*s(j)
C      call eigen (p,p,maxint,npar,pstar,-npar)
      CALL mneig(p,maxint,npar,maxint,pstar,epspdf,ifault)
      pmin = pstar(1)
      pmax = pstar(1)
      DO 215 ip= 2, npar
      IF (pstar(ip) .LT. pmin)  pmin = pstar(ip)
      IF (pstar(ip) .GT. pmax)  pmax = pstar(ip)
  215 CONTINUE
      pmax = max(abs(pmax), one)
      IF ((pmin .LE. zero .AND. lwarn) .OR.  isw(5) .GE. 2) THEN
         WRITE (isyswr,550)
         WRITE (isyswr,551) (pstar(ip),ip=1,npar)
      ENDIF
      IF (pmin .GT. epspdf*pmax)  GO TO 217
      IF (isw(2) .EQ. 3)  isw(2)=2
      padd = 1.0e-3*pmax - pmin
      DO 216 ip= 1, npar
      ndex = ip*(ip+1)/2
  216 vhmat(ndex) = vhmat(ndex) *(1.0 + padd)
      cstatu= 'NOT POSDEF'
      WRITE (chbuff,'(G12.5)') padd
      CALL mnwarn('W',cfrom,
     +   'MATRIX FORCED POS-DEF BY ADDING '//chbuff//' TO DIAGONAL.')
  217 CONTINUE
C
  550 FORMAT (' EIGENVALUES OF SECOND-DERIVATIVE MATRIX:' )
  551 FORMAT (7x,6e12.4)
      RETURN
      END