C*********************************************************************
C
C ******************************************************************
C **** ****
C **** THIS IS SUBROUTINE STEPIT ****
C **** ****
C ******************************************************************
C For further information, you should see the article below.
C Also, you should cite the following article in publications
C that make use of this subroutine.
C Chandler, J. P. (1969). STEPIT: Finds local minima of a smooth
C function of several parameters (CPA 312). Behavioral Science,14,81-82.
C**********************************************************************
SUBROUTINE STEPIT
C COPYRIGHT 1965 -- J. P. CHANDLER, PHYSICS DEPT., INDIANA UNIVERSITY.
C
C STEPIT 5.2 AUGUST 1967
C AVAILABLE FROM.... QUANTUM CHEMISTRY PROGRAM EXCHANGE
C I.U. CHEMISTRY DEPT., BLOOMINGTON, INDIANA.
C
C MODIFIED FOR FORTRAN G OS/360 MVT
C
C IMPLICIT REAL*8 (A-H,O-Z)
IMPLICIT DOUBLEPRECISION(A-H,O-Z)
DIMENSION VEC(50),TRIAL(50),XSAVE(50),CHI(50),DX(50),SECOND(2,2)
DIMENSION OLDVEC(50),SALVO(50),XOSC(50,15),CHIOSC(15),NFLAT(50)
C
COMMON/STPIT/CHISQ,X(50),XMAX(50),XMIN(50),DELTAX(50),ERR(50,50)
+,DELMIN(50),MASK(50),NV,NTRACE,MATRIX
COMMON/FASTER/JVARY
COMMON/ICALL/ICALLS,IMAX
NVMAX=50
MOSQUE=15
KW=61
RATIO=10.0
COLIN=0.99
NCOMP=5
ACK=2.0
SIGNIF=2.E8
HUGE=1.E37
C
JVARY=0
40 IF(NV)290,290,50
50 NACTIV=0
DO 150I=1,NV
IF(MASK(I))150,60,150
60 IF(SIGNIF*DABS(DELTAX(I))-DABS(X(I)))70,70,100
70 IF(X(I))90,80,90
80 DELTAX(I)=0.01
GOTO 100
90 DELTAX(I)=0.01*X(I)
100 IF(DELMIN(I))120,110,120
110 DELMIN(I)=DELTAX(I)/SIGNIF
120 IF(XMAX(I)-XMIN(I))130,130,140
130 XMAX(I)=HUGE
XMIN(I)=-HUGE
140 NACTIV=NACTIV+1
VACD=DMIN1(XMAX(I),X(I))
X(I)=DMAX1(XMIN(I),VACD)
150 CONTINUE
COMPAR=0.0
IF(NACTIV-1)160,190,180
160 DO 170J=1,NV
170 MASK(J)=0
GOTO 50
180 A=NACTIV
SUB=2.0/(A-1.0)
P=2.0*(1.0/DSQRT(A)/(1.0-0.5**SUB)-1.0)
COMPAR=DMIN1(9.99D-1,DABS((1.0-(1.0-COLIN)**SUB)*(1.0+P*(1.-
+COLIN))))
190 IF(NTRACE)280,200,200
200 WRITE(6,210)
210 FORMAT(
+'1ENTER SUBROUTINE STEPIT. COPYRIGHT 1965 J. P. CHANDLER, PHYSICS
+ DEPT., INDIANA UNIVERSITY.'//' INITIAL VALUES....'/)
WRITE(6,220)(MASK(J),J=1,NV)
220 FORMAT(/' MASK = '10(I6,6X)/(4X10I12))
WRITE(6,230)(X(J),J=1,NV)
230 FORMAT(/' X = '10E12.4/(10X10E12.4))
WRITE(6,240)(XMAX(J),J=1,NV)
240 FORMAT(/' XMAX = '10E12.4/(10X10E12.4))
WRITE(6,250)(XMIN(J),J=1,NV)
250 FORMAT(/' XMIN = '10E12.4/(10X10E12.4))
WRITE(6,260)(DELTAX(J),J=1,NV)
260 FORMAT(/' DELTAX = '10E12.4/(10X10E12.4))
WRITE(6,270)(DELMIN(J),J=1,NV)
270 FORMAT(/' DELMIN = '10E12.4/(10X10E12.4))
280 CALL FUNK
NF=1
JOCK=1
290 IF(NTRACE)320,300,300
300 WRITE(6,310)NV,NACTIV,MATRIX,NCOMP,RATIO,ACK,COLIN,COMPAR,CHISQ
310 FORMAT(//' 'I3,' VARIABLES,'I3,' ACTIVE.'10X8HMATRIX =I4,
+10X7H NCOMP=I2//' RATIO ='F5.1,10X5HACK =F5.1,10X7HCOLIN =F6.3,
+10X8HCOMPAR =F6.3///' CHISQ ='E15.8)
320 IF(NV)2150,2150,330
330 IF(NTRACE)360,360,340
340 WRITE(6,350)
350 FORMAT(//60(1X1H*)//10X37HTRACE MAP OF THE MINIMIZATION PROCESS//)
C
360 DO 370I=1,NV
DX(I)=DELTAX(I)
VEC(I)=0.
DO 370J=1,20
370 ERR(I,J)=0.
CHIOLD=CHISQ
NOSC=0
C
380 NCIRC=0
NZIP=0
C
C MAIN DO LOOP FOR CYCLING THROUGH THE VARIABLES.
C FIRST TRIAL STEP WITH EACH VARIABLE IS SEPARATE.
390 NACK=0
DO 1350I=1,NV
OLDVEC(I)=VEC(I)
VEC(I)=0.0
TRIAL(I)=0.0
IF(MASK(I))400,410,400
400 VEC(I)=-0.0
NFLAT(I)=1
GOTO 1350
410 NACK=NACK+1
SAVE=X(I)
IF(SIGNIF*DABS(DX(I))-DABS(X(I)))580,580,420
420 X(I)=SAVE+DX(I)
JVARY=0
IF(JOCK)440,440,430
430 JOCK=0
JVARY=I
440 NFLAG=1
IF(X(I)-XMIN(I))460,460,450
450 IF(X(I)-XMAX(I))470,460,460
460 NFLAG=NFLAG+3
GOTO 490
470 CALL FUNK
IF(ICALLS.LE.IMAX)GOTO10
RETURN
10 CONTINUE
NF=NF+1
JVARY=I
CHIME=CHISQ
IF(CHISQ-CHIOLD)620,480,490
480 NFLAG=NFLAG+1
490 X(I)=SAVE-DX(I)
IF(X(I)-XMIN(I))590,590,500
500 IF(X(I)-XMAX(I))510,590,590
510 CALL FUNK
NF=NF+1
JVARY=I
IF(CHISQ-CHIOLD)610,520,530
520 NFLAG=NFLAG+1
530 IF(NFLAG-3)540,580,590
540 IF((CHISQ-CHIME)*(CHIME-2.*CHIOLD+CHISQ))550,590,550
550 TRIAL(I)=.5*DX(I)*(CHISQ-CHIME)/(CHIME-2.*CHIOLD+CHISQ)
VEC(I)=TRIAL(I)/DABS(DX(I))
NFLAT(I)=0
X(I)=SAVE+TRIAL(I)
CALL FUNK
NF=NF+1
IF(CHISQ-CHIOLD)560,570,570
560 CHIOLD=CHISQ
JOCK=1
GOTO 600
570 TRIAL(I)=0.0
VEC(I)=0.0
GOTO 590
580 VEC(I)=-0.0
NFLAT(I)=1
590 X(I)=SAVE
600 NCIRC=NCIRC+1
IF(NCIRC-NACTIV)690,1430,1430
610 DX(I)=-DX(I)
C
C A LOWER VALUE HAS BEEN FOUND. HENCE THIS VARIABLE WILL CHANGE.
C
620 NCIRC=0
DEL=DX(I)
630 CHIME=CHIOLD
CHIOLD=CHISQ
VEC(I)=VEC(I)+DEL/DABS(DX(I))
NFLAT(I)=0
TRIAL(I)=TRIAL(I)+DEL
DEL=ACK*DEL
SAVE=X(I)
X(I)=SAVE+DEL
IF(X(I)-XMIN(I))680,680,640
640 IF(X(I)-XMAX(I))650,680,680
650 CALL FUNK
NF=NF+1
IF(CHISQ-CHIOLD)630,660,660
660 ZZZ=ACK*CHIME-(ACK+1.0)*CHIOLD+CHISQ
IF(ZZZ.EQ.0.)GOTO661
CINDER=(0.5/ACK)*(ACK**2*CHIME-(ACK**2-1.0)*CHIOLD-CHISQ)/ZZZ
GOTO 662
661 CINDER=0.0
662 X(I)=SAVE+CINDER*DEL
CALL FUNK
NF=NF+1
IF(CHISQ-CHIOLD)670,680,680
670 CHIOLD=CHISQ
TRIAL(I)=TRIAL(I)+CINDER*DEL
VEC(I)=VEC(I)+CINDER*DEL/DABS(DX(I))
GOTO 690
680 X(I)=SAVE
690 IF(NZIP-1)1340,700,700
700 IF(DABS(VEC(I))-ACK)750,710,710
710 DX(I)=ACK*DABS(DX(I))
VEC(I)=VEC(I)/ACK
OLDVEC(I)=OLDVEC(I)/ACK
DO 720J=1,MOSQUE
720 ERR(I,J)=ERR(I,J)/ACK
IF(NTRACE)750,750,730
730 WRITE(6,740)I,DX(I)
740 FORMAT(' STEP SIZE'I3,' INCREASED TO 'E12.5)
750 SUMO=0.0
SUMV=0.0
DO 760J=1,NV
SUMO=SUMO+OLDVEC(J)**2
760 SUMV=SUMV+VEC(J)**2
IF(SUMO*SUMV)1340,1340,770
770 SUMO=DSQRT(SUMO)
SUMV=DSQRT(SUMV)
COSINE=0.0
DO 780J=1,NV
780 COSINE=COSINE+(OLDVEC(J)/SUMO)*(VEC(J)/SUMV)
IF(NZIP-1)1340,790,800
790 IF(NACK-NACTIV)1340,820,820
800 IF(NACK-NACTIV)820,810,810
810 IF(NZIP-NCOMP)820,830,830
820 IF(COSINE-COMPAR)1340,830,830
C
C SIMON SAYS, TAKE AS MANY GIANT STEPS AS POSSIBLE...
C
830 IF(NTRACE)860,860,840
840 WRITE(6,850)CHIOLD,(VEC(J),J=1,I)
850 FORMAT(' CHISQ ='E15.8,5X14HNO. OF STEPS =10F9.2/(42X10F9.2))
860 NGIANT=0
NTRY=0
NRETRY=0
KL=1
NOSC=NOSC+1
IF(NOSC-MOSQUE)890,890,870
870 NOSC=MOSQUE
DO 880K=2,MOSQUE
CHIOSC(K-1)=CHIOSC(K)
DO 880J=1,NV
XOSC(J,K-1)=XOSC(J,K)
880 ERR(J,K-1)=ERR(J,K)
890 DO 900J=1,NV
XOSC(J,NOSC)=X(J)
900 ERR(J,NOSC)=VEC(J)/SUMV
CHIOSC(NOSC)=CHIOLD
IF(NOSC-3)960,910,910
C
C SEARCH FOR A PREVIOUS SUCCESSFUL GIANT STEP IN A DIRECTION MORE
C NEARLY PARALLEL TO THE DIRECTION OF THE PROPOSED STEP THAN WAS THE
C IMMEDIATELY PREVIOUS ONE.
C
910 COXCOM=0.0
DO 920J=1,NV
920 COXCOM=COXCOM+ERR(J,NOSC)*ERR(J,NOSC-1)
NAH=NOSC-2
930 NTRY=0
DO 950K=KL,NAH
NRETRY=NAH-K
COSINE=0.0
DO 940J=1,NV
940 COSINE=COSINE+ERR(J,NOSC)*ERR(J,K)
IF(COSINE-COXCOM)950,950,970
950 CONTINUE
960 CHIBAK=CHI(I)
GOTO 1020
970 NTRY=1
KL=K+1
IF(NTRACE)1000,1000,980
980 NT=NOSC-K
WRITE(6,990)NT
990 FORMAT(/1X8H********5X34HPOSSIBLEOSCILLATIONWITHPERIODI2,
+' DETECTED.')
1000 DO 1010J=1,NV
SALVO(J)=TRIAL(J)
1010 TRIAL(J)=(X(J)-XOSC(J,K))/ACK
CHIBAK=CHIOLD+(CHIOSC(K)-CHIOLD)/ACK
C
1020 DO 1040J=1,NV
XSAVE(J)=X(J)
TRIAL(J)=ACK*TRIAL(J)
IF(MASK(J))1040,1030,1040
1030 TRAAL=X(J)+TRIAL(J)
TRAIL=DMIN1(TRAAL,XMAX(J))
X(J)=DMAX1(TRAIL,XMIN(J))
1040 CONTINUE
JOCK=0
JVARY=0
CALL FUNK
NF=NF+1
IF(CHISQ-CHIOLD)1050,1080,1080
1050 CHIBAK=CHIOLD
CHIOLD=CHISQ
NGIANT=NGIANT+1
IF(NTRACE)1020,1020,1060
1060 WRITE(6,1070)CHISQ,(X(J),J=1,NV)
1070 FORMAT(' CHISQ ='E15.8/' X(I)....'/(10(1XE12.5)))
GOTO 1020
C
1080 IF(NRETRY)1100,1100,1090
1090 IF(NGIANT)1150,1150,1100
1100 CINDER=(0.5/ACK)*(ACK**2*CHIBAK-(ACK**2-1.0)*CHIOLD-CHISQ)/(ACK*
+CHIBAK-(ACK+1.0)*CHIOLD+CHISQ)
DO 1120J=1,NV
IF(MASK(J))1120,1110,1120
1110 CANDER=XSAVE(J)+CINDER*TRIAL(J)
XIJ=DMIN1(CANDER,XMAX(J))
X(J)=DMAX1(XIJ,XMIN(J))
1120 CONTINUE
JOCK=0
JVARY=0
CALL FUNK
NF=NF+1
IF(CHISQ-CHIOLD)1280,1130,1130
1130 IF(NGIANT)1170,1140,1170
1140 IF(NTRY)1150,1170,1150
1150 DO 1160J=1,NV
TRIAL(J)=SALVO(J)
1160 X(J)=XSAVE(J)
GOTO 1190
1170 DO 1180J=1,NV
TRIAL(J)=TRIAL(J)/ACK
1180 X(J)=XSAVE(J)
1190 IF(NTRACE)1240,1240,1200
1200 WRITE(6,1210)CHIOLD,NGIANT
1210 FORMAT(/' CHISQ ='E15.8,' AFTER'I3,' GIANT STEPS.')
WRITE(6,1220)(X(J),J=1,NV)
1220 FORMAT(' X(I)....'/(10(1XE12.5)))
WRITE(6,1230)
1230 FORMAT(/)
1240 IF(NGIANT)1250,1250,1310
1250 IF(NRETRY)1260,1260,930
1260 IF(NTRY)1270,1330,1270
1270 NTRY=0
GOTO 960
C
1280 CHIOLD=CHISQ
JOCK=1
IF(NTRACE)1310,1310,1290
1290 STEPS=FLOAT(NGIANT)+CINDER
WRITE(6,1300)CHIOLD,STEPS
1300 FORMAT(/' CHISQ ='E15.8,' AFTER'F6.1,' GIANT STEPS.')
WRITE(6,1220)(X(J),J=1,NV)
WRITE(6,1230)
1310 IF(NTRY)1320,380,1320
1320 NOSC=0
GOTO 380
1330 NOSC=MAX0(NOSC-1,0)
1340 CHI(I)=CHIOLD
1350 CONTINUE
C
C ANOTHER CYCLE THROUGH THE VARIABLES HAS BEEN COMPLETED.
C PRINT ANOTHER LINE OF TRACES.
C
IF(NTRACE)1370,1370,1360
1360 WRITE(6,850)CHIOLD,(VEC(J),J=1,NV)
1370 CONTINUE
1380 IF(NZIP)1420,1390,1420
1390 IF(NTRACE)1420,1420,1400
1400 WRITE(6,1220)(X(J),J=1,NV)
WRITE(6,1410)
1410 FORMAT(' ')
1420 NZIP=NZIP+1
GOTO 390
C
C A MINIMUM HAS BEEN FOUND. PRINT THE REMAINING TRACES.
C
1430 IF(NTRACE)1450,1450,1440
1440 WRITE(6,850)CHIOLD,(VEC(J),J=1,I)
1450 IF(NTRACE)1470,1470,1460
1460 WRITE(6,1220)(X(J),J=1,NV)
WRITE(6,1230)
C
C DECREASE THE SIZE OF THE STEPS FOR ALL VARIABLES.
C
1470 CONTINUE
1480 NOSC=0
NGATE=1
DO 1520J=1,NV
IF(MASK(J))1520,1490,1520
1490 IF(NFLAT(J))1500,1500,1520
1500 IF(DABS(DX(J))-DABS(DELMIN(J)))1520,1520,1510
1510 NGATE=0
1520 DX(J)=DX(J)/RATIO
IF(NGATE)1530,1530,1600
1530 IF(NTRACE)1570,1570,1540
1540 WRITE(6,1550)(DX(J),J=1,NV)
1550 FORMAT(60(1X1H*)//' STEP SIZES REDUCED TO....'//(10(1XE12.5)))
WRITE(6,1560)
1560 FORMAT(//)
1570 GOTO 380
1580 WRITE(6,1590)(DX(J),J=1,NV)
1590 FORMAT(' OPERAT. TERM.'//10(1X,E12.5))
C
1600 CHISQ=CHIOLD
IF(NTRACE)1630,1610,1610
1610 WRITE(6,1620)NF
1620 FORMAT(//1X,I5,' FUNCTION COMPUTATIONS ')
1630 CONTINUE
1640 IF(IABS(MATRIX-100)-50)1650,1650,2160
1650 IF(NACTIV-NV)2160,1660,2160
C
C COMPUTE THE STANDARD ERRORS AND THE CORRELATIONS.
C
1660 FAC=RATIO**(MATRIX-100)
DO 1680I=1,NV
DX(I)=DABS(FAC*DX(I))
XSAVE(I)=X(I)
JVARY=0
DO 1670J=1,2
X(I)=XSAVE(I)+DX(I)
CALL FUNK
NF=NF+1
JVARY=I
SECOND(1,J)=CHISQ
1670 DX(I)=-DX(I)
X(I)=XSAVE(I)
1680 ERR(I,I)=(SECOND(1,1)-2.0*CHIOLD+SECOND(1,2))/DX(I)**2
DO 1710I=2,NV
IM=I-1
DO 1710J=1,IM
DO 1700K=1,2
X(I)=XSAVE(I)+DX(I)
JVARY=0
DO 1690L=1,2
X(J)=XSAVE(J)+DX(J)
CALL FUNK
NF=NF+1
JVARY=J
SECOND(K,L)=CHISQ
X(J)=XSAVE(J)
1690 DX(J)=-DX(J)
X(I)=XSAVE(I)
1700 DX(I)=-DX(I)
ERR(I,J)=0.25*(SECOND(1,1)-SECOND(1,2)-SECOND(2,1)+SECOND(2,2))/
+DABS(DX(I)*DX(J))
1710 ERR(J,I)=ERR(I,J)
IF(NTRACE)1770,1720,1720
1720 WRITE(6,1730)
1730 FORMAT('1SIZES OF INCREMENTS TO BE USED BELOW....')
WRITE(6,1740)(DX(J),J=1,NV)
1740 FORMAT(/(10(1XE12.5)))
WRITE(6,1750)
1750 FORMAT(/////' MATRIX OF THE SECOND PARTIAL DERIVATIVES....'/)
DO 1760I=1,NV
1760 WRITE(6,1740)(ERR(I,J),J=1,I)
1770 DO 1780I=1,NV
DO 1780J=1,I
IF(ERR(I,J))1780,1790,1780
1780 CONTINUE
GOTO 1810
1790 WRITE(6,1800)
1800 FORMAT(////
+' THE ABOVE MATRIX CONTAINS ONE OR MORE ZEROES. A LARGER VALUE OF
+ -MATRIX- SHOULD BE TRIED, TO SEE IF THEY ARE LEGITIMATE. ')
C
C INVERT THE MATRIX USING SYMINV2 (COMM. OF THE A.C.M. 6, P. 67).
C
1810 DET=1.0
DETLOG=0.
DO 1820J=1,NV
1820 SALVO(J)=1.0
DO 1970I=1,NV
BIGAJJ=0.0
DO 1850J=1,NV
IF(SALVO(J))1830,1850,1830
1830 IF(DABS(ERR(J,J))-BIGAJJ)1850,1850,1840
1840 BIGAJJ=DABS(ERR(J,J))
K=J
1850 CONTINUE
IF(BIGAJJ)1870,1860,1870
1860 DET=0.0
GOTO 1980
1870 SALVO(K)=0.0
IF(DET.GT.1.E35.OR.DET.LT.1.E-35)GOTO1871
GOTO 1875
1871 WRITE(6,1872)
1872 FORMAT(///,'DETERMINANT COMPUTATION WILL CAUSE OVERFLOW.')
GOTO 2160
1875 DET=DET*ERR(K,K)
DETLOG=DETLOG+DLOG(DABS(ERR(K,K)))/2.303
TRIAL(K)=1.0/ERR(K,K)
ERR(K,K)=0.0
XSAVE(K)=1.0
M=K-1
IF(M)1910,1910,1880
1880 DO 1900J=1,M
XSAVE(J)=ERR(K,J)
TRIAL(J)=ERR(K,J)*TRIAL(K)
IF(SALVO(J))1860,1900,1890
1890 TRIAL(J)=-TRIAL(J)
1900 ERR(K,J)=0.0
1910 M=K+1
IF(M-NV)1920,1920,1960
1920 DO 1950J=M,NV
XSAVE(J)=ERR(J,K)
IF(SALVO(J))1860,1930,1940
1930 XSAVE(J)=-XSAVE(J)
1940 TRIAL(J)=-ERR(J,K)*TRIAL(K)
1950 ERR(J,K)=0.0
1960 DO 1970J=1,NV
DO 1970K=J,NV
1970 ERR(K,J)=ERR(K,J)+XSAVE(J)*TRIAL(K)
IF(DET)2000,1980,2020
1980 WRITE(6,1990)
1990 FORMAT(////
+' ERROR MATRIX IS SINGULAR. -MATRIX- SHOULD PROBABLY BE INCREASED
+. '/////)
GOTO 2150
2000 WRITE(6,2010)
2010 FORMAT(////
+' ERROR MATRIX IS NEGATIVE DEFINITE. -MATRIX- SHOULD PROBABLY BE
+DECREASED.')
2020 IF(NTRACE)2050,2030,2030
2030 WRITE(6,2040)DET,DETLOG
2040 FORMAT(////' DETERMINANT OF ABOVE MATRIX = 'E12.5, E12.5)
C +10X14HLOG10F(DET)=E12.5)
C ****** THE ABOVE LINE WAS COMMENTED OUT TO MAKE THE
C ****** COMPILIER HAPPY. IT BELONGS WITH LINE 2040
C ***************************************************
C
2050 DO 2090I=1,NV
DO 2060J=1,I
ERR(I,J)=ERR(I,J)*2.0
2060 ERR(J,I)=ERR(I,J)
IF(ERR(I,I))2070,2070,2090
2070 WRITE(6,2080)ERR(I,I)
2080 FORMAT(///' NEGATIVE OR ZERO MEAN SQUARE ERROR ENCOUNTERED...',
+3X,E15.8/' -MATRIX- SHOULD PROBABLY BE DECREASED.'///)
2090 XSAVE(I)=DSIGN(DSQRT(DABS(ERR(I,I))),ERR(I,I))
IF(NTRACE)2160,2100,2100
2100 WRITE(6,2110)
2110 FORMAT(/////' STANDARD ERRORS....')
WRITE(6,1740)(XSAVE(J),J=1,NV)
WRITE(6,2120)
2120 FORMAT(/////' LOWER TRIANGLE OF THE CORRELATION MATRIX....'/)
DO 2140I=2,NV
IM=I-1
DO 2130J=1,IM
2130 TRIAL(J)=ERR(I,J)/DABS(XSAVE(I)*XSAVE(J))
2140 WRITE(6,1740)(TRIAL(J),J=1,IM)
2150 WRITE(6,1620)NF
C
2160 CONTINUE
2190 JVARY=0
CALL FUNK
IF(NTRACE)2230,2200,2200
2200 WRITE(6,2210)(X(J),J=1,NV)
2210 FORMAT(///10X24HFINAL VALUES OF X(I)....//(7(1XE16.9)))
WRITE(6,2220)CHISQ
2220 FORMAT(//' FINAL VALUE OF CHISQ = 'E15.8//)
2230 RETURN
END