MINICERN/mathlib/gen/divon/delslv.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:03:25  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE DELSLV(N,FMAJOR,FMINOR,LMAX,FRACT,X,XLOW,XUP,VOL,NCUT,N
  11      1CDIM,ICUT,DELPLS,DELNEG,REGTOL,FTOL,FORIG,DFORIG,FNEW,FNLIN,FNLROW
  12      2,DIAGJ,SPDIAG,SOL,Z,NFCNT)
  13       INTEGER N, NCUT, NCDIM, NFCNT
  14       INTEGER ICUT(NCDIM)
  15       DOUBLE PRECISION FMAJOR, FMINOR, FRACT, REGTOL, FTOL
  16       DOUBLE PRECISION X(N), XUP(N), XLOW(N), DELPLS(N), DELNEG(N)
  17       DOUBLE PRECISION FORIG(NCDIM), FNLIN(NCDIM)
  18       DOUBLE PRECISION DFORIG(NCDIM), FNEW(NCDIM), FNLROW(NCDIM)
  19       DOUBLE PRECISION DIAGJ(NCDIM), SPDIAG(NCDIM), SOL(NCDIM), Z(N)
  20       LOGICAL LMAX
  21       INTEGER I, IAB, IBACK, II, ISAVE, ITRY,
  22      1        NCUTM1, NNEAR
  23       DOUBLE PRECISION BIG, DEL, DELMAX, DELMIN, DFNEW, DFUN,
  24      1        FDIF, FGAM, FNORM, FNRMNW, FOMX, FRAT, FZ, GAMMA,
  25      2        GAMNEW, PROD, RATGAM, REGINV, SINGTL, TSTVAL,
  26      3        VAL, VOL, XMULT, YDI, YDMIN
  27 CMM   INTEGER IABS
  28       DATA SINGTL/ 1.0D-4/
  29       DATA BIG/ 1.0D+10/
  30       NFCNT=0
  31       FDIF=FMAJOR-FMINOR
  32       PROD=1.0D+0
  33       DO 10 I=1,N
  34       PROD=PROD*(DELPLS(I)+DELNEG(I))
  35  10   CONTINUE
  36       GAMMA=PROD/VOL
  37       FGAM=GAMMA*FMAJOR+(1.0D+0-GAMMA)*FMINOR
  38       DO 20 I=1,N
  39       Z(I)=X(I)
  40  20   CONTINUE
  41       DO 50 I=1,NCUT
  42       II=ICUT(I)
  43       IAB=ABS(II)
  44       IF(II.LT.0) GOTO 30
  45       Z(IAB)=X(IAB)+DELPLS(IAB)
  46       GOTO 40
  47  30   Z(IAB)=X(IAB)-DELNEG(IAB)
  48  40   FNEW(I)=DFUN(N,Z)
  49       Z(IAB)=X(IAB)
  50       NFCNT=NFCNT+1
  51  50   CONTINUE
  52 #if defined(CERNLIB_IBM)||defined(CERNLIB_SINGLE)
  53  60   YDMIN=1.0D+40
  54 #endif
  55 #if (!defined(CERNLIB_IBM))&&(defined(CERNLIB_DOUBLE))
  56    60 YDMIN= 1.0D+30
  57 #endif
  58       NNEAR=0
  59       DO 70 I=1,NCUT
  60       YDI=ABS(FMAJOR-FNEW(I))
  61       IF(YDI.GT.YDMIN) GOTO 70
  62       YDMIN=YDI
  63       NNEAR=I
  64  70   CONTINUE
  65       IF(NNEAR.EQ.0) RETURN
  66       IF((LMAX.AND.FNEW(NNEAR).LT.FGAM).OR.
  67      1(.NOT.LMAX.AND.FNEW(NNEAR).GT.FGAM)) GOTO 130
  68       ISAVE=ICUT(NNEAR)
  69       IF(NNEAR.EQ.NCUT.OR.NCUT.EQ.1) GOTO 90
  70       NCUTM1=NCUT-1
  71       DO 80 I=NNEAR,NCUTM1
  72       ICUT(I)=ICUT(I+1)
  73       FNEW(I)=FNEW(I+1)
  74  80   CONTINUE
  75  90   NCUT=NCUT-1
  76       IF(ISAVE.LT.0) GOTO 100
  77       DELPLS(ISAVE)=XUP(ISAVE)-X(ISAVE)
  78       GOTO 110
  79  100  IAB=ABS(ISAVE)
  80       DELNEG(IAB)=X(IAB)-XLOW(IAB)
  81  110  PROD=1.0D+0
  82       DO 120 I=1,N
  83       PROD=PROD*(DELPLS(I)+DELNEG(I))
  84  120  CONTINUE
  85       GAMMA=PROD/VOL
  86       FGAM=GAMMA*FMAJOR+(1.0D+0-GAMMA)*FMINOR
  87       IF(NCUT.EQ.0) RETURN
  88       GOTO 60
  89  130  DO 160 I=1,NCUT
  90       II=ICUT(I)
  91       IAB=ABS(II)
  92       IF(II.LT.0) GOTO 140
  93       DEL=DELPLS(IAB)
  94       GOTO 150
  95  140  DEL=DELNEG(IAB)
  96  150  FORIG(I)=FNEW(I)
  97       DFORIG(I)=(FNEW(I)-FMAJOR)/DEL
  98  160  CONTINUE
  99       CALL FEQN(NCUT,FORIG,FGAM,FNLIN)
 100       CALL RLEN(NCUT,FNLIN,FNORM)
 101  170  DO 180 I=1,NCUT
 102       FNLIN(I)=-FNLIN(I)
 103  180  CONTINUE
 104       IF(NCUT.EQ.1) GOTO 200
 105       DIAGJ(1)=DFORIG(1)
 106       SPDIAG(1)=-DFORIG(2)
 107       NCUTM1=NCUT-1
 108       DO 190 I=1,NCUTM1
 109       DIAGJ(I)=DFORIG(I)
 110       SPDIAG(I)=-DFORIG(I+1)
 111  190  CONTINUE
 112  200  DO 210 I=1,NCUT
 113       II=ICUT(I)
 114       II=ABS(II)
 115       FNLROW(I)=-GAMMA*FDIF/(DELPLS(II)+DELNEG(II))
 116  210  CONTINUE
 117       FNLROW(1)=DFORIG(1)+FNLROW(1)
 118       IF(NCUT.EQ.1) GOTO 230
 119       DO 220 I=1,NCUTM1
 120       XMULT=0.0D+0
 121       IF(ABS(FNLROW(I)).LT.BIG*ABS(DIAGJ(I))) XMULT=FNLROW(I)/DIAGJ(I
 122      1)
 123       FNLROW(I+1)=FNLROW(I+1)-XMULT*SPDIAG(I)
 124       FNLIN(NCUT)=FNLIN(NCUT)-XMULT*FNLIN(I)
 125  220  CONTINUE
 126  230  SOL(NCUT)=FNLIN(NCUT)
 127       IF(ABS(FNLROW(NCUT)).LT.BIG*ABS(FNLIN(NCUT))) SOL(NCUT)=FNLIN(N
 128      1CUT)/FNLROW(NCUT)
 129       IF(NCUT.EQ.1) GOTO 250
 130       DO 240 I=2,NCUT
 131       IBACK=NCUT-I+1
 132       VAL=FNLIN(IBACK)-SOL(IBACK+1)*SPDIAG(IBACK)
 133       SOL(IBACK)=VAL
 134       IF(ABS(DIAGJ(IBACK)).LT.BIG*ABS(VAL)) SOL(IBACK)=VAL/DIAGJ(IBAC
 135      1K)
 136  240  CONTINUE
 137  250  ITRY=0
 138       DO 280 I=1,NCUT
 139       II=ICUT(I)
 140       IAB=ABS(II)
 141       IF(II.LT.0) GOTO 260
 142       DELMAX=FRACT*(XUP(IAB)-X(IAB)-DELPLS(IAB))
 143       DELMIN=-DELPLS(IAB)
 144       GOTO 270
 145  260  DELMAX=FRACT*(X(IAB)-XLOW(IAB)-DELNEG(IAB))
 146       DELMIN=-DELNEG(IAB)
 147  270  IF(SOL(I).GT.DELMAX) SOL(I)=0.75D+0*DELMAX
 148       IF(SOL(I).LT.DELMIN) SOL(I)=0.75D+0*DELMIN
 149  280  CONTINUE
 150  290  DO 310 I=1,NCUT
 151       II=ICUT(I)
 152       IAB=ABS(II)
 153       IF(II.LT.0) GOTO 300
 154       DELPLS(IAB)=DELPLS(IAB)+SOL(I)
 155       GOTO 310
 156  300  DELNEG(IAB)=DELNEG(IAB)+SOL(I)
 157  310  CONTINUE
 158       DO 340 I=1,NCUT
 159       II=ICUT(I)
 160       IAB=ABS(II)
 161       IF(II.LT.0) GOTO 320
 162       Z(IAB)=X(IAB)+DELPLS(IAB)
 163       GOTO 330
 164  320  Z(IAB)=X(IAB)-DELNEG(IAB)
 165  330  FZ=DFUN(N,Z)
 166       NFCNT=NFCNT+1
 167       FNEW(I)=FZ
 168       Z(IAB)=X(IAB)
 169  340  CONTINUE
 170       PROD=1.0D+0
 171       DO 350 I=1,N
 172       PROD=PROD*(DELPLS(I)+DELNEG(I))
 173  350  CONTINUE
 174       GAMNEW=PROD/VOL
 175       FGAM=GAMNEW*FMAJOR+(1.0D+0-GAMNEW)*FMINOR
 176       CALL FEQN(NCUT,FNEW,FGAM,FNLIN)
 177       CALL RLEN(NCUT,FNLIN,FNRMNW)
 178       IF(FNRMNW.GT.FNORM) GOTO 380
 179       FOMX=0.0D+0
 180       DO 360 I=1,NCUT
 181       IF(ABS(FNEW(I)).GT.FOMX) FOMX=ABS(FNEW(I))
 182       DFNEW=FNEW(I)-FORIG(I)
 183       TSTVAL=1.0D+0
 184       IF(ABS(SOL(I)).LT.BIG*ABS(DFNEW)) TSTVAL=DFNEW/SOL(I)
 185       IF(ABS(TSTVAL).LT.SINGTL*ABS(DFORIG(I))) TSTVAL=SINGTL*DFORIG(I
 186      1)
 187       DFORIG(I)=TSTVAL
 188       FORIG(I)=FNEW(I)
 189  360  CONTINUE
 190       FNORM=FNRMNW
 191       FOMX=MAX(FOMX,ABS(FGAM))
 192       FRAT=FNORM/(1.0D+0+FOMX)
 193       REGINV=1.0D+0/REGTOL
 194       RATGAM=GAMNEW/GAMMA
 195       GAMMA=GAMNEW
 196 C---  Activate to do debugging
 197 C     WRITE(6,420) FNORM,FOMX,FRAT
 198 C     WRITE(6,430) FGAM,RATGAM
 199       FGAM=GAMMA*FMAJOR+(1.0D+0-GAMMA)*FMINOR
 200       IF(RATGAM.GT.REGTOL.AND.RATGAM.LT.REGINV) GOTO 370
 201       IF(FRAT.LT.FTOL) GOTO 370
 202       GOTO 170
 203  370  RETURN
 204  380  ITRY=ITRY+1
 205       IF(ITRY.GT.2) RETURN
 206       DO 410 I=1,NCUT
 207       II=ICUT(I)
 208       IAB=ABS(II)
 209       IF(II.GE.0) THEN
 210          DELPLS(IAB)=DELPLS(IAB)-SOL(I)
 211       ELSE
 212          DELNEG(IAB)=DELNEG(IAB)-SOL(I)
 213       ENDIF
 214       SOL(I)=SOL(I)*0.25D+0
 215  410  CONTINUE
 216       GOTO 290
 217  420  FORMAT(' FNORM, FOMX, FRAT', 3(1PD15.5))
 218  430  FORMAT(' FGAM, RATGAM', 2(1PD15.5))
 219       END