MINICERN/mathlib/gen/divon/locsch.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:03:27  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10       SUBROUTINE LOCSCH (IIMIN,N,NFREE,IFREE,X,FX,XLOW,XHI,DELTA,SUFTOL,
  11      1NFCNT,Y,FY,Z,FZ,P)
  12       INTEGER IIMIN, N, NFREE, NFCNT
  13       INTEGER IFREE(N)
  14       DOUBLE PRECISION FX, DELTA, SUFTOL, FY, FZ
  15       DOUBLE PRECISION X(N), XLOW(N), XHI(N), Y(N), Z(N), P(N)
  16       INTEGER I, II, ILOC, ITEST
  17       DOUBLE PRECISION A, B, B1, D, DELX, DFUN, E, EPSMCH, ETA, FA,
  18      1                 FBEST, FSAV, FTEST, FU, FV, FW, GTEST1, GTEST2,
  19      2                 GTP, GU, OLDF, PE, PNORM, R, RR, RTEPS,
  20      3                 SCXBD, SFTBND, SMAX, SNMAX, SOPP, SPMAX, SS,
  21      4                 SSAV, STP, STPDN, STPNEG, STPPOS, STPUP, T,
  22      5                 TOL, U, XBEST, XLAMDA, XV, XW, YSAV
  23       DATA EPSMCH/2.22D-16/
  24       RTEPS=SQRT(EPSMCH)
  25       DELX=DELTA
  26       DO 10 I=1,N
  27       Y(I)=X(I)
  28       Z(I)=X(I)
  29  10   CONTINUE
  30       FZ=FX
  31       FY=FX
  32       SMAX=1.0D+30
  33       DO 40 I=1,NFREE
  34       II=IFREE(I)
  35       STPUP=XHI(II)-X(II)
  36       STPDN=X(II)-XLOW(II)
  37       IF(STPUP.LT.STPDN) GOTO 20
  38       P(I)=1.0D+0
  39       STP=STPUP
  40       GOTO 30
  41  20   P(I)=-1.0D+0
  42       STP=STPDN
  43  30   IF(STP.LT.SMAX) SMAX=STP
  44  40   CONTINUE
  45       SMAX=0.9D+0*SMAX
  46  50   IF(DELX.GT.SMAX) DELX=SMAX
  47       DO 60 I=1,NFREE
  48       II=IFREE(I)
  49       Y(II)=X(II)+DELX*P(I)
  50  60   CONTINUE
  51       FY=DFUN(N,Y)
  52       IF(IIMIN.EQ.2) FY=-FY
  53       NFCNT=NFCNT+1
  54       FTEST=SUFTOL*(1.0D+0+ABS(FX))
  55       IF(ABS(FX-FY).GT.FTEST.OR.DELX.EQ.SMAX) GOTO 70
  56       DELX=5.0D+0*DELX
  57       GOTO 50
  58  70   DELX=DELTA
  59       CALL ORTHVC(N,NFREE,IFREE,X,FX,XLOW,XHI,Y,FY,P,SPMAX,SNMAX)
  60       SMAX=0.9D+0*SPMAX
  61       SOPP=0.9D+0*SNMAX
  62       IF(SMAX.GE.SNMAX.OR.NFREE.EQ.1) GOTO 90
  63       DO 80 I=1,NFREE
  64       P(I)=-P(I)
  65  80   CONTINUE
  66       SSAV=SMAX
  67       SMAX=SOPP
  68       SOPP=SSAV
  69  90   IF(DELX.GT.SMAX) DELX=SMAX
  70  100  DO 110 I=1,NFREE
  71       II=IFREE(I)
  72       Z(II)=Y(II)+DELX*P(I)
  73  110  CONTINUE
  74       FZ=DFUN(N,Z)
  75       IF(IIMIN.EQ.2) FZ=-FZ
  76       NFCNT=NFCNT+1
  77       FTEST=SUFTOL*(1.0D+0+ABS(FY))
  78       IF(ABS(FY-FZ).GT.FTEST.OR.(5.0D+0*DELX).GT.SMAX) GOTO 120
  79       DELX=5.0D+0*DELX
  80       GOTO 100
  81  120  IF(FY.EQ.FZ) GOTO 210
  82       IF(FY.LT.FZ) GOTO 130
  83       GTP=(FZ-FY)/DELX
  84       XLAMDA=SMAX/0.9D+0
  85       U=MIN(2.0D+0*DELX,SMAX)
  86       GOTO 160
  87  130  DO 140 I=1,NFREE
  88       P(I)=-P(I)
  89  140  CONTINUE
  90       U=MIN(2.0D+0*DELX,SOPP)
  91       XLAMDA=SOPP/0.9D+0+DELX
  92       GTP=(FY-FZ)/DELX
  93       DO 150 I=1,N
  94       YSAV=Y(I)
  95       Y(I)=Z(I)
  96       Z(I)=YSAV
  97  150  CONTINUE
  98       FSAV=FY
  99       FY=FZ
 100       FZ=FSAV
 101  160  CALL RLEN(NFREE,P,PNORM)
 102       PE=PNORM+RTEPS
 103       ILOC=1
 104       FU=FY
 105       GU=GTP
 106       SFTBND=0.0D+0
 107       ETA=RTEPS
 108       T=RTEPS/PE
 109  170  CALL NEWPTQ(RTEPS,T,ETA,SFTBND,XLAMDA,U,FU,GU,XBEST,FBEST,XW,FW,X
 110      1V,FV,A,FA,B,OLDF,B1,SCXBD,E,D,RR,SS,GTEST1,GTEST2,TOL,ILOC,ITEST)
 111       IF(ITEST.NE.1) GOTO 190
 112       R=XBEST+U
 113       DO 180 I=1,NFREE
 114       II=IFREE(I)
 115       Z(II)=Y(II)+R*P(I)
 116  180  CONTINUE
 117       FU=DFUN(N,Z)
 118       IF(IIMIN.EQ.2) FU=-FU
 119       NFCNT=NFCNT+1
 120       FZ=FU
 121       GOTO 170
 122  190  IF(ITEST.NE.0) GOTO 210
 123       DO 200 I=1,NFREE
 124       II=IFREE(I)
 125       Z(II)=Y(II)+XBEST*P(I)
 126  200  CONTINUE
 127       FZ=FBEST
 128  210  IF(FZ.EQ.FX) RETURN
 129       DO 220 I=1,NFREE
 130       II=IFREE(I)
 131       P(I)=Z(II)-X(II)
 132  220  CONTINUE
 133       SPMAX=1.0D+30
 134       SNMAX=1.0D+30
 135       DO 230 I=1,NFREE
 136       II=IFREE(I)
 137       CALL MXSTEP(X(II),XLOW(II),XHI(II),P(I),STPPOS,STPNEG)
 138       IF(STPPOS.LT.SPMAX) SPMAX=STPPOS
 139       IF(STPNEG.LT.SNMAX) SNMAX=STPNEG
 140  230  CONTINUE
 141       CALL RLEN(NFREE,P,PNORM)
 142       GTP=FZ-FX
 143       U=MIN(2.0D+0,0.9D+0*SPMAX)
 144       XLAMDA=SPMAX
 145       IF(FZ.LT.FX) GOTO 270
 146       DELX=DELTA/(PNORM+RTEPS)
 147       IF(DELX.GT.0.9D+0*SNMAX) DELX=0.9D+0*SNMAX
 148       DO 240 I=1,NFREE
 149       II=IFREE(I)
 150       Y(II)=X(II)-DELX*P(I)
 151  240  CONTINUE
 152       FY=DFUN(N,Y)
 153       IF(IIMIN.EQ.2) FY=-FY
 154       NFCNT=NFCNT+1
 155       IF(FY.LT.FX) GOTO 250
 156       IF(DELX.LT.1.0D+0) GTP=(FX-FY)/DELX
 157       GOTO 270
 158  250  GTP=(FY-FX)/DELX
 159       DO 260 I=1,NFREE
 160       P(I)=-P(I)
 161  260  CONTINUE
 162       U=MIN(2.0D+0*DELX,0.9D+0*SNMAX)
 163       XLAMDA=SNMAX
 164  270  ILOC=1
 165       PE=PNORM+RTEPS
 166       FU=FX
 167       GU=GTP
 168       SFTBND=0.0D+0
 169       ETA=RTEPS
 170       T=RTEPS/PE
 171  280  CALL NEWPTQ(RTEPS,T,ETA,SFTBND,XLAMDA,U,FU,GU,XBEST,FBEST,XW,FW,X
 172      1V,FV,A,FA,B,OLDF,B1,SCXBD,E,D,RR,SS,GTEST1,GTEST2,TOL,ILOC,ITEST)
 173       IF(ITEST.NE.1) GOTO 300
 174       R=XBEST+U
 175       DO 290 I=1,NFREE
 176       II=IFREE(I)
 177       Y(II)=X(II)+R*P(I)
 178  290  CONTINUE
 179       FU=DFUN(N,Y)
 180       IF(IIMIN.EQ.2) FU=-FU
 181       NFCNT=NFCNT+1
 182       GOTO 280
 183  300  IF(ITEST.NE.0) RETURN
 184       DO 310 I=1,NFREE
 185       II=IFREE(I)
 186       Z(II)=X(II)+XBEST*P(I)
 187  310  CONTINUE
 188       FZ=FBEST
 189       RETURN
 190       END