1 SUBROUTINE PRELIM (N,NPT,X,XL,XU,RHOBEG,IPRINT,MAXFUN,XBASE, 2 1 XPT,FVAL,GOPT,HQ,PQ,BMAT,ZMAT,NDIM,SL,SU,NF,KOPT) 3 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 4 DIMENSION X(*),XL(*),XU(*),XBASE(*),XPT(NPT,*),FVAL(*),GOPT(*), 5 1 HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),SL(*),SU(*) 6C 7C The arguments N, NPT, X, XL, XU, RHOBEG, IPRINT and MAXFUN are the 8C same as the corresponding arguments in SUBROUTINE BOBYQA. 9C The arguments XBASE, XPT, FVAL, HQ, PQ, BMAT, ZMAT, NDIM, SL and SU 10C are the same as the corresponding arguments in BOBYQB, the elements 11C of SL and SU being set in BOBYQA. 12C GOPT is usually the gradient of the quadratic model at XOPT+XBASE, but 13C it is set by PRELIM to the gradient of the quadratic model at XBASE. 14C If XOPT is nonzero, BOBYQB will change it to its usual value later. 15C NF is maintaned as the number of calls of CALFUN so far. 16C KOPT will be such that the least calculated value of F so far is at 17C the point XPT(KOPT,.)+XBASE in the space of the variables. 18C 19C SUBROUTINE PRELIM sets the elements of XBASE, XPT, FVAL, GOPT, HQ, PQ, 20C BMAT and ZMAT for the first iteration, and it maintains the values of 21C NF and KOPT. The vector X is also changed by PRELIM. 22C 23C Set some constants. 24C 25 HALF=0.5D0 26 ONE=1.0D0 27 TWO=2.0D0 28 ZERO=0.0D0 29 RHOSQ=RHOBEG*RHOBEG 30 RECIP=ONE/RHOSQ 31 NP=N+1 32C 33C Set XBASE to the initial vector of variables, and set the initial 34C elements of XPT, BMAT, HQ, PQ and ZMAT to zero. 35C 36 DO 20 J=1,N 37 XBASE(J)=X(J) 38 DO 10 K=1,NPT 39 10 XPT(K,J)=ZERO 40 DO 20 I=1,NDIM 41 20 BMAT(I,J)=ZERO 42 DO 30 IH=1,(N*NP)/2 43 30 HQ(IH)=ZERO 44 DO 40 K=1,NPT 45 PQ(K)=ZERO 46 DO 40 J=1,NPT-NP 47 40 ZMAT(K,J)=ZERO 48C 49C Begin the initialization procedure. NF becomes one more than the number 50C of function values so far. The coordinates of the displacement of the 51C next initial interpolation point from XBASE are set in XPT(NF+1,.). 52C 53 NF=0 54 50 NFM=NF 55 NFX=NF-N 56 NF=NF+1 57 IF (NFM .LE. 2*N) THEN 58 IF (NFM .GE. 1 .AND. NFM .LE. N) THEN 59 STEPA=RHOBEG 60 IF (SU(NFM) .EQ. ZERO) STEPA=-STEPA 61 XPT(NF,NFM)=STEPA 62 ELSE IF (NFM .GT. N) THEN 63 STEPA=XPT(NF-N,NFX) 64 STEPB=-RHOBEG 65 IF (SL(NFX) .EQ. ZERO) STEPB=DMIN1(TWO*RHOBEG,SU(NFX)) 66 IF (SU(NFX) .EQ. ZERO) STEPB=DMAX1(-TWO*RHOBEG,SL(NFX)) 67 XPT(NF,NFX)=STEPB 68 END IF 69 ELSE 70 ITEMP=(NFM-NP)/N 71 JPT=NFM-ITEMP*N-N 72 IPT=JPT+ITEMP 73 IF (IPT .GT. N) THEN 74 ITEMP=JPT 75 JPT=IPT-N 76 IPT=ITEMP 77 END IF 78 XPT(NF,IPT)=XPT(IPT+1,IPT) 79 XPT(NF,JPT)=XPT(JPT+1,JPT) 80 END IF 81C 82C Calculate the next value of F. The least function value so far and 83C its index are required. 84C 85 DO 60 J=1,N 86 X(J)=DMIN1(DMAX1(XL(J),XBASE(J)+XPT(NF,J)),XU(J)) 87 IF (XPT(NF,J) .EQ. SL(J)) X(J)=XL(J) 88 IF (XPT(NF,J) .EQ. SU(J)) X(J)=XU(J) 89 60 CONTINUE 90 F = CALFUN (N,X,IPRINT) 91c$$$ IF (IPRINT .EQ. 3) THEN 92c$$$ PRINT 70, NF,F,(X(I),I=1,N) 93c$$$ 70 FORMAT (/4X,'Function number',I6,' F =',1PD18.10, 94c$$$ 1 ' The corresponding X is:'/(2X,5D15.6)) 95c$$$ END IF 96c$$$ CALL minqi3 (IPRINT, F, NF, N, X) 97 FVAL(NF)=F 98 IF (NF .EQ. 1) THEN 99 FBEG=F 100 KOPT=1 101 ELSE IF (F .LT. FVAL(KOPT)) THEN 102 KOPT=NF 103 END IF 104C 105C Set the nonzero initial elements of BMAT and the quadratic model in the 106C cases when NF is at most 2*N+1. If NF exceeds N+1, then the positions 107C of the NF-th and (NF-N)-th interpolation points may be switched, in 108C order that the function value at the first of them contributes to the 109C off-diagonal second derivative terms of the initial quadratic model. 110C 111 IF (NF .LE. 2*N+1) THEN 112 IF (NF .GE. 2 .AND. NF .LE. N+1) THEN 113 GOPT(NFM)=(F-FBEG)/STEPA 114 IF (NPT .LT. NF+N) THEN 115 BMAT(1,NFM)=-ONE/STEPA 116 BMAT(NF,NFM)=ONE/STEPA 117 BMAT(NPT+NFM,NFM)=-HALF*RHOSQ 118 END IF 119 ELSE IF (NF .GE. N+2) THEN 120 IH=(NFX*(NFX+1))/2 121 TEMP=(F-FBEG)/STEPB 122 DIFF=STEPB-STEPA 123 HQ(IH)=TWO*(TEMP-GOPT(NFX))/DIFF 124 GOPT(NFX)=(GOPT(NFX)*STEPB-TEMP*STEPA)/DIFF 125 IF (STEPA*STEPB .LT. ZERO) THEN 126 IF (F .LT. FVAL(NF-N)) THEN 127 FVAL(NF)=FVAL(NF-N) 128 FVAL(NF-N)=F 129 IF (KOPT .EQ. NF) KOPT=NF-N 130 XPT(NF-N,NFX)=STEPB 131 XPT(NF,NFX)=STEPA 132 END IF 133 END IF 134 BMAT(1,NFX)=-(STEPA+STEPB)/(STEPA*STEPB) 135 BMAT(NF,NFX)=-HALF/XPT(NF-N,NFX) 136 BMAT(NF-N,NFX)=-BMAT(1,NFX)-BMAT(NF,NFX) 137 ZMAT(1,NFX)=DSQRT(TWO)/(STEPA*STEPB) 138 ZMAT(NF,NFX)=DSQRT(HALF)/RHOSQ 139 ZMAT(NF-N,NFX)=-ZMAT(1,NFX)-ZMAT(NF,NFX) 140 END IF 141C 142C Set the off-diagonal second derivatives of the Lagrange functions and 143C the initial quadratic model. 144C 145 ELSE 146 IH=(IPT*(IPT-1))/2+JPT 147 ZMAT(1,NFX)=RECIP 148 ZMAT(NF,NFX)=RECIP 149 ZMAT(IPT+1,NFX)=-RECIP 150 ZMAT(JPT+1,NFX)=-RECIP 151 TEMP=XPT(NF,IPT)*XPT(NF,JPT) 152 HQ(IH)=(FBEG-FVAL(IPT+1)-FVAL(JPT+1)+F)/TEMP 153 END IF 154 IF (NF .LT. NPT .AND. NF .LT. MAXFUN) GOTO 50 155 RETURN 156 END 157 158