-
Notifications
You must be signed in to change notification settings - Fork 28
/
Copy pathndlfortran.f
2172 lines (2154 loc) · 73.5 KB
/
ndlfortran.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
C ----------------------------------------------------------------------
C This file contains the LBFGS algorithm and supporting routines
C
C ****************
C LBFGS SUBROUTINE
C ****************
C
SUBROUTINE LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
C
INTEGER N,M,IPRINT(2),IFLAG
DOUBLE PRECISION X(N),G(N),DIAG(N),W(N*(2*M+1)+2*M)
DOUBLE PRECISION F,EPS,XTOL
LOGICAL DIAGCO
C
C LIMITED MEMORY BFGS METHOD FOR LARGE SCALE OPTIMIZATION
C JORGE NOCEDAL
C *** July 1990 ***
C
C
C This subroutine solves the unconstrained minimization problem
C
C min F(x), x= (x1,x2,...,xN),
C
C using the limited memory BFGS method. The routine is especially
C effective on problems involving a large number of variables. In
C a typical iteration of this method an approximation Hk to the
C inverse of the Hessian is obtained by applying M BFGS updates to
C a diagonal matrix Hk0, using information from the previous M steps.
C The user specifies the number M, which determines the amount of
C storage required by the routine. The user may also provide the
C diagonal matrices Hk0 if not satisfied with the default choice.
C The algorithm is described in "On the limited memory BFGS method
C for large scale optimization", by D. Liu and J. Nocedal,
C Mathematical Programming B 45 (1989) 503-528.
C
C The user is required to calculate the function value F and its
C gradient G. In order to allow the user complete control over
C these computations, reverse communication is used. The routine
C must be called repeatedly under the control of the parameter
C IFLAG.
C
C The steplength is determined at each iteration by means of the
C line search routine MCVSRCH, which is a slight modification of
C the routine CSRCH written by More' and Thuente.
C
C The calling statement is
C
C CALL LBFGS(N,M,X,F,G,DIAGCO,DIAG,IPRINT,EPS,XTOL,W,IFLAG)
C
C where
C
C N is an INTEGER variable that must be set by the user to the
C number of variables. It is not altered by the routine.
C Restriction: N>0.
C
C M is an INTEGER variable that must be set by the user to
C the number of corrections used in the BFGS update. It
C is not altered by the routine. Values of M less than 3 are
C not recommended; large values of M will result in excessive
C computing time. 3<= M <=7 is recommended. Restriction: M>0.
C
C X is a DOUBLE PRECISION array of length N. On initial entry
C it must be set by the user to the values of the initial
C estimate of the solution vector. On exit with IFLAG=0, it
C contains the values of the variables at the best point
C found (usually a solution).
C
C F is a DOUBLE PRECISION variable. Before initial entry and on
C a re-entry with IFLAG=1, it must be set by the user to
C contain the value of the function F at the point X.
C
C G is a DOUBLE PRECISION array of length N. Before initial
C entry and on a re-entry with IFLAG=1, it must be set by
C the user to contain the components of the gradient G at
C the point X.
C
C DIAGCO is a LOGICAL variable that must be set to .TRUE. if the
C user wishes to provide the diagonal matrix Hk0 at each
C iteration. Otherwise it should be set to .FALSE., in which
C case LBFGS will use a default value described below. If
C DIAGCO is set to .TRUE. the routine will return at each
C iteration of the algorithm with IFLAG=2, and the diagonal
C matrix Hk0 must be provided in the array DIAG.
C
C
C DIAG is a DOUBLE PRECISION array of length N. If DIAGCO=.TRUE.,
C then on initial entry or on re-entry with IFLAG=2, DIAG
C it must be set by the user to contain the values of the
C diagonal matrix Hk0. Restriction: all elements of DIAG
C must be positive.
C
C IPRINT is an INTEGER array of length two which must be set by the
C user.
C
C IPRINT(1) specifies the frequency of the output:
C IPRINT(1) < 0 : no output is generated,
C IPRINT(1) = 0 : output only at first and last iteration,
C IPRINT(1) > 0 : output every IPRINT(1) iterations.
C
C IPRINT(2) specifies the type of output generated:
C IPRINT(2) = 0 : iteration count, number of function
C evaluations, function value, norm of the
C gradient, and steplength,
C IPRINT(2) = 1 : same as IPRINT(2)=0, plus vector of
C variables and gradient vector at the
C initial point,
C IPRINT(2) = 2 : same as IPRINT(2)=1, plus vector of
C variables,
C IPRINT(2) = 3 : same as IPRINT(2)=2, plus gradient vector.
C
C
C EPS is a positive DOUBLE PRECISION variable that must be set by
C the user, and determines the accuracy with which the solution
C is to be found. The subroutine terminates when
C
C ||G|| < EPS max(1,||X||),
C
C where ||.|| denotes the Euclidean norm.
C
C XTOL is a positive DOUBLE PRECISION variable that must be set by
C the user to an estimate of the machine precision (e.g.
C 10**(-16) on a SUN station 3/60). The line search routine will
C terminate if the relative width of the interval of uncertainty
C is less than XTOL.
C
C W is a DOUBLE PRECISION array of length N(2M+1)+2M used as
C workspace for LBFGS. This array must not be altered by the
C user.
C
C IFLAG is an INTEGER variable that must be set to 0 on initial entry
C to the subroutine. A return with IFLAG<0 indicates an error,
C and IFLAG=0 indicates that the routine has terminated without
C detecting errors. On a return with IFLAG=1, the user must
C evaluate the function F and gradient G. On a return with
C IFLAG=2, the user must provide the diagonal matrix Hk0.
C
C The following negative values of IFLAG, detecting an error,
C are possible:
C
C IFLAG=-1 The line search routine MCSRCH failed. The
C parameter INFO provides more detailed information
C (see also the documentation of MCSRCH):
C
C INFO = 0 IMPROPER INPUT PARAMETERS.
C
C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF
C UNCERTAINTY IS AT MOST XTOL.
C
C INFO = 3 MORE THAN 20 FUNCTION EVALUATIONS WERE
C REQUIRED AT THE PRESENT ITERATION.
C
C INFO = 4 THE STEP IS TOO SMALL.
C
C INFO = 5 THE STEP IS TOO LARGE.
C
C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
C THERE MAY NOT BE A STEP WHICH SATISFIES
C THE SUFFICIENT DECREASE AND CURVATURE
C CONDITIONS. TOLERANCES MAY BE TOO SMALL.
C
C
C IFLAG=-2 The i-th diagonal element of the diagonal inverse
C Hessian approximation, given in DIAG, is not
C positive.
C
C IFLAG=-3 Improper input parameters for LBFGS (N or M are
C not positive).
C
C
C
C ON THE DRIVER:
C
C The program that calls LBFGS must contain the declaration:
C
C EXTERNAL LB2
C
C LB2 is a BLOCK DATA that defines the default values of several
C parameters described in the COMMON section.
C
C
C
C COMMON:
C
C The subroutine contains one common area, which the user may wish to
C reference:
C
COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
C
C MP is an INTEGER variable with default value 6. It is used as the
C unit number for the printing of the monitoring information
C controlled by IPRINT.
C
C LP is an INTEGER variable with default value 6. It is used as the
C unit number for the printing of error messages. This printing
C may be suppressed by setting LP to a non-positive value.
C
C GTOL is a DOUBLE PRECISION variable with default value 0.9, which
C controls the accuracy of the line search routine MCSRCH. If the
C function and gradient evaluations are inexpensive with respect
C to the cost of the iteration (which is sometimes the case when
C solving very large problems) it may be advantageous to set GTOL
C to a small value. A typical small value is 0.1. Restriction:
C GTOL should be greater than 1.D-04.
C
C STPMIN and STPMAX are non-negative DOUBLE PRECISION variables which
C specify lower and uper bounds for the step in the line search.
C Their default values are 1.D-20 and 1.D+20, respectively. These
C values need not be modified unless the exponents are too large
C for the machine being used, or unless the problem is extremely
C badly scaled (in which case the exponents should be increased).
C
C
C MACHINE DEPENDENCIES
C
C The only variables that are machine-dependent are XTOL,
C STPMIN and STPMAX.
C
C
C GENERAL INFORMATION
C
C Other routines called directly: DAXPY, DDOT, LB1, MCSRCH
C
C Input/Output : No input; diagnostic messages on unit MP and
C error messages on unit LP.
C
C
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C
DOUBLE PRECISION GTOL,ONE,ZERO,GNORM,DDOT,STP1,FTOL,STPMIN,
. STPMAX,STP,YS,YY,SQ,YR,BETA,XNORM
INTEGER MP,LP,ITER,NFUN,POINT,ISPT,IYPT,MAXFEV,INFO,
. BOUND,NPT,CP,I,NFEV,INMC,IYCN,ISCN
LOGICAL FINISH
C
SAVE
DATA ONE,ZERO/1.0D+0,0.0D+0/
C
C INITIALIZE
C ----------
C
IF(IFLAG.EQ.0) GO TO 10
GO TO (172,100) IFLAG
10 ITER= 0
IF(N.LE.0.OR.M.LE.0) GO TO 196
IF(GTOL.LE.1.D-04) THEN
IF(LP.GT.0) WRITE(LP,245)
GTOL=9.D-01
ENDIF
NFUN= 1
POINT= 0
FINISH= .FALSE.
IF(DIAGCO) THEN
DO 30 I=1,N
30 IF (DIAG(I).LE.ZERO) GO TO 195
ELSE
DO 40 I=1,N
40 DIAG(I)= 1.0D0
ENDIF
C
C THE WORK VECTOR W IS DIVIDED AS FOLLOWS:
C ---------------------------------------
C THE FIRST N LOCATIONS ARE USED TO STORE THE GRADIENT AND
C OTHER TEMPORARY INFORMATION.
C LOCATIONS (N+1)...(N+M) STORE THE SCALARS RHO.
C LOCATIONS (N+M+1)...(N+2M) STORE THE NUMBERS ALPHA USED
C IN THE FORMULA THAT COMPUTES H*G.
C LOCATIONS (N+2M+1)...(N+2M+NM) STORE THE LAST M SEARCH
C STEPS.
C LOCATIONS (N+2M+NM+1)...(N+2M+2NM) STORE THE LAST M
C GRADIENT DIFFERENCES.
C
C THE SEARCH STEPS AND GRADIENT DIFFERENCES ARE STORED IN A
C CIRCULAR ORDER CONTROLLED BY THE PARAMETER POINT.
C
ISPT= N+2*M
IYPT= ISPT+N*M
DO 50 I=1,N
50 W(ISPT+I)= -G(I)*DIAG(I)
GNORM= DSQRT(DDOT(N,G,1,G,1))
STP1= ONE/GNORM
C
C PARAMETERS FOR LINE SEARCH ROUTINE
C
FTOL= 1.0D-4
MAXFEV= 20
C
IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
* GNORM,N,M,X,F,G,STP,FINISH)
C
C --------------------
C MAIN ITERATION LOOP
C --------------------
C
80 ITER= ITER+1
INFO=0
BOUND=ITER-1
IF(ITER.EQ.1) GO TO 165
IF (ITER .GT. M)BOUND=M
C
YS= DDOT(N,W(IYPT+NPT+1),1,W(ISPT+NPT+1),1)
IF(.NOT.DIAGCO) THEN
YY= DDOT(N,W(IYPT+NPT+1),1,W(IYPT+NPT+1),1)
DO 90 I=1,N
90 DIAG(I)= YS/YY
ELSE
IFLAG=2
RETURN
ENDIF
100 CONTINUE
IF(DIAGCO) THEN
DO 110 I=1,N
110 IF (DIAG(I).LE.ZERO) GO TO 195
ENDIF
C
C COMPUTE -H*G USING THE FORMULA GIVEN IN: Nocedal, J. 1980,
C "Updating quasi-Newton matrices with limited storage",
C Mathematics of Computation, Vol.24, No.151, pp. 773-782.
C ---------------------------------------------------------
C
CP= POINT
IF (POINT.EQ.0) CP=M
W(N+CP)= ONE/YS
DO 112 I=1,N
112 W(I)= -G(I)
CP= POINT
DO 125 I= 1,BOUND
CP=CP-1
IF (CP.EQ. -1)CP=M-1
SQ= DDOT(N,W(ISPT+CP*N+1),1,W,1)
INMC=N+M+CP+1
IYCN=IYPT+CP*N
W(INMC)= W(N+CP+1)*SQ
CALL DAXPY(N,-W(INMC),W(IYCN+1),1,W,1)
125 CONTINUE
C
DO 130 I=1,N
130 W(I)=DIAG(I)*W(I)
C
DO 145 I=1,BOUND
YR= DDOT(N,W(IYPT+CP*N+1),1,W,1)
BETA= W(N+CP+1)*YR
INMC=N+M+CP+1
BETA= W(INMC)-BETA
ISCN=ISPT+CP*N
CALL DAXPY(N,BETA,W(ISCN+1),1,W,1)
CP=CP+1
IF (CP.EQ.M)CP=0
145 CONTINUE
C
C STORE THE NEW SEARCH DIRECTION
C ------------------------------
C
DO 160 I=1,N
160 W(ISPT+POINT*N+I)= W(I)
C
C OBTAIN THE ONE-DIMENSIONAL MINIMIZER OF THE FUNCTION
C BY USING THE LINE SEARCH ROUTINE MCSRCH
C ----------------------------------------------------
165 NFEV=0
STP=ONE
IF (ITER.EQ.1) STP=STP1
DO 170 I=1,N
170 W(I)=G(I)
172 CONTINUE
CALL MCSRCH(N,X,F,G,W(ISPT+POINT*N+1),STP,FTOL,
* XTOL,MAXFEV,INFO,NFEV,DIAG)
IF (INFO .EQ. -1) THEN
IFLAG=1
RETURN
ENDIF
IF (INFO .NE. 1) GO TO 190
NFUN= NFUN + NFEV
C
C COMPUTE THE NEW STEP AND GRADIENT CHANGE
C -----------------------------------------
C
NPT=POINT*N
DO 175 I=1,N
W(ISPT+NPT+I)= STP*W(ISPT+NPT+I)
175 W(IYPT+NPT+I)= G(I)-W(I)
POINT=POINT+1
IF (POINT.EQ.M)POINT=0
C
C TERMINATION TEST
C ----------------
C
GNORM= DSQRT(DDOT(N,G,1,G,1))
XNORM= DSQRT(DDOT(N,X,1,X,1))
XNORM= DMAX1(1.0D0,XNORM)
IF (GNORM/XNORM .LE. EPS) FINISH=.TRUE.
C
IF(IPRINT(1).GE.0) CALL LB1(IPRINT,ITER,NFUN,
* GNORM,N,M,X,F,G,STP,FINISH)
IF (FINISH) THEN
IFLAG=0
RETURN
ENDIF
GO TO 80
C
C ------------------------------------------------------------
C END OF MAIN ITERATION LOOP. ERROR EXITS.
C ------------------------------------------------------------
C
190 IFLAG=-1
IF(LP.GT.0) WRITE(LP,200) INFO
RETURN
195 IFLAG=-2
IF(LP.GT.0) WRITE(LP,235) I
RETURN
196 IFLAG= -3
IF(LP.GT.0) WRITE(LP,240)
C
C FORMATS
C -------
C
200 FORMAT(/' IFLAG= -1 ',/' LINE SEARCH FAILED. SEE',
. ' DOCUMENTATION OF ROUTINE MCSRCH',/' ERROR RETURN',
. ' OF LINE SEARCH: INFO= ',I2,/
. ' POSSIBLE CAUSES: FUNCTION OR GRADIENT ARE INCORRECT',/,
. ' OR INCORRECT TOLERANCES')
235 FORMAT(/' IFLAG= -2',/' THE',I5,'-TH DIAGONAL ELEMENT OF THE',/,
. ' INVERSE HESSIAN APPROXIMATION IS NOT POSITIVE')
240 FORMAT(/' IFLAG= -3',/' IMPROPER INPUT PARAMETERS (N OR M',
. ' ARE NOT POSITIVE)')
245 FORMAT(/' GTOL IS LESS THAN OR EQUAL TO 1.D-04',
. / ' IT HAS BEEN RESET TO 9.D-01')
RETURN
END
C
C LAST LINE OF SUBROUTINE LBFGS
C
C
SUBROUTINE LB1(IPRINT,ITER,NFUN,
* GNORM,N,M,X,F,G,STP,FINISH)
C
C -------------------------------------------------------------
C THIS ROUTINE PRINTS MONITORING INFORMATION. THE FREQUENCY AND
C AMOUNT OF OUTPUT ARE CONTROLLED BY IPRINT.
C -------------------------------------------------------------
C
INTEGER IPRINT(2),ITER,NFUN,LP,MP,N,M
DOUBLE PRECISION X(N),G(N),F,GNORM,STP,GTOL,STPMIN,STPMAX
LOGICAL FINISH
COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
C
IF (ITER.EQ.0)THEN
WRITE(MP,10)
WRITE(MP,20) N,M
WRITE(MP,30)F,GNORM
IF (IPRINT(2).GE.1)THEN
WRITE(MP,40)
WRITE(MP,50) (X(I),I=1,N)
WRITE(MP,60)
WRITE(MP,50) (G(I),I=1,N)
ENDIF
WRITE(MP,10)
WRITE(MP,70)
ELSE
IF ((IPRINT(1).EQ.0).AND.(ITER.NE.1.AND..NOT.FINISH))RETURN
IF (IPRINT(1).NE.0)THEN
IF(MOD(ITER-1,IPRINT(1)).EQ.0.OR.FINISH)THEN
IF(IPRINT(2).GT.1.AND.ITER.GT.1) WRITE(MP,70)
WRITE(MP,80)ITER,NFUN,F,GNORM,STP
ELSE
RETURN
ENDIF
ELSE
IF( IPRINT(2).GT.1.AND.FINISH) WRITE(MP,70)
WRITE(MP,80)ITER,NFUN,F,GNORM,STP
ENDIF
IF (IPRINT(2).EQ.2.OR.IPRINT(2).EQ.3)THEN
IF (FINISH)THEN
WRITE(MP,90)
ELSE
WRITE(MP,40)
ENDIF
WRITE(MP,50)(X(I),I=1,N)
IF (IPRINT(2).EQ.3)THEN
WRITE(MP,60)
WRITE(MP,50)(G(I),I=1,N)
ENDIF
ENDIF
IF (FINISH) WRITE(MP,100)
ENDIF
C
10 FORMAT('*************************************************')
20 FORMAT(' N=',I5,' NUMBER OF CORRECTIONS=',I2,
. /, ' INITIAL VALUES')
30 FORMAT(' F= ',1PD10.3,' GNORM= ',1PD10.3)
40 FORMAT(' VECTOR X= ')
50 FORMAT(6(2X,1PD10.3))
60 FORMAT(' GRADIENT VECTOR G= ')
70 FORMAT(/' I NFN',4X,'FUNC',8X,'GNORM',7X,'STEPLENGTH'/)
80 FORMAT(2(I4,1X),3X,3(1PD10.3,2X))
90 FORMAT(' FINAL POINT X= ')
100 FORMAT(/' THE MINIMIZATION TERMINATED WITHOUT DETECTING ERRORS.',
. /' IFLAG = 0')
C
RETURN
END
C ******
C
C
C ----------------------------------------------------------
C DATA
C ----------------------------------------------------------
C
BLOCK DATA LB2
INTEGER LP,MP
DOUBLE PRECISION GTOL,STPMIN,STPMAX
COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
DATA MP,LP,GTOL,STPMIN,STPMAX/6,6,9.0D-01,1.0D-20,1.0D+20/
END
C
C
C ----------------------------------------------------------
C
subroutine daxpy(n,da,dx,incx,dy,incy)
c
c constant times a vector plus a vector.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
double precision dx(1),dy(1),da
integer i,incx,incy,ix,iy,m,mp1,n
c
if(n.le.0)return
if (da .eq. 0.0d0) return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dy(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
10 continue
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,4)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dy(i) = dy(i) + da*dx(i)
30 continue
if( n .lt. 4 ) return
40 mp1 = m + 1
do 50 i = mp1,n,4
dy(i) = dy(i) + da*dx(i)
dy(i + 1) = dy(i + 1) + da*dx(i + 1)
dy(i + 2) = dy(i + 2) + da*dx(i + 2)
dy(i + 3) = dy(i + 3) + da*dx(i + 3)
50 continue
return
end
C
C
C ----------------------------------------------------------
C
double precision function ddot(n,dx,incx,dy,incy)
c
c forms the dot product of two vectors.
c uses unrolled loops for increments equal to one.
c jack dongarra, linpack, 3/11/78.
c
double precision dx(1),dy(1),dtemp
integer i,incx,incy,ix,iy,m,mp1,n
c
ddot = 0.0d0
dtemp = 0.0d0
if(n.le.0)return
if(incx.eq.1.and.incy.eq.1)go to 20
c
c code for unequal increments or equal increments
c not equal to 1
c
ix = 1
iy = 1
if(incx.lt.0)ix = (-n+1)*incx + 1
if(incy.lt.0)iy = (-n+1)*incy + 1
do 10 i = 1,n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
10 continue
ddot = dtemp
return
c
c code for both increments equal to 1
c
c
c clean-up loop
c
20 m = mod(n,5)
if( m .eq. 0 ) go to 40
do 30 i = 1,m
dtemp = dtemp + dx(i)*dy(i)
30 continue
if( n .lt. 5 ) go to 60
40 mp1 = m + 1
do 50 i = mp1,n,5
dtemp = dtemp + dx(i)*dy(i) + dx(i + 1)*dy(i + 1) +
* dx(i + 2)*dy(i + 2) + dx(i + 3)*dy(i + 3) + dx(i + 4)*dy(i + 4)
50 continue
60 ddot = dtemp
return
end
C ------------------------------------------------------------------
C
C **************************
C LINE SEARCH ROUTINE MCSRCH
C **************************
C
SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL,MAXFEV,INFO,NFEV,WA)
INTEGER N,MAXFEV,INFO,NFEV
DOUBLE PRECISION F,STP,FTOL,GTOL,XTOL,STPMIN,STPMAX
DOUBLE PRECISION X(N),G(N),S(N),WA(N)
COMMON /LB3/MP,LP,GTOL,STPMIN,STPMAX
SAVE
C
C SUBROUTINE MCSRCH
C
C A slight modification of the subroutine CSRCH of More' and Thuente.
C The changes are to allow reverse communication, and do not affect
C the performance of the routine.
C
C THE PURPOSE OF MCSRCH IS TO FIND A STEP WHICH SATISFIES
C A SUFFICIENT DECREASE CONDITION AND A CURVATURE CONDITION.
C
C AT EACH STAGE THE SUBROUTINE UPDATES AN INTERVAL OF
C UNCERTAINTY WITH ENDPOINTS STX AND STY. THE INTERVAL OF
C UNCERTAINTY IS INITIALLY CHOSEN SO THAT IT CONTAINS A
C MINIMIZER OF THE MODIFIED FUNCTION
C
C F(X+STP*S) - F(X) - FTOL*STP*(GRADF(X)'S).
C
C IF A STEP IS OBTAINED FOR WHICH THE MODIFIED FUNCTION
C HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE DERIVATIVE,
C THEN THE INTERVAL OF UNCERTAINTY IS CHOSEN SO THAT IT
C CONTAINS A MINIMIZER OF F(X+STP*S).
C
C THE ALGORITHM IS DESIGNED TO FIND A STEP WHICH SATISFIES
C THE SUFFICIENT DECREASE CONDITION
C
C F(X+STP*S) .LE. F(X) + FTOL*STP*(GRADF(X)'S),
C
C AND THE CURVATURE CONDITION
C
C ABS(GRADF(X+STP*S)'S)) .LE. GTOL*ABS(GRADF(X)'S).
C
C IF FTOL IS LESS THAN GTOL AND IF, FOR EXAMPLE, THE FUNCTION
C IS BOUNDED BELOW, THEN THERE IS ALWAYS A STEP WHICH SATISFIES
C BOTH CONDITIONS. IF NO STEP CAN BE FOUND WHICH SATISFIES BOTH
C CONDITIONS, THEN THE ALGORITHM USUALLY STOPS WHEN ROUNDING
C ERRORS PREVENT FURTHER PROGRESS. IN THIS CASE STP ONLY
C SATISFIES THE SUFFICIENT DECREASE CONDITION.
C
C THE SUBROUTINE STATEMENT IS
C
C SUBROUTINE MCSRCH(N,X,F,G,S,STP,FTOL,XTOL, MAXFEV,INFO,NFEV,WA)
C WHERE
C
C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C OF VARIABLES.
C
C X IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
C BASE POINT FOR THE LINE SEARCH. ON OUTPUT IT CONTAINS
C X + STP*S.
C
C F IS A VARIABLE. ON INPUT IT MUST CONTAIN THE VALUE OF F
C AT X. ON OUTPUT IT CONTAINS THE VALUE OF F AT X + STP*S.
C
C G IS AN ARRAY OF LENGTH N. ON INPUT IT MUST CONTAIN THE
C GRADIENT OF F AT X. ON OUTPUT IT CONTAINS THE GRADIENT
C OF F AT X + STP*S.
C
C S IS AN INPUT ARRAY OF LENGTH N WHICH SPECIFIES THE
C SEARCH DIRECTION.
C
C STP IS A NONNEGATIVE VARIABLE. ON INPUT STP CONTAINS AN
C INITIAL ESTIMATE OF A SATISFACTORY STEP. ON OUTPUT
C STP CONTAINS THE FINAL ESTIMATE.
C
C FTOL AND GTOL ARE NONNEGATIVE INPUT VARIABLES. (In this reverse
C communication implementation GTOL is defined in a COMMON
C statement.) TERMINATION OCCURS WHEN THE SUFFICIENT DECREASE
C CONDITION AND THE DIRECTIONAL DERIVATIVE CONDITION ARE
C SATISFIED.
C
C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION OCCURS
C WHEN THE RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C IS AT MOST XTOL.
C
C STPMIN AND STPMAX ARE NONNEGATIVE INPUT VARIABLES WHICH
C SPECIFY LOWER AND UPPER BOUNDS FOR THE STEP. (In this reverse
C communication implementatin they are defined in a COMMON
C statement).
C
C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
C OCCURS WHEN THE NUMBER OF CALLS TO FCN IS AT LEAST
C MAXFEV BY THE END OF AN ITERATION.
C
C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C
C INFO = 0 IMPROPER INPUT PARAMETERS.
C
C INFO =-1 A RETURN IS MADE TO COMPUTE THE FUNCTION AND GRADIENT.
C
C INFO = 1 THE SUFFICIENT DECREASE CONDITION AND THE
C DIRECTIONAL DERIVATIVE CONDITION HOLD.
C
C INFO = 2 RELATIVE WIDTH OF THE INTERVAL OF UNCERTAINTY
C IS AT MOST XTOL.
C
C INFO = 3 NUMBER OF CALLS TO FCN HAS REACHED MAXFEV.
C
C INFO = 4 THE STEP IS AT THE LOWER BOUND STPMIN.
C
C INFO = 5 THE STEP IS AT THE UPPER BOUND STPMAX.
C
C INFO = 6 ROUNDING ERRORS PREVENT FURTHER PROGRESS.
C THERE MAY NOT BE A STEP WHICH SATISFIES THE
C SUFFICIENT DECREASE AND CURVATURE CONDITIONS.
C TOLERANCES MAY BE TOO SMALL.
C
C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
C CALLS TO FCN.
C
C WA IS A WORK ARRAY OF LENGTH N.
C
C SUBPROGRAMS CALLED
C
C MCSTEP
C
C FORTRAN-SUPPLIED...ABS,MAX,MIN
C
C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
C JORGE J. MORE', DAVID J. THUENTE
C
C **********
INTEGER INFOC,J
LOGICAL BRACKT,STAGE1
DOUBLE PRECISION DG,DGM,DGINIT,DGTEST,DGX,DGXM,DGY,DGYM,
* FINIT,FTEST1,FM,FX,FXM,FY,FYM,P5,P66,STX,STY,
* STMIN,STMAX,WIDTH,WIDTH1,XTRAPF,ZERO
DATA P5,P66,XTRAPF,ZERO /0.5D0,0.66D0,4.0D0,0.0D0/
IF(INFO.EQ.-1) GO TO 45
INFOC = 1
C
C CHECK THE INPUT PARAMETERS FOR ERRORS.
C
IF (N .LE. 0 .OR. STP .LE. ZERO .OR. FTOL .LT. ZERO .OR.
* GTOL .LT. ZERO .OR. XTOL .LT. ZERO .OR. STPMIN .LT. ZERO
* .OR. STPMAX .LT. STPMIN .OR. MAXFEV .LE. 0) RETURN
C
C COMPUTE THE INITIAL GRADIENT IN THE SEARCH DIRECTION
C AND CHECK THAT S IS A DESCENT DIRECTION.
C
DGINIT = ZERO
DO 10 J = 1, N
DGINIT = DGINIT + G(J)*S(J)
10 CONTINUE
IF (DGINIT .GE. ZERO) then
write(LP,15)
15 FORMAT(/' THE SEARCH DIRECTION IS NOT A DESCENT DIRECTION')
RETURN
ENDIF
C
C INITIALIZE LOCAL VARIABLES.
C
BRACKT = .FALSE.
STAGE1 = .TRUE.
NFEV = 0
FINIT = F
DGTEST = FTOL*DGINIT
WIDTH = STPMAX - STPMIN
WIDTH1 = WIDTH/P5
DO 20 J = 1, N
WA(J) = X(J)
20 CONTINUE
C
C THE VARIABLES STX, FX, DGX CONTAIN THE VALUES OF THE STEP,
C FUNCTION, AND DIRECTIONAL DERIVATIVE AT THE BEST STEP.
C THE VARIABLES STY, FY, DGY CONTAIN THE VALUE OF THE STEP,
C FUNCTION, AND DERIVATIVE AT THE OTHER ENDPOINT OF
C THE INTERVAL OF UNCERTAINTY.
C THE VARIABLES STP, F, DG CONTAIN THE VALUES OF THE STEP,
C FUNCTION, AND DERIVATIVE AT THE CURRENT STEP.
C
STX = ZERO
FX = FINIT
DGX = DGINIT
STY = ZERO
FY = FINIT
DGY = DGINIT
C
C START OF ITERATION.
C
30 CONTINUE
C
C SET THE MINIMUM AND MAXIMUM STEPS TO CORRESPOND
C TO THE PRESENT INTERVAL OF UNCERTAINTY.
C
IF (BRACKT) THEN
STMIN = MIN(STX,STY)
STMAX = MAX(STX,STY)
ELSE
STMIN = STX
STMAX = STP + XTRAPF*(STP - STX)
END IF
C
C FORCE THE STEP TO BE WITHIN THE BOUNDS STPMAX AND STPMIN.
C
STP = MAX(STP,STPMIN)
STP = MIN(STP,STPMAX)
C
C IF AN UNUSUAL TERMINATION IS TO OCCUR THEN LET
C STP BE THE LOWEST POINT OBTAINED SO FAR.
C
IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
* .OR. NFEV .GE. MAXFEV-1 .OR. INFOC .EQ. 0
* .OR. (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX)) STP = STX
C
C EVALUATE THE FUNCTION AND GRADIENT AT STP
C AND COMPUTE THE DIRECTIONAL DERIVATIVE.
C We return to main program to obtain F and G.
C
DO 40 J = 1, N
X(J) = WA(J) + STP*S(J)
40 CONTINUE
INFO=-1
RETURN
C
45 INFO=0
NFEV = NFEV + 1
DG = ZERO
DO 50 J = 1, N
DG = DG + G(J)*S(J)
50 CONTINUE
FTEST1 = FINIT + STP*DGTEST
C
C TEST FOR CONVERGENCE.
C
IF ((BRACKT .AND. (STP .LE. STMIN .OR. STP .GE. STMAX))
* .OR. INFOC .EQ. 0) INFO = 6
IF (STP .EQ. STPMAX .AND.
* F .LE. FTEST1 .AND. DG .LE. DGTEST) INFO = 5
IF (STP .EQ. STPMIN .AND.
* (F .GT. FTEST1 .OR. DG .GE. DGTEST)) INFO = 4
IF (NFEV .GE. MAXFEV) INFO = 3
IF (BRACKT .AND. STMAX-STMIN .LE. XTOL*STMAX) INFO = 2
IF (F .LE. FTEST1 .AND. ABS(DG) .LE. GTOL*(-DGINIT)) INFO = 1
C
C CHECK FOR TERMINATION.
C
IF (INFO .NE. 0) RETURN
C
C IN THE FIRST STAGE WE SEEK A STEP FOR WHICH THE MODIFIED
C FUNCTION HAS A NONPOSITIVE VALUE AND NONNEGATIVE DERIVATIVE.
C
IF (STAGE1 .AND. F .LE. FTEST1 .AND.
* DG .GE. MIN(FTOL,GTOL)*DGINIT) STAGE1 = .FALSE.
C
C A MODIFIED FUNCTION IS USED TO PREDICT THE STEP ONLY IF
C WE HAVE NOT OBTAINED A STEP FOR WHICH THE MODIFIED
C FUNCTION HAS A NONPOSITIVE FUNCTION VALUE AND NONNEGATIVE
C DERIVATIVE, AND IF A LOWER FUNCTION VALUE HAS BEEN
C OBTAINED BUT THE DECREASE IS NOT SUFFICIENT.
C
IF (STAGE1 .AND. F .LE. FX .AND. F .GT. FTEST1) THEN
C
C DEFINE THE MODIFIED FUNCTION AND DERIVATIVE VALUES.
C
FM = F - STP*DGTEST
FXM = FX - STX*DGTEST
FYM = FY - STY*DGTEST
DGM = DG - DGTEST
DGXM = DGX - DGTEST
DGYM = DGY - DGTEST
C
C CALL CSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
C AND TO COMPUTE THE NEW STEP.
C
CALL MCSTEP(STX,FXM,DGXM,STY,FYM,DGYM,STP,FM,DGM,
* BRACKT,STMIN,STMAX,INFOC)
C
C RESET THE FUNCTION AND GRADIENT VALUES FOR F.
C
FX = FXM + STX*DGTEST
FY = FYM + STY*DGTEST
DGX = DGXM + DGTEST
DGY = DGYM + DGTEST
ELSE
C
C CALL MCSTEP TO UPDATE THE INTERVAL OF UNCERTAINTY
C AND TO COMPUTE THE NEW STEP.
C
CALL MCSTEP(STX,FX,DGX,STY,FY,DGY,STP,F,DG,
* BRACKT,STMIN,STMAX,INFOC)
END IF
C
C FORCE A SUFFICIENT DECREASE IN THE SIZE OF THE
C INTERVAL OF UNCERTAINTY.
C
IF (BRACKT) THEN
IF (ABS(STY-STX) .GE. P66*WIDTH1)
* STP = STX + P5*(STY - STX)
WIDTH1 = WIDTH
WIDTH = ABS(STY-STX)
END IF
C
C END OF ITERATION.
C
GO TO 30
C
C LAST LINE OF SUBROUTINE MCSRCH.
C
END
SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
* STPMIN,STPMAX,INFO)
INTEGER INFO
DOUBLE PRECISION STX,FX,DX,STY,FY,DY,STP,FP,DP,STPMIN,STPMAX
LOGICAL BRACKT,BOUND
C
C SUBROUTINE MCSTEP
C
C THE PURPOSE OF MCSTEP IS TO COMPUTE A SAFEGUARDED STEP FOR
C A LINESEARCH AND TO UPDATE AN INTERVAL OF UNCERTAINTY FOR
C A MINIMIZER OF THE FUNCTION.
C
C THE PARAMETER STX CONTAINS THE STEP WITH THE LEAST FUNCTION
C VALUE. THE PARAMETER STP CONTAINS THE CURRENT STEP. IT IS
C ASSUMED THAT THE DERIVATIVE AT STX IS NEGATIVE IN THE
C DIRECTION OF THE STEP. IF BRACKT IS SET TRUE THEN A
C MINIMIZER HAS BEEN BRACKETED IN AN INTERVAL OF UNCERTAINTY
C WITH ENDPOINTS STX AND STY.
C
C THE SUBROUTINE STATEMENT IS
C
C SUBROUTINE MCSTEP(STX,FX,DX,STY,FY,DY,STP,FP,DP,BRACKT,
C STPMIN,STPMAX,INFO)
C
C WHERE
C
C STX, FX, AND DX ARE VARIABLES WHICH SPECIFY THE STEP,
C THE FUNCTION, AND THE DERIVATIVE AT THE BEST STEP OBTAINED
C SO FAR. THE DERIVATIVE MUST BE NEGATIVE IN THE DIRECTION
C OF THE STEP, THAT IS, DX AND STP-STX MUST HAVE OPPOSITE
C SIGNS. ON OUTPUT THESE PARAMETERS ARE UPDATED APPROPRIATELY.
C
C STY, FY, AND DY ARE VARIABLES WHICH SPECIFY THE STEP,
C THE FUNCTION, AND THE DERIVATIVE AT THE OTHER ENDPOINT OF
C THE INTERVAL OF UNCERTAINTY. ON OUTPUT THESE PARAMETERS ARE
C UPDATED APPROPRIATELY.
C
C STP, FP, AND DP ARE VARIABLES WHICH SPECIFY THE STEP,
C THE FUNCTION, AND THE DERIVATIVE AT THE CURRENT STEP.
C IF BRACKT IS SET TRUE THEN ON INPUT STP MUST BE
C BETWEEN STX AND STY. ON OUTPUT STP IS SET TO THE NEW STEP.
C
C BRACKT IS A LOGICAL VARIABLE WHICH SPECIFIES IF A MINIMIZER
C HAS BEEN BRACKETED. IF THE MINIMIZER HAS NOT BEEN BRACKETED
C THEN ON INPUT BRACKT MUST BE SET FALSE. IF THE MINIMIZER
C IS BRACKETED THEN ON OUTPUT BRACKT IS SET TRUE.
C
C STPMIN AND STPMAX ARE INPUT VARIABLES WHICH SPECIFY LOWER
C AND UPPER BOUNDS FOR THE STEP.
C
C INFO IS AN INTEGER OUTPUT VARIABLE SET AS FOLLOWS:
C IF INFO = 1,2,3,4,5, THEN THE STEP HAS BEEN COMPUTED
C ACCORDING TO ONE OF THE FIVE CASES BELOW. OTHERWISE
C INFO = 0, AND THIS INDICATES IMPROPER INPUT PARAMETERS.
C
C SUBPROGRAMS CALLED
C
C FORTRAN-SUPPLIED ... ABS,MAX,MIN,SQRT
C
C ARGONNE NATIONAL LABORATORY. MINPACK PROJECT. JUNE 1983
C JORGE J. MORE', DAVID J. THUENTE
C
DOUBLE PRECISION GAMMA,P,Q,R,S,SGND,STPC,STPF,STPQ,THETA
INFO = 0
C
C CHECK THE INPUT PARAMETERS FOR ERRORS.
C
IF ((BRACKT .AND. (STP .LE. MIN(STX,STY) .OR.
* STP .GE. MAX(STX,STY))) .OR.
* DX*(STP-STX) .GE. 0.0 .OR. STPMAX .LT. STPMIN) RETURN
C
C DETERMINE IF THE DERIVATIVES HAVE OPPOSITE SIGN.
C