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