forked from compxco/genray
-
Notifications
You must be signed in to change notification settings - Fork 0
/
dmnf.f
4863 lines (4837 loc) · 164 KB
/
dmnf.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
SUBROUTINE DMNF(N, D, X, CALCF, IV, LIV, LV, V,
1 UIPARM, URPARM, UFPARM)
C
C *** MINIMIZE GENERAL UNCONSTRAINED OBJECTIVE FUNCTION USING
C *** FINITE-DIFFERENCE GRADIENTS AND SECANT HESSIAN APPROXIMATIONS.
C
INTEGER N, LIV, LV
INTEGER IV(LIV), UIPARM(1)
DOUBLE PRECISION D(N), X(N), V(LV), URPARM(1)
C DIMENSION V(77 + N*(N+17)/2), UIPARM(*), URPARM(*)
EXTERNAL CALCF, UFPARM
C
C *** PURPOSE ***
C
C THIS ROUTINE INTERACTS WITH SUBROUTINE DRMNF IN AN ATTEMPT
C TO FIND AN N-VECTOR X* THAT MINIMIZES THE (UNCONSTRAINED)
C OBJECTIVE FUNCTION COMPUTED BY CALCF. (OFTEN THE X* FOUND IS
C A LOCAL MINIMIZER RATHER THAN A GLOBAL ONE.)
C
C *** PARAMETERS ***
C
C THE PARAMETERS FOR DMNF ARE THE SAME AS THOSE FOR DMNG
C (WHICH SEE), EXCEPT THAT CALCG IS OMITTED. INSTEAD OF CALLING
C CALCG TO OBTAIN THE GRADIENT OF THE OBJECTIVE FUNCTION AT X,
C DMNF CALLS DS7GRD, WHICH COMPUTES AN APPROXIMATION TO THE
C GRADIENT BY FINITE (FORWARD AND CENTRAL) DIFFERENCES USING THE
C METHOD OF REF. 1. THE FOLLOWING INPUT COMPONENT IS OF INTEREST
C IN THIS REGARD (AND IS NOT DESCRIBED IN DMNG).
C
C V(ETA0)..... V(42) IS AN ESTIMATED BOUND ON THE RELATIVE ERROR IN THE
C OBJECTIVE FUNCTION VALUE COMPUTED BY CALCF...
C (TRUE VALUE) = (COMPUTED VALUE) * (1 + E),
C WHERE ABS(E) .LE. V(ETA0). DEFAULT = MACHEP * 10**3,
C WHERE MACHEP IS THE UNIT ROUNDOFF.
C
C THE OUTPUT VALUES IV(NFCALL) AND IV(NGCALL) HAVE DIFFERENT
C MEANINGS FOR DMNF THAN FOR DMNG...
C
C IV(NFCALL)... IV(6) IS THE NUMBER OF CALLS SO FAR MADE ON CALCF (I.E.,
C FUNCTION EVALUATIONS) EXCLUDING THOSE MADE ONLY FOR
C COMPUTING GRADIENTS. THE INPUT VALUE IV(MXFCAL) IS A
C LIMIT ON IV(NFCALL).
C IV(NGCALL)... IV(30) IS THE NUMBER OF FUNCTION EVALUATIONS MADE ONLY
C FOR COMPUTING GRADIENTS. THE TOTAL NUMBER OF FUNCTION
C EVALUATIONS IS THUS IV(NFCALL) + IV(NGCALL).
C
C *** REFERENCE ***
C
C 1. STEWART, G.W. (1967), A MODIFICATION OF DAVIDON*S MINIMIZATION
C METHOD TO ACCEPT DIFFERENCE APPROXIMATIONS OF DERIVATIVES,
C J. ASSOC. COMPUT. MACH. 14, PP. 72-83.
C.
C *** GENERAL ***
C
C CODED BY DAVID M. GAY (WINTER 1980). REVISED SEPT. 1982.
C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C SUPPORTED IN PART BY THE NATIONAL SCIENCE FOUNDATION UNDER
C GRANTS MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989,
C AND MCS-7906671.
C
C
C---------------------------- DECLARATIONS ---------------------------
C
EXTERNAL DRMNF
C
C DRMNF.... OVERSEES COMPUTATION OF FINITE-DIFFERENCE GRADIENT AND
C CALLS DRMNG TO CARRY OUT DMNG ALGORITHM.
C
INTEGER NF
DOUBLE PRECISION FX
C
C *** SUBSCRIPTS FOR IV ***
C
INTEGER NFCALL, TOOBIG
C
C/6
C DATA NFCALL/6/, TOOBIG/2/
C/7
PARAMETER (NFCALL=6, TOOBIG=2)
C/
C
C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
C
write(*,*)'dmnf before drmnf'
10 CALL DRMNF(D, FX, IV, LIV, LV, N, V, X)
write(*,*)'dmnf after drmnf'
IF (IV(1) .GT. 2) GO TO 999
C
C *** COMPUTE FUNCTION ***
C
NF = IV(NFCALL)
CALL CALCF(N, X, NF, FX, UIPARM, URPARM, UFPARM)
IF (NF .LE. 0) IV(TOOBIG) = 1
GO TO 10
C
C
999 RETURN
C *** LAST CARD OF DMNF FOLLOWS ***
END
SUBROUTINE I1MCRY(A, A1, B, C, D)
**** SPECIAL COMPUTATION FOR OLD CRAY MACHINES ****
INTEGER A, A1, B, C, D
A1 = 16777216*B + C
A = 16777216*A1 + D
END
SUBROUTINE DA7SST(IV, LIV, LV, V)
C
C *** ASSESS CANDIDATE STEP (***SOL VERSION 2.3) ***
C
INTEGER LIV, LV
INTEGER IV(LIV)
DOUBLE PRECISION V(LV)
C
C *** PURPOSE ***
C
C THIS SUBROUTINE IS CALLED BY AN UNCONSTRAINED MINIMIZATION
C ROUTINE TO ASSESS THE NEXT CANDIDATE STEP. IT MAY RECOMMEND ONE
C OF SEVERAL COURSES OF ACTION, SUCH AS ACCEPTING THE STEP, RECOM-
C PUTING IT USING THE SAME OR A NEW QUADRATIC MODEL, OR HALTING DUE
C TO CONVERGENCE OR FALSE CONVERGENCE. SEE THE RETURN CODE LISTING
C BELOW.
C
C-------------------------- PARAMETER USAGE --------------------------
C
C IV (I/O) INTEGER PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION
C BELOW OF IV VALUES REFERENCED.
C LIV (IN) LENGTH OF IV ARRAY.
C LV (IN) LENGTH OF V ARRAY.
C V (I/O) REAL PARAMETER AND SCRATCH VECTOR -- SEE DESCRIPTION
C BELOW OF V VALUES REFERENCED.
C
C *** IV VALUES REFERENCED ***
C
C IV(IRC) (I/O) ON INPUT FOR THE FIRST STEP TRIED IN A NEW ITERATION,
C IV(IRC) SHOULD BE SET TO 3 OR 4 (THE VALUE TO WHICH IT IS
C SET WHEN STEP IS DEFINITELY TO BE ACCEPTED). ON INPUT
C AFTER STEP HAS BEEN RECOMPUTED, IV(IRC) SHOULD BE
C UNCHANGED SINCE THE PREVIOUS RETURN OF DA7SST.
C ON OUTPUT, IV(IRC) IS A RETURN CODE HAVING ONE OF THE
C FOLLOWING VALUES...
C 1 = SWITCH MODELS OR TRY SMALLER STEP.
C 2 = SWITCH MODELS OR ACCEPT STEP.
C 3 = ACCEPT STEP AND DETERMINE V(RADFAC) BY GRADIENT
C TESTS.
C 4 = ACCEPT STEP, V(RADFAC) HAS BEEN DETERMINED.
C 5 = RECOMPUTE STEP (USING THE SAME MODEL).
C 6 = RECOMPUTE STEP WITH RADIUS = V(LMAXS) BUT DO NOT
C EVALUATE THE OBJECTIVE FUNCTION.
C 7 = X-CONVERGENCE (SEE V(XCTOL)).
C 8 = RELATIVE FUNCTION CONVERGENCE (SEE V(RFCTOL)).
C 9 = BOTH X- AND RELATIVE FUNCTION CONVERGENCE.
C 10 = ABSOLUTE FUNCTION CONVERGENCE (SEE V(AFCTOL)).
C 11 = SINGULAR CONVERGENCE (SEE V(LMAXS)).
C 12 = FALSE CONVERGENCE (SEE V(XFTOL)).
C 13 = IV(IRC) WAS OUT OF RANGE ON INPUT.
C RETURN CODE I HAS PRECEDENCE OVER I+1 FOR I = 9, 10, 11.
C IV(MLSTGD) (I/O) SAVED VALUE OF IV(MODEL).
C IV(MODEL) (I/O) ON INPUT, IV(MODEL) SHOULD BE AN INTEGER IDENTIFYING
C THE CURRENT QUADRATIC MODEL OF THE OBJECTIVE FUNCTION.
C IF A PREVIOUS STEP YIELDED A BETTER FUNCTION REDUCTION,
C THEN IV(MODEL) WILL BE SET TO IV(MLSTGD) ON OUTPUT.
C IV(NFCALL) (IN) INVOCATION COUNT FOR THE OBJECTIVE FUNCTION.
C IV(NFGCAL) (I/O) VALUE OF IV(NFCALL) AT STEP THAT GAVE THE BIGGEST
C FUNCTION REDUCTION THIS ITERATION. IV(NFGCAL) REMAINS
C UNCHANGED UNTIL A FUNCTION REDUCTION IS OBTAINED.
C IV(RADINC) (I/O) THE NUMBER OF RADIUS INCREASES (OR MINUS THE NUMBER
C OF DECREASES) SO FAR THIS ITERATION.
C IV(RESTOR) (OUT) SET TO 1 IF V(F) HAS BEEN RESTORED AND X SHOULD BE
C RESTORED TO ITS INITIAL VALUE, TO 2 IF X SHOULD BE SAVED,
C TO 3 IF X SHOULD BE RESTORED FROM THE SAVED VALUE, AND TO
C 0 OTHERWISE.
C IV(STAGE) (I/O) COUNT OF THE NUMBER OF MODELS TRIED SO FAR IN THE
C CURRENT ITERATION.
C IV(STGLIM) (IN) MAXIMUM NUMBER OF MODELS TO CONSIDER.
C IV(SWITCH) (OUT) SET TO 0 UNLESS A NEW MODEL IS BEING TRIED AND IT
C GIVES A SMALLER FUNCTION VALUE THAN THE PREVIOUS MODEL,
C IN WHICH CASE DA7SST SETS IV(SWITCH) = 1.
C IV(TOOBIG) (I/O) IS NONZERO ON INPUT IF STEP WAS TOO BIG (E.G., IF
C IT WOULD CAUSE OVERFLOW). IT IS SET TO 0 ON RETURN.
C IV(XIRC) (I/O) VALUE THAT IV(IRC) WOULD HAVE IN THE ABSENCE OF
C CONVERGENCE, FALSE CONVERGENCE, AND OVERSIZED STEPS.
C
C *** V VALUES REFERENCED ***
C
C V(AFCTOL) (IN) ABSOLUTE FUNCTION CONVERGENCE TOLERANCE. IF THE
C ABSOLUTE VALUE OF THE CURRENT FUNCTION VALUE V(F) IS LESS
C THAN V(AFCTOL) AND DA7SST DOES NOT RETURN WITH
C IV(IRC) = 11, THEN DA7SST RETURNS WITH IV(IRC) = 10.
C V(DECFAC) (IN) FACTOR BY WHICH TO DECREASE RADIUS WHEN IV(TOOBIG) IS
C NONZERO.
C V(DSTNRM) (IN) THE 2-NORM OF D*STEP.
C V(DSTSAV) (I/O) VALUE OF V(DSTNRM) ON SAVED STEP.
C V(DST0) (IN) THE 2-NORM OF D TIMES THE NEWTON STEP (WHEN DEFINED,
C I.E., FOR V(NREDUC) .GE. 0).
C V(F) (I/O) ON BOTH INPUT AND OUTPUT, V(F) IS THE OBJECTIVE FUNC-
C TION VALUE AT X. IF X IS RESTORED TO A PREVIOUS VALUE,
C THEN V(F) IS RESTORED TO THE CORRESPONDING VALUE.
C V(FDIF) (OUT) THE FUNCTION REDUCTION V(F0) - V(F) (FOR THE OUTPUT
C VALUE OF V(F) IF AN EARLIER STEP GAVE A BIGGER FUNCTION
C DECREASE, AND FOR THE INPUT VALUE OF V(F) OTHERWISE).
C V(FLSTGD) (I/O) SAVED VALUE OF V(F).
C V(F0) (IN) OBJECTIVE FUNCTION VALUE AT START OF ITERATION.
C V(GTSLST) (I/O) VALUE OF V(GTSTEP) ON SAVED STEP.
C V(GTSTEP) (IN) INNER PRODUCT BETWEEN STEP AND GRADIENT.
C V(INCFAC) (IN) MINIMUM FACTOR BY WHICH TO INCREASE RADIUS.
C V(LMAXS) (IN) MAXIMUM REASONABLE STEP SIZE (AND INITIAL STEP BOUND).
C IF THE ACTUAL FUNCTION DECREASE IS NO MORE THAN TWICE
C WHAT WAS PREDICTED, IF A RETURN WITH IV(IRC) = 7, 8, OR 9
C DOES NOT OCCUR, IF V(DSTNRM) .GT. V(LMAXS) OR THE CURRENT
C STEP IS A NEWTON STEP, AND IF
C V(PREDUC) .LE. V(SCTOL) * ABS(V(F0)), THEN DA7SST RETURNS
C WITH IV(IRC) = 11. IF SO DOING APPEARS WORTHWHILE, THEN
C DA7SST REPEATS THIS TEST (DISALLOWING A FULL NEWTON STEP)
C WITH V(PREDUC) COMPUTED FOR A STEP OF LENGTH V(LMAXS)
C (BY A RETURN WITH IV(IRC) = 6).
C V(NREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR
C NEWTON STEP. IF DA7SST IS CALLED WITH IV(IRC) = 6, I.E.,
C IF V(PREDUC) HAS BEEN COMPUTED WITH RADIUS = V(LMAXS) FOR
C USE IN THE SINGULAR CONVERGENCE TEST, THEN V(NREDUC) IS
C SET TO -V(PREDUC) BEFORE THE LATTER IS RESTORED.
C V(PLSTGD) (I/O) VALUE OF V(PREDUC) ON SAVED STEP.
C V(PREDUC) (I/O) FUNCTION REDUCTION PREDICTED BY QUADRATIC MODEL FOR
C CURRENT STEP.
C V(RADFAC) (OUT) FACTOR TO BE USED IN DETERMINING THE NEW RADIUS,
C WHICH SHOULD BE V(RADFAC)*DST, WHERE DST IS EITHER THE
C OUTPUT VALUE OF V(DSTNRM) OR THE 2-NORM OF
C DIAG(NEWD)*STEP FOR THE OUTPUT VALUE OF STEP AND THE
C UPDATED VERSION, NEWD, OF THE SCALE VECTOR D. FOR
C IV(IRC) = 3, V(RADFAC) = 1.0 IS RETURNED.
C V(RDFCMN) (IN) MINIMUM VALUE FOR V(RADFAC) IN TERMS OF THE INPUT
C VALUE OF V(DSTNRM) -- SUGGESTED VALUE = 0.1.
C V(RDFCMX) (IN) MAXIMUM VALUE FOR V(RADFAC) -- SUGGESTED VALUE = 4.0.
C V(RELDX) (IN) SCALED RELATIVE CHANGE IN X CAUSED BY STEP, COMPUTED
C (E.G.) BY FUNCTION DRLDST AS
C MAX (D(I)*ABS(X(I)-X0(I)), 1 .LE. I .LE. P) /
C MAX (D(I)*(ABS(X(I))+ABS(X0(I))), 1 .LE. I .LE. P).
C V(RFCTOL) (IN) RELATIVE FUNCTION CONVERGENCE TOLERANCE. IF THE
C ACTUAL FUNCTION REDUCTION IS AT MOST TWICE WHAT WAS PRE-
C DICTED AND V(NREDUC) .LE. V(RFCTOL)*ABS(V(F0)), THEN
C DA7SST RETURNS WITH IV(IRC) = 8 OR 9.
C V(SCTOL) (IN) SINGULAR CONVERGENCE TOLERANCE -- SEE V(LMAXS).
C V(STPPAR) (IN) MARQUARDT PARAMETER -- 0 MEANS FULL NEWTON STEP.
C V(TUNER1) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION
C REDUCTION WAS MUCH LESS THAN EXPECTED. SUGGESTED
C VALUE = 0.1.
C V(TUNER2) (IN) TUNING CONSTANT USED TO DECIDE IF THE FUNCTION
C REDUCTION WAS LARGE ENOUGH TO ACCEPT STEP. SUGGESTED
C VALUE = 10**-4.
C V(TUNER3) (IN) TUNING CONSTANT USED TO DECIDE IF THE RADIUS
C SHOULD BE INCREASED. SUGGESTED VALUE = 0.75.
C V(XCTOL) (IN) X-CONVERGENCE CRITERION. IF STEP IS A NEWTON STEP
C (V(STPPAR) = 0) HAVING V(RELDX) .LE. V(XCTOL) AND GIVING
C AT MOST TWICE THE PREDICTED FUNCTION DECREASE, THEN
C DA7SST RETURNS IV(IRC) = 7 OR 9.
C V(XFTOL) (IN) FALSE CONVERGENCE TOLERANCE. IF STEP GAVE NO OR ONLY
C A SMALL FUNCTION DECREASE AND V(RELDX) .LE. V(XFTOL),
C THEN DA7SST RETURNS WITH IV(IRC) = 12.
C
C------------------------------- NOTES -------------------------------
C
C *** APPLICATION AND USAGE RESTRICTIONS ***
C
C THIS ROUTINE IS CALLED AS PART OF THE NL2SOL (NONLINEAR
C LEAST-SQUARES) PACKAGE. IT MAY BE USED IN ANY UNCONSTRAINED
C MINIMIZATION SOLVER THAT USES DOGLEG, GOLDFELD-QUANDT-TROTTER,
C OR LEVENBERG-MARQUARDT STEPS.
C
C *** ALGORITHM NOTES ***
C
C SEE (1) FOR FURTHER DISCUSSION OF THE ASSESSING AND MODEL
C SWITCHING STRATEGIES. WHILE NL2SOL CONSIDERS ONLY TWO MODELS,
C DA7SST IS DESIGNED TO HANDLE ANY NUMBER OF MODELS.
C
C *** USAGE NOTES ***
C
C ON THE FIRST CALL OF AN ITERATION, ONLY THE I/O VARIABLES
C STEP, X, IV(IRC), IV(MODEL), V(F), V(DSTNRM), V(GTSTEP), AND
C V(PREDUC) NEED HAVE BEEN INITIALIZED. BETWEEN CALLS, NO I/O
C VALUES EXCEPT STEP, X, IV(MODEL), V(F) AND THE STOPPING TOLER-
C ANCES SHOULD BE CHANGED.
C AFTER A RETURN FOR CONVERGENCE OR FALSE CONVERGENCE, ONE CAN
C CHANGE THE STOPPING TOLERANCES AND CALL DA7SST AGAIN, IN WHICH
C CASE THE STOPPING TESTS WILL BE REPEATED.
C
C *** REFERENCES ***
C
C (1) DENNIS, J.E., JR., GAY, D.M., AND WELSCH, R.E. (1981),
C AN ADAPTIVE NONLINEAR LEAST-SQUARES ALGORITHM,
C ACM TRANS. MATH. SOFTWARE, VOL. 7, NO. 3.
C
C (2) POWELL, M.J.D. (1970) A FORTRAN SUBROUTINE FOR SOLVING
C SYSTEMS OF NONLINEAR ALGEBRAIC EQUATIONS, IN NUMERICAL
C METHODS FOR NONLINEAR ALGEBRAIC EQUATIONS, EDITED BY
C P. RABINOWITZ, GORDON AND BREACH, LONDON.
C
C *** HISTORY ***
C
C JOHN DENNIS DESIGNED MUCH OF THIS ROUTINE, STARTING WITH
C IDEAS IN (2). ROY WELSCH SUGGESTED THE MODEL SWITCHING STRATEGY.
C DAVID GAY AND STEPHEN PETERS CAST THIS SUBROUTINE INTO A MORE
C PORTABLE FORM (WINTER 1977), AND DAVID GAY CAST IT INTO ITS
C PRESENT FORM (FALL 1978), WITH MINOR CHANGES TO THE SINGULAR
C CONVERGENCE TEST IN MAY, 1984 (TO DEAL WITH FULL NEWTON STEPS).
C
C *** GENERAL ***
C
C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH
C SUPPORTED BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS
C MCS-7600324, DCR75-10143, 76-14311DSS, MCS76-11989, AND
C MCS-7906671.
C
C------------------------ EXTERNAL QUANTITIES ------------------------
C
C *** NO EXTERNAL FUNCTIONS AND SUBROUTINES ***
C
C-------------------------- LOCAL VARIABLES --------------------------
C
LOGICAL GOODX
INTEGER I, NFC
DOUBLE PRECISION EMAX, EMAXS, GTS, RFAC1, XMAX
DOUBLE PRECISION HALF, ONE, ONEP2, TWO, ZERO
C
C *** SUBSCRIPTS FOR IV AND V ***
C
INTEGER AFCTOL, DECFAC, DSTNRM, DSTSAV, DST0, F, FDIF, FLSTGD, F0,
1 GTSLST, GTSTEP, INCFAC, IRC, LMAXS, MLSTGD, MODEL, NFCALL,
2 NFGCAL, NREDUC, PLSTGD, PREDUC, RADFAC, RADINC, RDFCMN,
3 RDFCMX, RELDX, RESTOR, RFCTOL, SCTOL, STAGE, STGLIM,
4 STPPAR, SWITCH, TOOBIG, TUNER1, TUNER2, TUNER3, XCTOL,
5 XFTOL, XIRC
C
C *** DATA INITIALIZATIONS ***
C
C/6
C DATA HALF/0.5D+0/, ONE/1.D+0/, ONEP2/1.2D+0/, TWO/2.D+0/,
C 1 ZERO/0.D+0/
C/7
PARAMETER (HALF=0.5D+0, ONE=1.D+0, ONEP2=1.2D+0, TWO=2.D+0,
1 ZERO=0.D+0)
C/
C
C/6
C DATA IRC/29/, MLSTGD/32/, MODEL/5/, NFCALL/6/, NFGCAL/7/,
C 1 RADINC/8/, RESTOR/9/, STAGE/10/, STGLIM/11/, SWITCH/12/,
C 2 TOOBIG/2/, XIRC/13/
C/7
PARAMETER (IRC=29, MLSTGD=32, MODEL=5, NFCALL=6, NFGCAL=7,
1 RADINC=8, RESTOR=9, STAGE=10, STGLIM=11, SWITCH=12,
2 TOOBIG=2, XIRC=13)
C/
C/6
C DATA AFCTOL/31/, DECFAC/22/, DSTNRM/2/, DST0/3/, DSTSAV/18/,
C 1 F/10/, FDIF/11/, FLSTGD/12/, F0/13/, GTSLST/14/, GTSTEP/4/,
C 2 INCFAC/23/, LMAXS/36/, NREDUC/6/, PLSTGD/15/, PREDUC/7/,
C 3 RADFAC/16/, RDFCMN/24/, RDFCMX/25/, RELDX/17/, RFCTOL/32/,
C 4 SCTOL/37/, STPPAR/5/, TUNER1/26/, TUNER2/27/, TUNER3/28/,
C 5 XCTOL/33/, XFTOL/34/
C/7
PARAMETER (AFCTOL=31, DECFAC=22, DSTNRM=2, DST0=3, DSTSAV=18,
1 F=10, FDIF=11, FLSTGD=12, F0=13, GTSLST=14, GTSTEP=4,
2 INCFAC=23, LMAXS=36, NREDUC=6, PLSTGD=15, PREDUC=7,
3 RADFAC=16, RDFCMN=24, RDFCMX=25, RELDX=17, RFCTOL=32,
4 SCTOL=37, STPPAR=5, TUNER1=26, TUNER2=27, TUNER3=28,
5 XCTOL=33, XFTOL=34)
C/
C
C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
C
NFC = IV(NFCALL)
IV(SWITCH) = 0
IV(RESTOR) = 0
RFAC1 = ONE
GOODX = .TRUE.
I = IV(IRC)
IF (I .GE. 1 .AND. I .LE. 12)
1 GO TO (20,30,10,10,40,280,220,220,220,220,220,170), I
IV(IRC) = 13
GO TO 999
C
C *** INITIALIZE FOR NEW ITERATION ***
C
10 IV(STAGE) = 1
IV(RADINC) = 0
V(FLSTGD) = V(F0)
IF (IV(TOOBIG) .EQ. 0) GO TO 110
IV(STAGE) = -1
IV(XIRC) = I
GO TO 60
C
C *** STEP WAS RECOMPUTED WITH NEW MODEL OR SMALLER RADIUS ***
C *** FIRST DECIDE WHICH ***
C
20 IF (IV(MODEL) .NE. IV(MLSTGD)) GO TO 30
C *** OLD MODEL RETAINED, SMALLER RADIUS TRIED ***
C *** DO NOT CONSIDER ANY MORE NEW MODELS THIS ITERATION ***
IV(STAGE) = IV(STGLIM)
IV(RADINC) = -1
GO TO 110
C
C *** A NEW MODEL IS BEING TRIED. DECIDE WHETHER TO KEEP IT. ***
C
30 IV(STAGE) = IV(STAGE) + 1
C
C *** NOW WE ADD THE POSSIBILITY THAT STEP WAS RECOMPUTED WITH ***
C *** THE SAME MODEL, PERHAPS BECAUSE OF AN OVERSIZED STEP. ***
C
40 IF (IV(STAGE) .GT. 0) GO TO 50
C
C *** STEP WAS RECOMPUTED BECAUSE IT WAS TOO BIG. ***
C
IF (IV(TOOBIG) .NE. 0) GO TO 60
C
C *** RESTORE IV(STAGE) AND PICK UP WHERE WE LEFT OFF. ***
C
IV(STAGE) = -IV(STAGE)
I = IV(XIRC)
GO TO (20, 30, 110, 110, 70), I
C
50 IF (IV(TOOBIG) .EQ. 0) GO TO 70
C
C *** HANDLE OVERSIZE STEP ***
C
IV(TOOBIG) = 0
IF (IV(RADINC) .GT. 0) GO TO 80
IV(STAGE) = -IV(STAGE)
IV(XIRC) = IV(IRC)
C
60 IV(TOOBIG) = 0
V(RADFAC) = V(DECFAC)
IV(RADINC) = IV(RADINC) - 1
IV(IRC) = 5
IV(RESTOR) = 1
V(F) = V(FLSTGD)
GO TO 999
C
70 IF (V(F) .LT. V(FLSTGD)) GO TO 110
C
C *** THE NEW STEP IS A LOSER. RESTORE OLD MODEL. ***
C
IF (IV(MODEL) .EQ. IV(MLSTGD)) GO TO 80
IV(MODEL) = IV(MLSTGD)
IV(SWITCH) = 1
C
C *** RESTORE STEP, ETC. ONLY IF A PREVIOUS STEP DECREASED V(F).
C
80 IF (V(FLSTGD) .GE. V(F0)) GO TO 110
IV(RESTOR) = 1
V(F) = V(FLSTGD)
V(PREDUC) = V(PLSTGD)
V(GTSTEP) = V(GTSLST)
IF (IV(SWITCH) .EQ. 0) RFAC1 = V(DSTNRM) / V(DSTSAV)
V(DSTNRM) = V(DSTSAV)
NFC = IV(NFGCAL)
GOODX = .FALSE.
C
110 V(FDIF) = V(F0) - V(F)
IF (V(FDIF) .GT. V(TUNER2) * V(PREDUC)) GO TO 140
IF (IV(RADINC) .GT. 0) GO TO 140
C
C *** NO (OR ONLY A TRIVIAL) FUNCTION DECREASE
C *** -- SO TRY NEW MODEL OR SMALLER RADIUS
C
IF (V(F) .LT. V(F0)) GO TO 120
IV(MLSTGD) = IV(MODEL)
V(FLSTGD) = V(F)
V(F) = V(F0)
IV(RESTOR) = 1
GO TO 130
120 IV(NFGCAL) = NFC
130 IV(IRC) = 1
IF (IV(STAGE) .LT. IV(STGLIM)) GO TO 160
IV(IRC) = 5
IV(RADINC) = IV(RADINC) - 1
GO TO 160
C
C *** NONTRIVIAL FUNCTION DECREASE ACHIEVED ***
C
140 IV(NFGCAL) = NFC
RFAC1 = ONE
V(DSTSAV) = V(DSTNRM)
IF (V(FDIF) .GT. V(PREDUC)*V(TUNER1)) GO TO 190
C
C *** DECREASE WAS MUCH LESS THAN PREDICTED -- EITHER CHANGE MODELS
C *** OR ACCEPT STEP WITH DECREASED RADIUS.
C
IF (IV(STAGE) .GE. IV(STGLIM)) GO TO 150
C *** CONSIDER SWITCHING MODELS ***
IV(IRC) = 2
GO TO 160
C
C *** ACCEPT STEP WITH DECREASED RADIUS ***
C
150 IV(IRC) = 4
C
C *** SET V(RADFAC) TO FLETCHER*S DECREASE FACTOR ***
C
160 IV(XIRC) = IV(IRC)
EMAX = V(GTSTEP) + V(FDIF)
V(RADFAC) = HALF * RFAC1
IF (EMAX .LT. V(GTSTEP)) V(RADFAC) = RFAC1 * DMAX1(V(RDFCMN),
1 HALF * V(GTSTEP)/EMAX)
C
C *** DO FALSE CONVERGENCE TEST ***
C
170 IF (V(RELDX) .LE. V(XFTOL)) GO TO 180
IV(IRC) = IV(XIRC)
IF (V(F) .LT. V(F0)) GO TO 200
GO TO 230
C
180 IV(IRC) = 12
GO TO 240
C
C *** HANDLE GOOD FUNCTION DECREASE ***
C
190 IF (V(FDIF) .LT. (-V(TUNER3) * V(GTSTEP))) GO TO 210
C
C *** INCREASING RADIUS LOOKS WORTHWHILE. SEE IF WE JUST
C *** RECOMPUTED STEP WITH A DECREASED RADIUS OR RESTORED STEP
C *** AFTER RECOMPUTING IT WITH A LARGER RADIUS.
C
IF (IV(RADINC) .LT. 0) GO TO 210
IF (IV(RESTOR) .EQ. 1) GO TO 210
C
C *** WE DID NOT. TRY A LONGER STEP UNLESS THIS WAS A NEWTON
C *** STEP.
C
V(RADFAC) = V(RDFCMX)
GTS = V(GTSTEP)
IF (V(FDIF) .LT. (HALF/V(RADFAC) - ONE) * GTS)
1 V(RADFAC) = DMAX1(V(INCFAC), HALF*GTS/(GTS + V(FDIF)))
IV(IRC) = 4
IF (V(STPPAR) .EQ. ZERO) GO TO 230
IF (V(DST0) .GE. ZERO .AND. (V(DST0) .LT. TWO*V(DSTNRM)
1 .OR. V(NREDUC) .LT. ONEP2*V(FDIF))) GO TO 230
C *** STEP WAS NOT A NEWTON STEP. RECOMPUTE IT WITH
C *** A LARGER RADIUS.
IV(IRC) = 5
IV(RADINC) = IV(RADINC) + 1
C
C *** SAVE VALUES CORRESPONDING TO GOOD STEP ***
C
200 V(FLSTGD) = V(F)
IV(MLSTGD) = IV(MODEL)
IF (IV(RESTOR) .NE. 1) IV(RESTOR) = 2
V(DSTSAV) = V(DSTNRM)
IV(NFGCAL) = NFC
V(PLSTGD) = V(PREDUC)
V(GTSLST) = V(GTSTEP)
GO TO 230
C
C *** ACCEPT STEP WITH RADIUS UNCHANGED ***
C
210 V(RADFAC) = ONE
IV(IRC) = 3
GO TO 230
C
C *** COME HERE FOR A RESTART AFTER CONVERGENCE ***
C
220 IV(IRC) = IV(XIRC)
IF (V(DSTSAV) .GE. ZERO) GO TO 240
IV(IRC) = 12
GO TO 240
C
C *** PERFORM CONVERGENCE TESTS ***
C
230 IV(XIRC) = IV(IRC)
240 IF (IV(RESTOR) .EQ. 1 .AND. V(FLSTGD) .LT. V(F0)) IV(RESTOR) = 3
IF (DABS(V(F)) .LT. V(AFCTOL)) IV(IRC) = 10
IF (HALF * V(FDIF) .GT. V(PREDUC)) GO TO 999
EMAX = V(RFCTOL) * DABS(V(F0))
EMAXS = V(SCTOL) * DABS(V(F0))
IF (V(PREDUC) .LE. EMAXS .AND. (V(DSTNRM) .GT. V(LMAXS) .OR.
1 V(STPPAR) .EQ. ZERO)) IV(IRC) = 11
IF (V(DST0) .LT. ZERO) GO TO 250
I = 0
IF ((V(NREDUC) .GT. ZERO .AND. V(NREDUC) .LE. EMAX) .OR.
1 (V(NREDUC) .EQ. ZERO. AND. V(PREDUC) .EQ. ZERO)) I = 2
IF (V(STPPAR) .EQ. ZERO .AND. V(RELDX) .LE. V(XCTOL)
1 .AND. GOODX) I = I + 1
IF (I .GT. 0) IV(IRC) = I + 6
C
C *** CONSIDER RECOMPUTING STEP OF LENGTH V(LMAXS) FOR SINGULAR
C *** CONVERGENCE TEST.
C
250 IF (IV(IRC) .GT. 5 .AND. IV(IRC) .NE. 12) GO TO 999
IF (V(STPPAR) .EQ. ZERO) GO TO 999
IF (V(DSTNRM) .GT. V(LMAXS)) GO TO 260
IF (V(PREDUC) .GE. EMAXS) GO TO 999
IF (V(DST0) .LE. ZERO) GO TO 270
IF (HALF * V(DST0) .LE. V(LMAXS)) GO TO 999
GO TO 270
260 IF (HALF * V(DSTNRM) .LE. V(LMAXS)) GO TO 999
XMAX = V(LMAXS) / V(DSTNRM)
IF (XMAX * (TWO - XMAX) * V(PREDUC) .GE. EMAXS) GO TO 999
270 IF (V(NREDUC) .LT. ZERO) GO TO 290
C
C *** RECOMPUTE V(PREDUC) FOR USE IN SINGULAR CONVERGENCE TEST ***
C
V(GTSLST) = V(GTSTEP)
V(DSTSAV) = V(DSTNRM)
IF (IV(IRC) .EQ. 12) V(DSTSAV) = -V(DSTSAV)
V(PLSTGD) = V(PREDUC)
I = IV(RESTOR)
IV(RESTOR) = 2
IF (I .EQ. 3) IV(RESTOR) = 0
IV(IRC) = 6
GO TO 999
C
C *** PERFORM SINGULAR CONVERGENCE TEST WITH RECOMPUTED V(PREDUC) ***
C
280 V(GTSTEP) = V(GTSLST)
V(DSTNRM) = DABS(V(DSTSAV))
IV(IRC) = IV(XIRC)
IF (V(DSTSAV) .LE. ZERO) IV(IRC) = 12
V(NREDUC) = -V(PREDUC)
V(PREDUC) = V(PLSTGD)
IV(RESTOR) = 3
290 IF (-V(NREDUC) .LE. V(SCTOL) * DABS(V(F0))) IV(IRC) = 11
C
999 RETURN
C
C *** LAST LINE OF DA7SST FOLLOWS ***
END
SUBROUTINE DD7DOG(DIG, LV, N, NWTSTP, STEP, V)
C
C *** COMPUTE DOUBLE DOGLEG STEP ***
C
C *** PARAMETER DECLARATIONS ***
C
INTEGER LV, N
DOUBLE PRECISION DIG(N), NWTSTP(N), STEP(N), V(LV)
C
C *** PURPOSE ***
C
C THIS SUBROUTINE COMPUTES A CANDIDATE STEP (FOR USE IN AN UNCON-
C STRAINED MINIMIZATION CODE) BY THE DOUBLE DOGLEG ALGORITHM OF
C DENNIS AND MEI (REF. 1), WHICH IS A VARIATION ON POWELL*S DOGLEG
C SCHEME (REF. 2, P. 95).
C
C-------------------------- PARAMETER USAGE --------------------------
C
C DIG (INPUT) DIAG(D)**-2 * G -- SEE ALGORITHM NOTES.
C G (INPUT) THE CURRENT GRADIENT VECTOR.
C LV (INPUT) LENGTH OF V.
C N (INPUT) NUMBER OF COMPONENTS IN DIG, G, NWTSTP, AND STEP.
C NWTSTP (INPUT) NEGATIVE NEWTON STEP -- SEE ALGORITHM NOTES.
C STEP (OUTPUT) THE COMPUTED STEP.
C V (I/O) VALUES ARRAY, THE FOLLOWING COMPONENTS OF WHICH ARE
C USED HERE...
C V(BIAS) (INPUT) BIAS FOR RELAXED NEWTON STEP, WHICH IS V(BIAS) OF
C THE WAY FROM THE FULL NEWTON TO THE FULLY RELAXED NEWTON
C STEP. RECOMMENDED VALUE = 0.8 .
C V(DGNORM) (INPUT) 2-NORM OF DIAG(D)**-1 * G -- SEE ALGORITHM NOTES.
C V(DSTNRM) (OUTPUT) 2-NORM OF DIAG(D) * STEP, WHICH IS V(RADIUS)
C UNLESS V(STPPAR) = 0 -- SEE ALGORITHM NOTES.
C V(DST0) (INPUT) 2-NORM OF DIAG(D) * NWTSTP -- SEE ALGORITHM NOTES.
C V(GRDFAC) (OUTPUT) THE COEFFICIENT OF DIG IN THE STEP RETURNED --
C STEP(I) = V(GRDFAC)*DIG(I) + V(NWTFAC)*NWTSTP(I).
C V(GTHG) (INPUT) SQUARE-ROOT OF (DIG**T) * (HESSIAN) * DIG -- SEE
C ALGORITHM NOTES.
C V(GTSTEP) (OUTPUT) INNER PRODUCT BETWEEN G AND STEP.
C V(NREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE FULL NEWTON
C STEP.
C V(NWTFAC) (OUTPUT) THE COEFFICIENT OF NWTSTP IN THE STEP RETURNED --
C SEE V(GRDFAC) ABOVE.
C V(PREDUC) (OUTPUT) FUNCTION REDUCTION PREDICTED FOR THE STEP RETURNED.
C V(RADIUS) (INPUT) THE TRUST REGION RADIUS. D TIMES THE STEP RETURNED
C HAS 2-NORM V(RADIUS) UNLESS V(STPPAR) = 0.
C V(STPPAR) (OUTPUT) CODE TELLING HOW STEP WAS COMPUTED... 0 MEANS A
C FULL NEWTON STEP. BETWEEN 0 AND 1 MEANS V(STPPAR) OF THE
C WAY FROM THE NEWTON TO THE RELAXED NEWTON STEP. BETWEEN
C 1 AND 2 MEANS A TRUE DOUBLE DOGLEG STEP, V(STPPAR) - 1 OF
C THE WAY FROM THE RELAXED NEWTON TO THE CAUCHY STEP.
C GREATER THAN 2 MEANS 1 / (V(STPPAR) - 1) TIMES THE CAUCHY
C STEP.
C
C------------------------------- NOTES -------------------------------
C
C *** ALGORITHM NOTES ***
C
C LET G AND H BE THE CURRENT GRADIENT AND HESSIAN APPROXIMA-
C TION RESPECTIVELY AND LET D BE THE CURRENT SCALE VECTOR. THIS
C ROUTINE ASSUMES DIG = DIAG(D)**-2 * G AND NWTSTP = H**-1 * G.
C THE STEP COMPUTED IS THE SAME ONE WOULD GET BY REPLACING G AND H
C BY DIAG(D)**-1 * G AND DIAG(D)**-1 * H * DIAG(D)**-1,
C COMPUTING STEP, AND TRANSLATING STEP BACK TO THE ORIGINAL
C VARIABLES, I.E., PREMULTIPLYING IT BY DIAG(D)**-1.
C
C *** REFERENCES ***
C
C 1. DENNIS, J.E., AND MEI, H.H.W. (1979), TWO NEW UNCONSTRAINED OPTI-
C MIZATION ALGORITHMS WHICH USE FUNCTION AND GRADIENT
C VALUES, J. OPTIM. THEORY APPLIC. 28, PP. 453-482.
C 2. POWELL, M.J.D. (1970), A HYBRID METHOD FOR NON-LINEAR EQUATIONS,
C IN NUMERICAL METHODS FOR NON-LINEAR EQUATIONS, EDITED BY
C P. RABINOWITZ, GORDON AND BREACH, LONDON.
C
C *** GENERAL ***
C
C CODED BY DAVID M. GAY.
C THIS SUBROUTINE WAS WRITTEN IN CONNECTION WITH RESEARCH SUPPORTED
C BY THE NATIONAL SCIENCE FOUNDATION UNDER GRANTS MCS-7600324 AND
C MCS-7906671.
C
C------------------------ EXTERNAL QUANTITIES ------------------------
C
C *** INTRINSIC FUNCTIONS ***
C/+
DOUBLE PRECISION DSQRT
C/
C-------------------------- LOCAL VARIABLES --------------------------
C
INTEGER I
DOUBLE PRECISION CFACT, CNORM, CTRNWT, GHINVG, FEMNSQ, GNORM,
1 NWTNRM, RELAX, RLAMBD, T, T1, T2
DOUBLE PRECISION HALF, ONE, TWO, ZERO
C
C *** V SUBSCRIPTS ***
C
INTEGER BIAS, DGNORM, DSTNRM, DST0, GRDFAC, GTHG, GTSTEP,
1 NREDUC, NWTFAC, PREDUC, RADIUS, STPPAR
C
C *** DATA INITIALIZATIONS ***
C
C/6
C DATA HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/, ZERO/0.D+0/
C/7
PARAMETER (HALF=0.5D+0, ONE=1.D+0, TWO=2.D+0, ZERO=0.D+0)
C/
C
C/6
C DATA BIAS/43/, DGNORM/1/, DSTNRM/2/, DST0/3/, GRDFAC/45/,
C 1 GTHG/44/, GTSTEP/4/, NREDUC/6/, NWTFAC/46/, PREDUC/7/,
C 2 RADIUS/8/, STPPAR/5/
C/7
PARAMETER (BIAS=43, DGNORM=1, DSTNRM=2, DST0=3, GRDFAC=45,
1 GTHG=44, GTSTEP=4, NREDUC=6, NWTFAC=46, PREDUC=7,
2 RADIUS=8, STPPAR=5)
C/
C
C+++++++++++++++++++++++++++++++ BODY ++++++++++++++++++++++++++++++++
C
NWTNRM = V(DST0)
RLAMBD = ONE
IF (NWTNRM .GT. ZERO) RLAMBD = V(RADIUS) / NWTNRM
GNORM = V(DGNORM)
GHINVG = TWO * V(NREDUC)
V(GRDFAC) = ZERO
V(NWTFAC) = ZERO
IF (RLAMBD .LT. ONE) GO TO 30
C
C *** THE NEWTON STEP IS INSIDE THE TRUST REGION ***
C
V(STPPAR) = ZERO
V(DSTNRM) = NWTNRM
V(GTSTEP) = -GHINVG
V(PREDUC) = V(NREDUC)
V(NWTFAC) = -ONE
DO 20 I = 1, N
20 STEP(I) = -NWTSTP(I)
GO TO 999
C
30 V(DSTNRM) = V(RADIUS)
CFACT = (GNORM / V(GTHG))**2
C *** CAUCHY STEP = -CFACT * G.
CNORM = GNORM * CFACT
RELAX = ONE - V(BIAS) * (ONE - GNORM*CNORM/GHINVG)
IF (RLAMBD .LT. RELAX) GO TO 50
C
C *** STEP IS BETWEEN RELAXED NEWTON AND FULL NEWTON STEPS ***
C
V(STPPAR) = ONE - (RLAMBD - RELAX) / (ONE - RELAX)
T = -RLAMBD
V(GTSTEP) = T * GHINVG
V(PREDUC) = RLAMBD * (ONE - HALF*RLAMBD) * GHINVG
V(NWTFAC) = T
DO 40 I = 1, N
40 STEP(I) = T * NWTSTP(I)
GO TO 999
C
50 IF (CNORM .LT. V(RADIUS)) GO TO 70
C
C *** THE CAUCHY STEP LIES OUTSIDE THE TRUST REGION --
C *** STEP = SCALED CAUCHY STEP ***
C
T = -V(RADIUS) / GNORM
V(GRDFAC) = T
V(STPPAR) = ONE + CNORM / V(RADIUS)
V(GTSTEP) = -V(RADIUS) * GNORM
V(PREDUC) = V(RADIUS)*(GNORM - HALF*V(RADIUS)*(V(GTHG)/GNORM)**2)
DO 60 I = 1, N
60 STEP(I) = T * DIG(I)
GO TO 999
C
C *** COMPUTE DOGLEG STEP BETWEEN CAUCHY AND RELAXED NEWTON ***
C *** FEMUR = RELAXED NEWTON STEP MINUS CAUCHY STEP ***
C
70 CTRNWT = CFACT * RELAX * GHINVG / GNORM
C *** CTRNWT = INNER PROD. OF CAUCHY AND RELAXED NEWTON STEPS,
C *** SCALED BY GNORM**-1.
T1 = CTRNWT - GNORM*CFACT**2
C *** T1 = INNER PROD. OF FEMUR AND CAUCHY STEP, SCALED BY
C *** GNORM**-1.
T2 = V(RADIUS)*(V(RADIUS)/GNORM) - GNORM*CFACT**2
T = RELAX * NWTNRM
FEMNSQ = (T/GNORM)*T - CTRNWT - T1
C *** FEMNSQ = SQUARE OF 2-NORM OF FEMUR, SCALED BY GNORM**-1.
T = T2 / (T1 + DSQRT(T1**2 + FEMNSQ*T2))
C *** DOGLEG STEP = CAUCHY STEP + T * FEMUR.
T1 = (T - ONE) * CFACT
V(GRDFAC) = T1
T2 = -T * RELAX
V(NWTFAC) = T2
V(STPPAR) = TWO - T
V(GTSTEP) = T1*GNORM**2 + T2*GHINVG
V(PREDUC) = -T1*GNORM * ((T2 + ONE)*GNORM)
1 - T2 * (ONE + HALF*T2)*GHINVG
2 - HALF * (V(GTHG)*T1)**2
DO 80 I = 1, N
80 STEP(I) = T1*DIG(I) + T2*NWTSTP(I)
C
999 RETURN
C *** LAST LINE OF DD7DOG FOLLOWS ***
END
DOUBLE PRECISION FUNCTION DD7TPR(P, X, Y)
C
C *** RETURN THE INNER PRODUCT OF THE P-VECTORS X AND Y. ***
C
INTEGER P
DOUBLE PRECISION X(P), Y(P)
C
INTEGER I
DOUBLE PRECISION ONE, SQTETA, T, ZERO
DOUBLE PRECISION DR7MDC
EXTERNAL DR7MDC
C
C *** DR7MDC(2) RETURNS A MACHINE-DEPENDENT CONSTANT, SQTETA, WHICH
C *** IS SLIGHTLY LARGER THAN THE SMALLEST POSITIVE NUMBER THAT
C *** CAN BE SQUARED WITHOUT UNDERFLOWING.
C
C/6
C DATA ONE/1.D+0/, SQTETA/0.D+0/, ZERO/0.D+0/
C/7
PARAMETER (ONE=1.D+0, ZERO=0.D+0)
DATA SQTETA/0.D+0/
C/
C
DD7TPR = ZERO
IF (P .LE. 0) GO TO 999
IF (SQTETA .EQ. ZERO) SQTETA = DR7MDC(2)
DO 20 I = 1, P
T = DMAX1(DABS(X(I)), DABS(Y(I)))
IF (T .GT. ONE) GO TO 10
IF (T .LT. SQTETA) GO TO 20
T = (X(I)/SQTETA)*Y(I)
IF (DABS(T) .LT. SQTETA) GO TO 20
10 DD7TPR = DD7TPR + X(I)*Y(I)
20 CONTINUE
C
999 RETURN
C *** LAST LINE OF DD7TPR FOLLOWS ***
END
SUBROUTINE DITSUM(D, G, IV, LIV, LV, P, V, X)
C
C *** PRINT ITERATION SUMMARY FOR ***SOL (VERSION 2.3) ***
C
C *** PARAMETER DECLARATIONS ***
C
INTEGER LIV, LV, P
INTEGER IV(LIV)
DOUBLE PRECISION D(P), G(P), V(LV), X(P)
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C *** LOCAL VARIABLES ***
C
INTEGER ALG, I, IV1, M, NF, NG, OL, PU
C/6S
C REAL MODEL1(6), MODEL2(6)
C/7S
CHARACTER*4 MODEL1(6), MODEL2(6)
C/
DOUBLE PRECISION NRELDF, OLDF, PRELDF, RELDF, ZERO
C
C *** NO EXTERNAL FUNCTIONS OR SUBROUTINES ***
C
C *** SUBSCRIPTS FOR IV AND V ***
C
INTEGER ALGSAV, DSTNRM, F, FDIF, F0, NEEDHD, NFCALL, NFCOV, NGCOV,
1 NGCALL, NITER, NREDUC, OUTLEV, PREDUC, PRNTIT, PRUNIT,
2 RELDX, SOLPRT, STATPR, STPPAR, SUSED, X0PRT
C
C *** IV SUBSCRIPT VALUES ***
C
C/6
C DATA ALGSAV/51/, NEEDHD/36/, NFCALL/6/, NFCOV/52/, NGCALL/30/,
C 1 NGCOV/53/, NITER/31/, OUTLEV/19/, PRNTIT/39/, PRUNIT/21/,
C 2 SOLPRT/22/, STATPR/23/, SUSED/64/, X0PRT/24/
C/7
PARAMETER (ALGSAV=51, NEEDHD=36, NFCALL=6, NFCOV=52, NGCALL=30,
1 NGCOV=53, NITER=31, OUTLEV=19, PRNTIT=39, PRUNIT=21,
2 SOLPRT=22, STATPR=23, SUSED=64, X0PRT=24)
C/
C
C *** V SUBSCRIPT VALUES ***
C
C/6
C DATA DSTNRM/2/, F/10/, F0/13/, FDIF/11/, NREDUC/6/, PREDUC/7/,
C 1 RELDX/17/, STPPAR/5/
C/7
PARAMETER (DSTNRM=2, F=10, F0=13, FDIF=11, NREDUC=6, PREDUC=7,
1 RELDX=17, STPPAR=5)
C/
C
C/6
C DATA ZERO/0.D+0/
C/7
PARAMETER (ZERO=0.D+0)
C/
C/6S
C DATA MODEL1(1)/4H /, MODEL1(2)/4H /, MODEL1(3)/4H /,
C 1 MODEL1(4)/4H /, MODEL1(5)/4H G /, MODEL1(6)/4H S /,
C 2 MODEL2(1)/4H G /, MODEL2(2)/4H S /, MODEL2(3)/4HG-S /,
C 3 MODEL2(4)/4HS-G /, MODEL2(5)/4H-S-G/, MODEL2(6)/4H-G-S/
C/7S
DATA MODEL1/' ',' ',' ',' ',' G ',' S '/,
1 MODEL2/' G ',' S ','G-S ','S-G ','-S-G','-G-S'/
C/
C
C------------------------------- BODY --------------------------------
C
PU = IV(PRUNIT)
IF (PU .EQ. 0) GO TO 999
IV1 = IV(1)
IF (IV1 .GT. 62) IV1 = IV1 - 51
OL = IV(OUTLEV)
ALG = MOD(IV(ALGSAV)-1,2) + 1
IF (IV1 .LT. 2 .OR. IV1 .GT. 15) GO TO 370
IF (IV1 .GE. 12) GO TO 120
IF (IV1 .EQ. 2 .AND. IV(NITER) .EQ. 0) GO TO 390
IF (OL .EQ. 0) GO TO 120
IF (IV1 .GE. 10 .AND. IV(PRNTIT) .EQ. 0) GO TO 120
IF (IV1 .GT. 2) GO TO 10
IV(PRNTIT) = IV(PRNTIT) + 1
IF (IV(PRNTIT) .LT. IABS(OL)) GO TO 999
10 NF = IV(NFCALL) - IABS(IV(NFCOV))
IV(PRNTIT) = 0
RELDF = ZERO
PRELDF = ZERO
OLDF = DMAX1(DABS(V(F0)), DABS(V(F)))
IF (OLDF .LE. ZERO) GO TO 20
RELDF = V(FDIF) / OLDF
PRELDF = V(PREDUC) / OLDF
20 IF (OL .GT. 0) GO TO 60
C
C *** PRINT SHORT SUMMARY LINE ***
C
IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,30)
30 FORMAT(/10H IT NF,6X,1HF,7X,5HRELDF,3X,6HPRELDF,3X,5HRELDX,
1 2X,13HMODEL STPPAR)
IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,40)
40 FORMAT(/11H IT NF,7X,1HF,8X,5HRELDF,4X,6HPRELDF,4X,5HRELDX,
1 3X,6HSTPPAR)
IV(NEEDHD) = 0
IF (ALG .EQ. 2) GO TO 50
M = IV(SUSED)
WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
1 MODEL1(M), MODEL2(M), V(STPPAR)
GO TO 120
C
50 WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
1 V(STPPAR)
GO TO 120
C
C *** PRINT LONG SUMMARY LINE ***
C
60 IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 1) WRITE(PU,70)
70 FORMAT(/11H IT NF,6X,1HF,7X,5HRELDF,3X,6HPRELDF,3X,5HRELDX,
1 2X,13HMODEL STPPAR,2X,6HD*STEP,2X,7HNPRELDF)
IF (IV(NEEDHD) .EQ. 1 .AND. ALG .EQ. 2) WRITE(PU,80)
80 FORMAT(/11H IT NF,7X,1HF,8X,5HRELDF,4X,6HPRELDF,4X,5HRELDX,
1 3X,6HSTPPAR,3X,6HD*STEP,3X,7HNPRELDF)
IV(NEEDHD) = 0
NRELDF = ZERO
IF (OLDF .GT. ZERO) NRELDF = V(NREDUC) / OLDF
IF (ALG .EQ. 2) GO TO 90
M = IV(SUSED)
WRITE(PU,100) IV(NITER), NF, V(F), RELDF, PRELDF, V(RELDX),
1 MODEL1(M), MODEL2(M), V(STPPAR), V(DSTNRM), NRELDF
GO TO 120
C
90 WRITE(PU,110) IV(NITER), NF, V(F), RELDF, PRELDF,
1 V(RELDX), V(STPPAR), V(DSTNRM), NRELDF
100 FORMAT(I6,I5,D10.3,2D9.2,D8.1,A3,A4,2D8.1,D9.2)