1      SUBROUTINE NLLSQ(X,N)
2      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3      INCLUDE 'SIZES'
4      COMMON /KEYWRD/ KEYWRD
5      CHARACTER*241 KEYWRD
6      DIMENSION X(*)
7      COMMON /MESAGE/ IFLEPO,IITER
8************************************************************************
9*
10*  NLLSQ IS A NON-DERIVATIVE, NONLINEAR LEAST-SQUARES MINIMIZER. IT USES
11*        BARTEL'S PROCEDURE TO MINIMIZE A FUNCTION WHICH IS A SUM OF
12*        SQUARES.
13*
14*    ON INPUT N    = NUMBER OF UNKNOWNS
15*             X    = PARAMETERS OF FUNCTION TO BE MINIMIZED.
16*
17*    ON EXIT  X    = OPTIMIZED PARAMETERS.
18*
19*    THE FUNCTION TO BE MINIMIZED IS "COMPFG". COMPFG MUST HAVE THE
20*    CALLING SEQUENCE
21*                  CALL COMPFG(XPARAM,.TRUE.,ESCF,.TRUE.,EFS,.TRUE.)
22*                  SSQ=DOT(EFS,EFS,N)
23*    WHERE   EFS  IS A VECTOR WHICH  COMPFG  FILLS WITH THE N INDIVIDUAL
24*                 COMPONENTS OF THE ERROR FUNCTION AT THE POINT X
25*            SSQ IS THE VALUE OF THE SUM OF THE  EFS  SQUARED.
26*    IN THIS FORMULATION OF NLLSQ M AND N ARE THE SAME.
27*    THE PRECISE DEFINITIONS OF THESE TWO QUANTITIES IS:
28*
29*     N = NUMBER OF PARAMETERS TO BE OPTIMIZED.
30*     M = NUMBER OF REFERENCE FUNCTIONS. M MUST BE GREATER THEN, OR
31*         EQUAL TO, N
32************************************************************************
33C     Q = ORTHOGONAL MATRIX   (M BY M)
34C     R = RIGHT-TRIANGULAR MATRIX   (M BY N)
35C     MXCNT(1) = MAX ALLOW OVERALL FUN EVALS
36C     MXCNT(2) = MAX ALLOW NO OF FNC EVALS PER LIN SEARCH
37C     TOLS1 = RELATIVE TOLERANCE ON X OVERALL
38C     TOLS2 = ABSOLUTE TOLERANCE ON X OVERALL
39C     TOLS5 = RELATIVE TOLERANCE ON X FOR LINEAR SEARCHES
40C     TOLS6 = ABSOLUTE TOLERANCE ON X FOR LINEAR SEARCHES
41C     NRST = NUMBER OF CYCLES BETWEEN SIDESTEPS
42C     **********
43C ***** Modified by Jiro Toyoda at 1994-05-25 *****
44C     COMMON /TIME  / TIME0
45      COMMON /TIMEC / TIME0
46C ***************************** at 1994-05-25 *****
47      COMMON /NLLSQI/ NCOUNT
48      COMMON /NUMSCF/ NSCF
49      DIMENSION Y(MAXPAR), EFS(MAXPAR), P(MAXPAR)
50      COMMON /LAST  / LAST
51      COMMON /TIMDMP/ TLEFT, TDUMP
52      COMMON /NLLCOM/ Q(MAXPAR,MAXPAR),R(MAXPAR,MAXPAR*2)
53      COMMON /NLLCO2/ DDDUM(6),EFSLST(MAXPAR),XLAST(MAXPAR),IIIUM(7)
54      LOGICAL MIDDLE, RESFIL, MINPRT, LOG
55      SAVE IXSO
56      EQUIVALENCE ( IIIUM(2), ICYC),(IIIUM(3), IRST),
57     1(IIIUM(4),JRST),
58     2(DDDUM(2),ALF), (DDDUM(3),SSQ),(DDDUM(4), PN)
59      DATA IXSO/0/
60      MIDDLE=(INDEX(KEYWRD,'RESTART') .NE. 0)
61      LOG=(INDEX(KEYWRD,'NOLOG') .EQ. 0)
62      IFLEPO=10
63C*
64      M=N
65C*
66      TOL2=4.D-1
67      IF(INDEX(KEYWRD,'GNORM') .NE. 0) THEN
68         TOL2=READA(KEYWRD,INDEX(KEYWRD,'GNORM'))
69         IF(TOL2.LT.0.01D0.AND.INDEX(KEYWRD,' LET').EQ.0)THEN
70            WRITE(6,'(/,A)')'  GNORM HAS BEEN SET TOO LOW, RESET TO 0
71     1.01'
72            TOL2=0.01D0
73         ENDIF
74      ENDIF
75      LAST=0
76      TOLS1=1.D-12
77      TOLS2=1.D-10
78      TOLS5=1.D-6
79      TOLS6=1.D-3
80      NRST=4
81      TLAST=TLEFT
82      MINPRT=.TRUE.
83      RESFIL=.FALSE.
84      TLEFT=TLEFT-SECOND()+TIME0
85C     **********
86C     SET UP COUNTERS AND SWITCHES
87C     **********
88      NTO=N/6
89      IFRTL=0
90      NSST=0
91      IF(IXSO.EQ.0) IXSO=N
92      NP1 = N+1
93      NP2 = N+2
94      ICYC = 0
95      IRST = 0
96      JRST = 1
97      EPS =TOLS5
98      T = TOLS6
99C     **********
100C     GET STARTING-POINT FUNCTION VALUE
101C     SET UP ESTIMATE OF INITIAL LINE STEP
102C     **********
103      IF(MIDDLE) THEN
104         CALL PARSAV(0,N,M)
105         NSCF=IIIUM(1)
106         CLOSE(13)
107         NCOUNT=IIIUM(5)
108         DO 10 I=1,N
109   10    X(I)=XLAST(I)
110         TIME1=SECOND()
111         IF(INDEX(KEYWRD,'1SCF') .NE. 0) THEN
112            IFLEPO=13
113            LAST=1
114            RETURN
115         ENDIF
116         GOTO 60
117      ENDIF
118      CALL COMPFG(X,.TRUE.,ESCF,.TRUE.,EFSLST,.TRUE.)
119      SSQ=DOT(EFSLST,EFSLST,N)
120      NCOUNT = 1
121   20 CONTINUE
122      DO 40 I=1,M
123         DO 30 J=1,N
124            R(I,J) = 0.0D0
125            IF (I .EQ. J)  R(I,J)=1.0D0
126   30    CONTINUE
127         DO 40 J=I,M
128            Q(I,J) = 0.0D0
129            Q(J,I) = 0.0D0
130            IF (I .EQ. J)  Q(I,I)=1.0D0
131   40 CONTINUE
132      TEMP = 0.0D0
133      DO 50 I=1,N
134   50 TEMP = TEMP+X(I)**2
135      ALF = 100.0D0*(EPS*SQRT(TEMP)+T)
136C     **********
137C     MAIN LOOP
138C     **********
139      TIME1=SECOND()
140   60 CONTINUE
141C     **********
142C     UPDATE COUNTERS AND TEST FOR PRINTING THIS CYCLE
143C     **********
144      IFRTL=IFRTL+1
145      ICYC = ICYC+1
146      IRST = IRST+1
147C     **********
148C     SET  PRT,  THE LEVENBERG-MARQUARDT PARAMETER.
149C     **********
150      PRT = SQRT(SSQ)
151C     **********
152C     IF A SIDESTEP IS TO BE TAKEN, GO TO 31
153C     **********
154      IF (IRST .GE. NRST)  GO TO 210
155C     **********
156C     SOLVE THE SYSTEM    Q*R*P = -EFSLST    IN THE LEAST-SQUARES SENSE
157C     **********
158      NSST=0
159      DO 80 I=1,M
160         TEMP = 0.0D0
161         DO 70 J=1,M
162   70    TEMP = TEMP-Q(J,I)*EFSLST(J)
163   80 EFS(I) = TEMP
164      DO 90 J=1,N
165         JJ = NP1-J
166         DO 90 I=1,J
167            II = NP2-I
168   90 R(II,JJ) = R(I,J)
169      DO 160 I=1,N
170         I1 = I+1
171         Y(I) = PRT
172         EFSSS=0.0D0
173         IF (I .GE. N)  GO TO 110
174         DO 100 J=I1,N
175  100    Y(J) = 0.0D0
176  110    CONTINUE
177         DO 150 J=I,N
178            II = NP2-J
179            JJ = NP1-J
180            IF (ABS(Y(J)) .LT. ABS(R(II,JJ)))  GO TO 120
181            TEMP = Y(J)*SQRT(1.0D0+(R(II,JJ)/Y(J))**2)
182            GO TO 130
183  120       TEMP = R(II,JJ)*SQRT(1.0D0+(Y(J)/R(II,JJ))**2)
184  130       CONTINUE
185            SIN = R(II,JJ)/TEMP
186            COS = Y(J)/TEMP
187            R(II,JJ) = TEMP
188            TEMP = EFS(J)
189            EFS(J)=SIN*TEMP+COS*EFSSS
190            EFSSS=SIN*EFSSS-COS*TEMP
191            IF (J .GE. N)  GO TO 160
192            J1 = J+1
193            DO 140 K=J1,N
194               JJ = NP1-K
195               TEMP = R(II,JJ)
196               R(II,JJ) = SIN*TEMP+COS*Y(K)
197  140       Y(K) = SIN*Y(K)-COS*TEMP
198  150    CONTINUE
199  160 CONTINUE
200      P(N) = EFS(N)/R(2,1)
201      I = N
202  170 I = I-1
203      IF (I)  200,200,180
204  180 TEMP = EFS(I)
205      K = I+1
206      II = NP2-I
207      DO 190 J=K,N
208         JJ = NP1-J
209  190 TEMP = TEMP-R(II,JJ)*P(J)
210      JJ = NP1-I
211      P(I) = TEMP/R(II,JJ)
212      GO TO 170
213  200 CONTINUE
214      GO TO 230
215C     **********
216C     SIDESTEP SECTION
217C     **********
218  210 JRST = JRST+1
219      NSST=NSST+1
220      IF(NSST.GE.IXSO) GO TO 670
221      IF (JRST .GT. N)  JRST=2
222      IRST = 0
223C     **********
224C     PRODUCTION OF A VECTOR ORTHOGONAL TO THE LAST P-VECTOR
225C     **********
226      WORK = PN*(ABS(P(1))+PN)
227      TEMP = P(JRST)
228      P(1) = TEMP*(P(1)+SIGN(PN,P(1)))
229      DO 220 I=2,N
230  220 P(I) = TEMP*P(I)
231      P(JRST) = P(JRST)-WORK
232C     **********
233C     COMPUTE NORM AND NORM-SQUARE OF THE P-VECTOR
234C     **********
235  230 PNLAST = PN
236      PN=0.D0
237      PN2 = 0.0D0
238      DO 240 I=1,N
239         PN=PN+ABS(P(I))
240  240 PN2 = PN2+P(I)**2
241      IF(PN.LT.1.D-20) THEN
242         WRITE(6,'('' SYSTEM DOES NOT APPEAR TO BE OPTIMIZABLE.'',/
243     1,'' THIS CAN HAPPEN IF (A) IT WAS OPTIMIZED TO BEGIN WITH'',/
244     2,'' OR                 (B) IT IS NEITHER A GROUND NOR A'',
245     3'' TRANSITION STATE'')')
246         CALL GEOUT(1)
247         STOP
248      ENDIF
249      IF(PN2.LT.1.D-20)PN2=1.D-20
250      PN = SQRT(PN2)
251      IF(ALF.GT.1.D20)ALF=1.D20
252      IF(ICYC .GT. 1) THEN
253         ALF=ALF*1.D-20*PNLAST/PN
254         IF(ALF.GT.1.D10)        ALF=1.D10
255         ALF=ALF*1.D20
256      ENDIF
257      TTMP=ALF*PN
258      IF(TTMP.LT.0.0001D0) ALF=0.001D0/PN
259C     **********
260C     PRINTING SECTION
261C     **********
262C#      WRITE(6,501)TLEFT,ICYC,SSQ
263      DO 250 I=1,N
264         EFS(I)=X(I)
265  250 CONTINUE
266C     **********
267C     PERFORM LINE-MINIMIZATION FROM POINT X IN DIRECTION P OR -P
268C     **********
269      SSQLST = SSQ
270      DO 260 I=1,N
271         EFS(I)=0.D0
272  260 XLAST(I)=X(I)
273      CALL LOCMIN(M,X,N,P,SSQ,ALF,EFS,IERR,ESCF)
274      IF(SSQLST .LT. SSQ ) THEN
275         IF(IERR .EQ. 0)      SSQ=SSQLST
276         DO 270 I=1,N
277  270    X(I)=XLAST(I)
278         IRST=NRST
279         PN=PNLAST
280         TIME2=TIME1
281         TIME1=SECOND()
282         TCYCLE=TIME1-TIME2
283         TLEFT=TLEFT-TCYCLE
284         IF(TLEFT .GT. TCYCLE*2) GO TO 60
285         GOTO 630
286      ENDIF
287C     **********
288C     PRODUCE THE VECTOR   R*P
289C     **********
290      DO 290 I=1,N
291         TEMP = 0.0D0
292         DO 280 J=I,N
293  280    TEMP = TEMP+R(I,J)*P(J)
294  290 Y(I) = TEMP
295C     **********
296C     PRODUCE THE VECTOR ...
297C                  Y  =    (EFS-EFSLST-ALF*Q*R*P)/(ALF*(NORMSQUARE(P))
298C     COMPUTE NORM OF THIS VECTOR AS WELL
299C     **********
300      WORK = ALF*PN2
301      YN = 0.0D0
302      DO 310 I=1,M
303         TEMP = 0.0D0
304         DO 300 J=1,N
305  300    TEMP = TEMP+Q(I,J)*Y(J)
306         TEMP = (EFS(I)-EFSLST(I)-ALF*TEMP)
307         EFSLST(I) = EFS(I)
308         YN = YN+TEMP**2
309  310 EFS(I) = TEMP/WORK
310      YN = SQRT(YN)/WORK
311C     **********
312C     THE BROYDEN UPDATE   NEW MATRIX = OLD MATRIX + Y*(P-TRANS)
313C     HAS BEEN FORMED.  IT IS NOW NECESSARY TO UPDATE THE  QR DECOMP.
314C     FIRST LET    Y = (Q-TRANS)*Y.
315C     **********
316      DO 330 I=1,M
317         TEMP = 0.0D0
318         DO 320 J=1,M
319  320    TEMP = TEMP+Q(J,I)*EFS(J)
320  330 Y(I) = TEMP
321C     **********
322C     REDUCE THE VECTOR Y TO A MULTIPLE OF THE FIRST UNIT VECTOR USING
323C     A HOUSEHOLDER TRANSFORMATION FOR COMPONENTS N+1 THROUGH M AND
324C     ELEMENTARY ROTATIONS FOR THE FIRST N+1 COMPONENTS.  APPLY ALL
325C     TRANSFORMATIONS TRANSPOSED ON THE RIGHT TO THE MATRIX Q, AND
326C     APPLY THE ROTATIONS ON THE LEFT TO THE MATRIX R.
327C     THIS GIVES    (Q*(V-TRANS))*((V*R) + (V*Y)*(P-TRANS)),    WHERE
328C     V IS THE COMPOSITE OF THE TRANSFORMATIONS.  THE MATRIX
329C     ((V*R) + (V*Y)*(P-TRANS))    IS UPPER HESSENBERG.
330C     **********
331      IF (M .LE. NP1)  GO TO 390
332C
333C THE NEXT THREE LINES WERE INSERTED TO TRY TO GET ROUND OVERFLOW BUGS.
334C
335      CONST=1.D-12
336      DO 340 I=NP1,M
337  340 CONST=MAX(ABS(Y(NP1)),CONST)
338      YTAIL = 0.0D0
339      DO 350 I=NP1,M
340  350 YTAIL = YTAIL+(Y(I)/CONST)**2
341      YTAIL = SQRT(YTAIL)*CONST
342      BET = (1.0D25/YTAIL)/(YTAIL+ABS(Y(NP1)))
343      Y(NP1) = SIGN (YTAIL+ABS(Y(NP1)),Y(NP1))
344      DO 380 I=1,M
345         TMP = 0.0D0
346         DO 360 J=NP1,M
347  360    TMP = TMP+Q(I,J)*Y(J)*1.D-25
348         TMP = BET*TMP
349         DO 370 J=NP1,M
350  370    Q(I,J) = Q(I,J)-TMP*Y(J)
351  380 CONTINUE
352      Y(NP1) = YTAIL
353      I = NP1
354      GO TO 400
355  390 CONTINUE
356      I = M
357  400 CONTINUE
358  410 J = I
359      I = I-1
360      IF (I)  480,480,420
361  420 IF (Y(J))  430,410,430
362  430 IF (ABS(Y(I)) .LT. ABS(Y(J)))  GO TO 440
363      TEMP = ABS(Y(I))*SQRT(1.0D0+(Y(J)/Y(I))**2)
364      GO TO 450
365  440 TEMP = ABS(Y(J))*SQRT(1.0D0+(Y(I)/Y(J))**2)
366  450 COS = Y(I)/TEMP
367      SIN = Y(J)/TEMP
368      Y(I) = TEMP
369      DO 460 K=1,M
370         TEMP = COS*Q(K,I)+SIN*Q(K,J)
371         WORK = -SIN*Q(K,I)+COS*Q(K,J)
372         Q(K,I) = TEMP
373  460 Q(K,J) = WORK
374      IF (I .GT. N)  GO TO 410
375      R(J,I) = -SIN*R(I,I)
376      R(I,I) = COS*R(I,I)
377      IF (J .GT. N)  GO TO 410
378      DO 470 K=J,N
379         TEMP = COS*R(I,K)+SIN*R(J,K)
380         WORK = -SIN*R(I,K)+COS*R(J,K)
381         R(I,K) = TEMP
382  470 R(J,K) = WORK
383      GO TO 410
384  480 CONTINUE
385C     **********
386C     REDUCE THE UPPER-HESSENBERG MATRIX TO UPPER-TRIANGULAR FORM
387C     USING ELEMENTARY ROTATIONS.  APPLY THE SAME ROTATIONS, TRANSPOSED,
388C     ON THE RIGHT TO THE MATRIX  Q.
389C     **********
390      DO 490 K=1,N
391  490 R(1,K) = R(1,K)+YN*P(K)
392      JEND = NP1
393      IF (M .EQ. N)  JEND=N
394      DO 560 J=2,JEND
395         I = J-1
396         IF (R(J,I))  500,560,500
397  500    IF (ABS(R(I,I)) .LT. ABS(R(J,I)))  GO TO 510
398         TEMP = ABS(R(I,I))*SQRT(1.0D0+(R(J,I)/R(I,I))**2)
399         GO TO 520
400  510    TEMP = ABS(R(J,I))*SQRT(1.0D0+(R(I,I)/R(J,I))**2)
401  520    COS = R(I,I)/TEMP
402         SIN = R(J,I)/TEMP
403         R(I,I) = TEMP
404         IF (J .GT. N)  GO TO 540
405         DO 530 K=J,N
406            TEMP = COS*R(I,K)+SIN*R(J,K)
407            WORK = -SIN*R(I,K)+COS*R(J,K)
408            R(I,K) = TEMP
409  530    R(J,K) = WORK
410  540    DO 550 K=1,M
411            TEMP = COS*Q(K,I)+SIN*Q(K,J)
412            WORK = -SIN*Q(K,I)+COS*Q(K,J)
413            Q(K,I) = TEMP
414  550    Q(K,J) = WORK
415  560 CONTINUE
416C     **********
417C     CHECK THE STOPPING CRITERIA
418C     **********
419      TEMP = 0.0D0
420      DO 570 I=1,N
421  570 TEMP = TEMP+X(I)**2
422      TOLX = TOLS1*SQRT(TEMP)+TOLS2
423      IF (SQRT(ALF*PN2) .LE. TOLX)  GO TO 650
424      IF(SSQ.GE.2.D0*N) GO TO 590
425      DO 580 I=1,N
426C*****
427C     The stopping criterion is that no individual gradient be
428C         greater than TOL2
429C*****
430         IF(ABS(EFSLST(I)).GE.TOL2) GO TO 590
431  580 CONTINUE
432C#      WRITE(6,730) SSQ
433      GO TO 660
434  590 CONTINUE
435      TIME2=TIME1
436      TIME1=SECOND()
437      TCYCLE=TIME1-TIME2
438      TLEFT=TLEFT-TCYCLE
439      IF(RESFIL)THEN
440         WRITE(6,600)TLEFT,SQRT(SSQ),ESCF
441  600    FORMAT('  RESTART FILE WRITTEN,  TIME LEFT:',F9.1,
442     1' GRAD.:',F10.3,' HEAT:',G14.7)
443         RESFIL=.FALSE.
444      ELSE
445         IF(MINPRT) WRITE(6,610)ICYC,MIN(TCYCLE,9999.99D0),
446     1MIN(TLEFT,9999999.9D0),MIN(SQRT(SSQ),999999.999D0),ESCF
447         IF(LOG) WRITE(11,610)ICYC,MIN(TCYCLE,9999.99D0),
448     1MIN(TLEFT,9999999.9D0),MIN(SQRT(SSQ),999999.999D0),ESCF
449  610    FORMAT(' CYCLE:',I5,' TIME:',F6.1,' TIME LEFT:',F9.1,
450     1' GRAD.:',F10.3,' HEAT:',G14.7)
451      ENDIF
452      IF(TLAST-TLEFT.GT.TDUMP)THEN
453         TLAST=TLEFT
454         RESFIL=.TRUE.
455         DO 620 I=1,N
456  620    XLAST(I)=X(I)
457         IIIUM(1)=NSCF
458         CALL PARSAV(2,N,M)
459      ENDIF
460      IF(TLEFT .GT. TCYCLE*2) GO TO 60
461  630 IIIUM(5)=NCOUNT
462      DO 640 I=1,N
463  640 XLAST(I)=X(I)
464      IIIUM(1)=NSCF
465      CALL PARSAV(1,N,M)
466      IFLEPO=-1
467      RETURN
468  650 WRITE (6,760)  NCOUNT
469      GOTO 870
470  660 WRITE (6,770)  NCOUNT
471      GOTO 870
472  670 CONTINUE
473      WRITE(6,680) IXSO
474  680 FORMAT(1H ,5X,'ATTEMPT TO GO DOWNHILL IS UNSUCCESSFUL AFTER',I5,5X
475     1,'ORTHOGONAL SEARCHES')
476      GOTO 870
477C#  730 FORMAT(1H ,'FINAL GRADIENT =',F15.7)
478  690 FORMAT(1H ,3X,'ALF =',E12.4)
479  700 FORMAT(1H ,3X,'NCOUNT =',I5)
480  710 FORMAT(3X,'TIME LEFT:',F7.1,' CYCLE',I5,3X,'GNORM SQUARED IS'
481     1,F13.5)
482  720 FORMAT(4(5X,'X(',I2,') = ',E15.8))
483  730 FORMAT(4(5X,'P(',I2,') = ',E15.8))
484  740 FORMAT(5X,'R-MATRIX DIAGONAL ENTRIES ...')
485  750 FORMAT(6E13.3)
486  760 FORMAT('0TEST ON X SATISFIED, NUMBER OF FUNCTION CALLS = ',I5)
487  770 FORMAT('0TEST ON SSQ SATISFIED, NUMBER OF FUNCTION CALLS = ',I5)
488  780 FORMAT(' ///// NEXT CYCLE IS A SIDE-STEP ALONG THE ',I2,
489     1  '-TH NORMAL TO P')
490  790 FORMAT('0ALLOWED NUMBER OF FUNCTION CALLS EXCEEDED.'/
491     1  ' NUMBER OF FUNCTION CALLS WAS ',I5)
492  800 FORMAT('  L.-M. PARAMETER = ',E15.7,
493     1  '   SUMSQUARES CHANGE = ',E15.7)
494  810 FORMAT(1H )
495  820 FORMAT(1H )
496  830 FORMAT(1H ,3X,'I',7X,I2,9(10X,I2))
497  840 FORMAT(1H ,1X,'X(I)',1X,F10.5,2X,9(F10.5,2X))
498  850 FORMAT(1H ,1X,'G(I)',1X,F10.5,2X,9(F10.5,2X))
499  860 FORMAT(1H ,1X,'P(I)',1X,F10.5,2X,9(F10.5,2X))
500  870 LAST=1
501      RETURN
502      END
503