This file is indexed.

/usr/share/gnudatalanguage/mpfit/tnmin.pro is in gdl-mpfit 1.85+2017.01.03-1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

   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
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
1291
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
1307
1308
1309
1310
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394
1395
1396
1397
1398
1399
1400
1401
1402
1403
1404
1405
1406
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421
1422
1423
1424
1425
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
1796
1797
1798
1799
1800
1801
1802
1803
1804
1805
1806
1807
1808
1809
1810
1811
1812
1813
1814
1815
1816
1817
1818
1819
1820
1821
1822
1823
1824
1825
1826
1827
1828
1829
1830
1831
1832
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843
1844
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
1871
1872
1873
1874
1875
1876
1877
1878
1879
1880
1881
1882
1883
1884
1885
1886
1887
1888
1889
1890
1891
1892
1893
1894
1895
1896
1897
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908
1909
1910
1911
1912
1913
1914
1915
1916
1917
1918
1919
1920
1921
1922
1923
1924
1925
1926
1927
1928
1929
1930
1931
1932
1933
1934
1935
1936
1937
1938
1939
1940
1941
1942
1943
1944
1945
1946
1947
1948
1949
1950
1951
1952
1953
1954
1955
1956
1957
1958
1959
1960
1961
1962
1963
1964
1965
1966
1967
1968
1969
1970
1971
1972
1973
1974
1975
1976
1977
1978
1979
1980
1981
1982
1983
1984
1985
1986
1987
1988
1989
1990
1991
1992
1993
1994
1995
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017
2018
2019
2020
2021
2022
2023
2024
2025
2026
2027
2028
2029
2030
2031
2032
2033
2034
2035
2036
2037
2038
2039
2040
2041
2042
2043
2044
2045
2046
2047
2048
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069
2070
2071
2072
2073
2074
2075
2076
2077
2078
2079
2080
2081
2082
2083
2084
2085
2086
2087
2088
2089
2090
2091
2092
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112
2113
2114
2115
2116
2117
2118
2119
2120
2121
2122
2123
2124
2125
2126
2127
2128
2129
2130
2131
2132
2133
2134
2135
2136
2137
2138
2139
2140
2141
2142
2143
2144
2145
2146
2147
2148
2149
2150
2151
2152
2153
2154
2155
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169
2170
2171
2172
2173
2174
2175
2176
2177
2178
2179
2180
2181
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195
2196
2197
2198
2199
2200
2201
2202
2203
2204
2205
2206
2207
2208
2209
2210
2211
2212
2213
2214
2215
2216
2217
2218
2219
2220
2221
2222
2223
2224
2225
2226
2227
2228
2229
2230
2231
2232
2233
2234
2235
2236
2237
2238
2239
2240
2241
2242
2243
2244
2245
2246
2247
2248
2249
2250
2251
2252
2253
2254
2255
2256
2257
2258
2259
2260
2261
2262
2263
2264
2265
2266
2267
2268
2269
2270
2271
2272
2273
;+
; NAME:
;   TNMIN
;
; AUTHOR:
;   Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
;   craigm@lheamail.gsfc.nasa.gov
;   UPDATED VERSIONs can be found on my WEB PAGE:
;      http://cow.physics.wisc.edu/~craigm/idl/idl.html
;
; PURPOSE:
;   Performs function minimization (Truncated-Newton Method)
;
; MAJOR TOPICS:
;   Optimization and Minimization
;
; CALLING SEQUENCE:
;   parms = TNMIN(MYFUNCT, X, FUNCTARGS=fcnargs, NFEV=nfev,
;                 MAXITER=maxiter, ERRMSG=errmsg, NPRINT=nprint,
;                 QUIET=quiet, XTOL=xtol, STATUS=status,
;                 FGUESS=fguess, PARINFO=parinfo, BESTMIN=bestmin,
;                 ITERPROC=iterproc, ITERARGS=iterargs, niter=niter)
;
; DESCRIPTION:
;
;  TNMIN uses the Truncated-Newton method to minimize an arbitrary IDL
;  function with respect to a given set of free parameters.  The
;  user-supplied function must compute the gradient with respect to
;  each parameter.  TNMIN is based on TN.F (TNBC) by Stephen Nash.
;
;  If you want to solve a least-squares problem, to perform *curve*
;  *fitting*, then you will probably want to use the routines MPFIT,
;  MPFITFUN and MPFITEXPR.  Those routines are specifically optimized
;  for the least-squares problem.  TNMIN is suitable for constrained
;  and unconstrained optimization problems with a medium number of
;  variables.  Function *maximization* can be performed using the
;  MAXIMIZE keyword.
;
;  TNMIN is similar to MPFIT in that it allows parameters to be fixed,
;  simple bounding limits to be placed on parameter values, and
;  parameters to be tied to other parameters.  One major difference
;  between MPFIT and TNMIN is that TNMIN does not compute derivatives
;  automatically by default.  See PARINFO and AUTODERIVATIVE below for
;  more discussion and examples.
;
; USER FUNCTION
;
;  The user must define an IDL function which returns the desired
;  value as a single scalar.  The IDL function must accept a list of
;  numerical parameters, P.  Additionally, keyword parameters may be
;  used to pass more data or information to the user function, via the
;  FUNCTARGS keyword.
;
;  The user function should be declared in the following way:
;
;     FUNCTION MYFUNCT, p, dp [, keywords permitted ]
;       ; Parameter values are passed in "p"
;       f  = ....   ; Compute function value
;       dp = ....   ; Compute partial derivatives (optional)
;       return, f
;     END
;
;  The function *must* accept at least one argument, the parameter
;  list P.  The vector P is implicitly assumed to be a one-dimensional
;  array.  Users may pass additional information to the function by
;  composing and _EXTRA structure and passing it in the FUNCTARGS
;  keyword.
;
;  User functions may also indicate a fatal error condition using the
;  ERROR_CODE common block variable, as described below under the
;  TNMIN_ERROR common block definition (by setting ERROR_CODE to a
;  number between -15 and -1).
;
;  EXPLICIT vs. NUMERICAL DERIVATIVES
;
;  By default, the user must compute the function gradient components
;  explicitly using AUTODERIVATIVE=0.  As explained below, numerical
;  derivatives can also be calculated using AUTODERIVATIVE=1.
;
;  For explicit derivatives, the user should call TNMIN using the
;  default keyword value AUTODERIVATIVE=0.  [ This is different
;  behavior from MPFIT, where AUTODERIVATIVE=1 is the default. ] The
;  IDL user routine should compute the gradient of the function as a
;  one-dimensional array of values, one for each of the parameters.
;  They are passed back to TNMIN via "dp" as shown above.
;
;  The derivatives with respect to fixed parameters are ignored; zero
;  is an appropriate value to insert for those derivatives.  Upon
;  input to the user function, DP is set to a vector with the same
;  length as P, with a value of 1 for a parameter which is free, and a
;  value of zero for a parameter which is fixed (and hence no
;  derivative needs to be calculated).  This input vector may be
;  overwritten as needed.
;
;  For numerical derivatives, a finite differencing approximation is
;  used to estimate the gradient values.  Users can activate this
;  feature by passing the keyword AUTODERIVATIVE=1.  Fine control over
;  this behavior can be achieved using the STEP, RELSTEP and TNSIDE
;  fields of the PARINFO structure.
;
;  When finite differencing is used for computing derivatives (ie,
;  when AUTODERIVATIVE=1), the parameter DP is not passed.  Therefore
;  functions can use N_PARAMS() to indicate whether they must compute
;  the derivatives or not.  However there is no penalty (other than
;  computation time) for computing the gradient values and users may
;  switch between AUTODERIVATIVE=0 or =1 in order to test both
;  scenarios.
;
; CONSTRAINING PARAMETER VALUES WITH THE PARINFO KEYWORD
;
;  The behavior of TNMIN can be modified with respect to each
;  parameter to be optimized.  A parameter value can be fixed; simple
;  boundary constraints can be imposed; limitations on the parameter
;  changes can be imposed; properties of the automatic derivative can
;  be modified; and parameters can be tied to one another.
;
;  These properties are governed by the PARINFO structure, which is
;  passed as a keyword parameter to TNMIN.
;
;  PARINFO should be an array of structures, one for each parameter.
;  Each parameter is associated with one element of the array, in
;  numerical order.  The structure can have the following entries
;  (none are required):
;
;     .VALUE - the starting parameter value (but see the START_PARAMS
;              parameter for more information).
;
;     .FIXED - a boolean value, whether the parameter is to be held
;              fixed or not.  Fixed parameters are not varied by
;              TNMIN, but are passed on to MYFUNCT for evaluation.
;
;     .LIMITED - a two-element boolean array.  If the first/second
;                element is set, then the parameter is bounded on the
;                lower/upper side.  A parameter can be bounded on both
;                sides.  Both LIMITED and LIMITS must be given
;                together.
;
;     .LIMITS - a two-element float or double array.  Gives the
;               parameter limits on the lower and upper sides,
;               respectively.  Zero, one or two of these values can be
;               set, depending on the values of LIMITED.  Both LIMITED
;               and LIMITS must be given together.
;
;     .PARNAME - a string, giving the name of the parameter.  The
;                fitting code of TNMIN does not use this tag in any
;                way.
;
;     .STEP - the step size to be used in calculating the numerical
;             derivatives.  If set to zero, then the step size is
;             computed automatically.  Ignored when AUTODERIVATIVE=0.
;
;     .TNSIDE - the sidedness of the finite difference when computing
;               numerical derivatives.  This field can take four
;               values:
;
;                  0 - one-sided derivative computed automatically
;                  1 - one-sided derivative (f(x+h) - f(x)  )/h
;                 -1 - one-sided derivative (f(x)   - f(x-h))/h
;                  2 - two-sided derivative (f(x+h) - f(x-h))/(2*h)
;
;              Where H is the STEP parameter described above.  The
;              "automatic" one-sided derivative method will chose a
;              direction for the finite difference which does not
;              violate any constraints.  The other methods do not
;              perform this check.  The two-sided method is in
;              principle more precise, but requires twice as many
;              function evaluations.  Default: 0.
;
;     .TIED - a string expression which "ties" the parameter to other
;             free or fixed parameters.  Any expression involving
;             constants and the parameter array P are permitted.
;             Example: if parameter 2 is always to be twice parameter
;             1 then use the following: parinfo(2).tied = '2 * P(1)'.
;             Since they are totally constrained, tied parameters are
;             considered to be fixed; no errors are computed for them.
;             [ NOTE: the PARNAME can't be used in expressions. ]
;
;  Future modifications to the PARINFO structure, if any, will involve
;  adding structure tags beginning with the two letters "MP" or "TN".
;  Therefore programmers are urged to avoid using tags starting with
;  these two combinations of letters; otherwise they are free to
;  include their own fields within the PARINFO structure, and they
;  will be ignored.
;
;  PARINFO Example:
;  parinfo = replicate({value:0.D, fixed:0, limited:[0,0], $
;                       limits:[0.D,0]}, 5)
;  parinfo(0).fixed = 1
;  parinfo(4).limited(0) = 1
;  parinfo(4).limits(0)  = 50.D
;  parinfo(*).value = [5.7D, 2.2, 500., 1.5, 2000.]
;
;  A total of 5 parameters, with starting values of 5.7,
;  2.2, 500, 1.5, and 2000 are given.  The first parameter
;  is fixed at a value of 5.7, and the last parameter is
;  constrained to be above 50.
;
;
; INPUTS:
;
;   MYFUNCT - a string variable containing the name of the function to
;             be minimized (see USER FUNCTION above).  The IDL routine
;             should return the value of the function and optionally
;             its gradients.
;
;   X - An array of starting values for each of the parameters of the
;       model.
;
;       This parameter is optional if the PARINFO keyword is used.
;       See above.  The PARINFO keyword provides a mechanism to fix or
;       constrain individual parameters.  If both X and PARINFO are
;       passed, then the starting *value* is taken from X, but the
;       *constraints* are taken from PARINFO.
;
;
; RETURNS:
;
;   Returns the array of best-fit parameters.
;
;
; KEYWORD PARAMETERS:
;
;   AUTODERIVATIVE - If this is set, derivatives of the function will
;                    be computed automatically via a finite
;                    differencing procedure.  If not set, then MYFUNCT
;                    must provide the (explicit) derivatives.
;                    Default: 0 (explicit derivatives required)
;
;   BESTMIN - upon return, the best minimum function value that TNMIN
;             could find.
;
;   EPSABS - a nonnegative real variable.  Termination occurs when the
;            absolute error between consecutive iterates is at most
;            EPSABS.  Note that using EPSREL is normally preferable
;            over EPSABS, except in cases where ABS(F) is much larger
;            than 1 at the optimal point.  A value of zero indicates
;            the absolute error test is not applied.  If EPSABS is
;            specified, then both EPSREL and EPSABS tests are applied;
;            if either succeeds then termination occurs.
;            Default: 0 (i.e., only EPSREL is applied).
;
;   EPSREL - a nonnegative input variable. Termination occurs when the
;            relative error between two consecutive iterates is at
;            most EPSREL.  Therefore, EPSREL measures the relative
;            error desired in the function.  An additional, more
;            lenient, stopping condition can be applied by specifying
;            the EPSABS keyword.
;            Default: 100 * Machine Precision
;
;   ERRMSG - a string error or warning message is returned.
;
;   FGUESS - provides the routine a guess to the minimum value.
;            Default: 0
;
;   FUNCTARGS - A structure which contains the parameters to be passed
;               to the user-supplied function specified by MYFUNCT via
;               the _EXTRA mechanism.  This is the way you can pass
;               additional data to your user-supplied function without
;               using common blocks.
;
;               Consider the following example:
;                if FUNCTARGS = { XVAL:[1.D,2,3], YVAL:[1.D,4,9]}
;                then the user supplied function should be declared
;                like this:
;                FUNCTION MYFUNCT, P, XVAL=x, YVAL=y
;
;               By default, no extra parameters are passed to the
;               user-supplied function.
;
;   ITERARGS - The keyword arguments to be passed to ITERPROC via the
;              _EXTRA mechanism.  This should be a structure, and is
;              similar in operation to FUNCTARGS.
;              Default: no arguments are passed.
;
;   ITERDERIV - Intended to print function gradient information.  If
;               set, then the ITERDERIV keyword is set in each call to
;               ITERPROC.  In the default ITERPROC, parameter values
;               and gradient information are both printed when this
;               keyword is set.
;
;   ITERPROC - The name of a procedure to be called upon each NPRINT
;              iteration of the TNMIN routine.  It should be declared
;              in the following way:
;
;              PRO ITERPROC, MYFUNCT, p, iter, fnorm, FUNCTARGS=fcnargs, $
;                PARINFO=parinfo, QUIET=quiet, _EXTRA=extra
;                ; perform custom iteration update
;              END
;
;              ITERPROC must accept the _EXTRA keyword, in case of
;              future changes to the calling procedure.
;
;              MYFUNCT is the user-supplied function to be minimized,
;              P is the current set of model parameters, ITER is the
;              iteration number, and FUNCTARGS are the arguments to be
;              passed to MYFUNCT.  FNORM is should be the function
;              value.  QUIET is set when no textual output should be
;              printed.  See below for documentation of PARINFO.
;
;              In implementation, ITERPROC can perform updates to the
;              terminal or graphical user interface, to provide
;              feedback while the fit proceeds.  If the fit is to be
;              stopped for any reason, then ITERPROC should set the
;              common block variable ERROR_CODE to negative value
;              between -15 and -1 (see TNMIN_ERROR common block
;              below).  In principle, ITERPROC should probably not
;              modify the parameter values, because it may interfere
;              with the algorithm's stability.  In practice it is
;              allowed.
;
;              Default: an internal routine is used to print the
;                       parameter values.
;
;   MAXITER - The maximum number of iterations to perform.  If the
;             number is exceeded, then the STATUS value is set to 5
;             and TNMIN returns.
;             Default: 200 iterations
;
;   MAXIMIZE - If set, the function is maximized instead of
;              minimized.
;
;   MAXNFEV - The maximum number of function evaluations to perform.
;             If the number is exceeded, then the STATUS value is set
;             to -17 and TNMIN returns.  A value of zero indicates no
;             maximum.
;             Default: 0 (no maximum)
;
;   NFEV - upon return, the number of MYFUNCT function evaluations
;          performed.
;
;   NITER - upon return, number of iterations completed.
;
;   NPRINT - The frequency with which ITERPROC is called.  A value of
;            1 indicates that ITERPROC is called with every iteration,
;            while 2 indicates every other iteration, etc.
;            Default value: 1
;
;   PARINFO - Provides a mechanism for more sophisticated constraints
;             to be placed on parameter values.  When PARINFO is not
;             passed, then it is assumed that all parameters are free
;             and unconstrained.  Values in PARINFO are never modified
;             during a call to TNMIN.
;
;             See description above for the structure of PARINFO.
;
;             Default value:  all parameters are free and unconstrained.
;
;   QUIET - set this keyword when no textual output should be printed
;           by TNMIN
;
;   STATUS - an integer status code is returned.  All values greater
;            than zero can represent success (however STATUS EQ 5 may
;            indicate failure to converge).  Gaps in the numbering
;            system below are to maintain compatibility with MPFIT.
;            Upon return, STATUS can have one of the following values:
;
;        -18  a fatal internal error occurred during optimization.
;
;        -17  the maximum number of function evaluations has been
;             reached without convergence.
;
;        -16  a parameter or function value has become infinite or an
;             undefined number.  This is usually a consequence of
;             numerical overflow in the user's function, which must be
;             avoided.
;
;        -15 to -1
;             these are error codes that either MYFUNCT or ITERPROC
;             may return to terminate the fitting process (see
;             description of TNMIN_ERROR common below).  If either
;             MYFUNCT or ITERPROC set ERROR_CODE to a negative number,
;             then that number is returned in STATUS.  Values from -15
;             to -1 are reserved for the user functions and will not
;             clash with TNMIN.
;
;	   0  improper input parameters.
;
;	   1  convergence was reached.
;
;          2-4 (RESERVED)
;
;	   5  the maximum number of iterations has been reached
;
;          6-8 (RESERVED)
;
;
; EXAMPLE:
;
;   FUNCTION F, X, DF, _EXTRA=extra  ;; *** MUST ACCEPT KEYWORDS
;     F  = (X(0)-1)^2 + (X(1)+7)^2
;     DF = [ 2D * (X(0)-1), 2D * (X(1)+7) ] ; Gradient
;     RETURN, F
;   END
;
;   P = TNMIN('F', [0D, 0D], BESTMIN=F0)
;   Minimizes the function F(x0,x1) = (x0-1)^2 + (x1+7)^2.
;
;
; COMMON BLOCKS:
;
;   COMMON TNMIN_ERROR, ERROR_CODE
;
;     User routines may stop the fitting process at any time by
;     setting an error condition.  This condition may be set in either
;     the user's model computation routine (MYFUNCT), or in the
;     iteration procedure (ITERPROC).
;
;     To stop the fitting, the above common block must be declared,
;     and ERROR_CODE must be set to a negative number.  After the user
;     procedure or function returns, TNMIN checks the value of this
;     common block variable and exits immediately if the error
;     condition has been set.  By default the value of ERROR_CODE is
;     zero, indicating a successful function/procedure call.
;
;
; REFERENCES:
;
;   TRUNCATED-NEWTON METHOD, TN.F
;      Stephen G. Nash, Operations Research and Applied Statistics
;      Department
;      http://www.netlib.org/opt/tn
;
;   Nash, S. G. 1984, "Newton-Type Minimization via the Lanczos
;      Method," SIAM J. Numerical Analysis, 21, p. 770-778
;
;
; MODIFICATION HISTORY:
;   Derived from TN.F by Stephen Nash with many changes and additions,
;      to conform to MPFIT paradigm, 19 Jan 1999, CM
;   Changed web address to COW, 25 Sep 1999, CM
;   Alphabetized documented keyword parameters, 02 Oct 1999, CM
;   Changed to ERROR_CODE for error condition, 28 Jan 2000, CM
;   Continued and fairly major improvements (CM, 08 Jan 2001):
;      - calling of user procedure is now concentrated in TNMIN_CALL,
;        and arguments are reduced by storing a large number of them
;        in common blocks;
;      - finite differencing is done within TNMIN_CALL; added
;        AUTODERIVATIVE=1 so that finite differencing can be enabled,
;        both one and two sided;
;      - a new procedure to parse PARINFO fields, borrowed from MPFIT;
;        brought PARINFO keywords up to date with MPFIT;
;      - go through and check for float vs. double discrepancies;
;      - add explicit MAXIMIZE keyword, and support in TNMIN_CALL and
;        TNMIN_DEFITER to print the correct values in that case;
;        TNMIN_DEFITER now prints function gradient values if
;        requested;
;      - convert to common-based system of MPFIT for storing machine
;        constants; revert TNMIN_ENORM to simple sum of squares, at
;        least for now;
;      - remove limit on number of function evaluations, at least for
;        now, and until I can understand what happens when we do
;        numerical derivatives.
;   Further changes: more float vs double; disable TNMINSTEP for now;
;     experimented with convergence test in case of function
;     maximization, 11 Jan 2001, CM
;   TNMINSTEP is parsed but not enabled, 13 Mar 2001
;   Major code cleanups; internal docs; reduced commons, CM, 08 Apr
;     2001
;   Continued code cleanups; documentation; the STATUS keyword
;     actually means something, CM, 10 Apr 2001
;   Added reference to Nash paper, CM, 08 Feb 2002
;   Fixed 16-bit loop indices, D. Schelgel, 22 Apr 2003
;   Changed parens to square brackets because of conflicts with
;     limits function.  K. Tolbert, 23 Feb 2005
;   Some documentation clarifications, CM, 09 Nov 2007
;   Ensure that MY_FUNCT returns a scalar; make it more likely that
;     error messages get back out to the user; fatal CATCH'd error 
;     now returns STATUS = -18, CM, 17 Sep 2008
;   Fix TNMIN_CALL when parameters are TIEd (thanks to Alfred de
;     Wijn), CM, 22 Nov 2009
;   Remember to TIE the parameters before final return (thanks to
;     Michael Smith), CM, 20 Jan 2010
;
; TODO
;  - scale derivatives semi-automatically;
;  - ability to scale and offset parameters;
;
;  $Id: tnmin.pro,v 1.20 2016/05/19 16:08:08 cmarkwar Exp $
;-
; Copyright (C) 1998-2001,2002,2003,2007,2008,2009 Craig Markwardt
; This software is provided as is without any warranty whatsoever.
; Permission to use, copy and distribute unmodified copies for
; non-commercial purposes, and to modify and use for personal or
; internal use, is granted.  All other rights are reserved.
;-

;%% TRUNCATED-NEWTON METHOD:  SUBROUTINES
;   FOR OTHER MACHINES, MODIFY ROUTINE MCHPR1 (MACHINE EPSILON)
;   WRITTEN BY:  STEPHEN G. NASH
;                OPERATIONS RESEARCH AND APPLIED STATISTICS DEPT.
;                GEORGE MASON UNIVERSITY
;                FAIRFAX, VA 22030
;******************************************************************

;; Routine which declares functions and common blocks
pro tnmin_dummy
  forward_function tnmin_enorm, tnmin_step1, tnmin
  forward_function tnmin_call, tnmin_autoder
  common tnmin_error, error_code
  common tnmin_machar, tnmin_machar_vals
  common tnmin_config, tnmin_tnconfig
  common tnmin_fcnargs, tnmin_tnfcnargs
  common tnmin_work, lsk, lyk, ldiagb, lsr, lyr
  a = 1
  return
end

;; Following are machine constants that can be loaded once.  I have
;; found that bizarre underflow messages can be produced in each call
;; to MACHAR(), so this structure minimizes the number of calls to
;; one.
pro tnmin_setmachar, double=isdouble

  common tnmin_machar, tnmin_machar_vals

  ;; In earlier versions of IDL, MACHAR itself could produce a load of
  ;; error messages.  We try to mask some of that out here.
  if (!version.release) LT 5 then dummy = check_math(1, 1)

  mch = 0.
  mch = machar(double=keyword_set(isdouble))
  dmachep = mch.eps
  dmaxnum = mch.xmax
  dminnum = mch.xmin
  dmaxlog = alog(mch.xmax)
  dminlog = alog(mch.xmin)
  if keyword_set(isdouble) then $
    dmaxgam = 171.624376956302725D $
  else $
    dmaxgam = 171.624376956302725
  drdwarf = sqrt(dminnum*1.5) * 10
  drgiant = sqrt(dmaxnum) * 0.1

  tnmin_machar_vals = {machep: dmachep, maxnum: dmaxnum, minnum: dminnum, $
                       maxlog: dmaxlog, minlog: dminlog, maxgam: dmaxgam, $
                       rdwarf: drdwarf, rgiant: drgiant}

  if (!version.release) LT 5 then dummy = check_math(0, 0)

  return
end

;; Procedure to parse the parameter values in PARINFO
pro tnmin_parinfo, parinfo, tnames, tag, values, default=def, status=status, $
                   n_param=n

  status = 0
  if n_elements(n) EQ 0 then n = n_elements(parinfo)

  if n EQ 0 then begin
      if n_elements(def) EQ 0 then return
      values = def
      status = 1
      return
  endif

  if n_elements(parinfo) EQ 0 then goto, DO_DEFAULT
  if n_elements(tnames) EQ 0 then tnames = tag_names(parinfo)
  wh = where(tnames EQ tag, ct)

  if ct EQ 0 then begin
      DO_DEFAULT:
      if n_elements(def) EQ 0 then return
      values = make_array(n, value=def(0))
      values(0) = def
  endif else begin
      values = parinfo.(wh(0))
  endelse

  status = 1
  return
end

;; Procedure to tie one parameter to another.
pro tnmin_tie, p, _ptied
  _wh = where(_ptied NE '', _ct)
 if _ct EQ 0 then return
  for _i = 0L, _ct-1 do begin
      _cmd = 'p('+strtrim(_wh(_i),2)+') = '+_ptied(_wh(_i))
      _err = execute(_cmd)
      if _err EQ 0 then begin
          message, 'ERROR: Tied expression "'+_cmd+'" failed.'
          return
      endif
  endfor
end

function tnmin_autoder, fcn, x, dx, dside=dside

  common tnmin_machar, machvals
  common tnmin_config, tnconfig

  MACHEP0 = machvals.machep
  DWARF   = machvals.minnum
  if n_elements(dside) NE n_elements(x) then dside = tnconfig.dside

  eps = sqrt(MACHEP0)
  h = eps * (1. + abs(x))

  ;; if STEP is given, use that
  wh = where(tnconfig.step GT 0, ct)
  if ct GT 0 then h(wh) = tnconfig.step(wh)

  ;; if relative step is given, use that
  wh = where(tnconfig.dstep GT 0, ct)
  if ct GT 0 then h(wh) = abs(tnconfig.dstep(wh)*x(wh))

  ;; In case any of the step values are zero
  wh = where(h EQ 0, ct)
  if ct GT 0 then h(wh) = eps

  ;; Reverse the sign of the step if we are up against the parameter
  ;; limit, or if the user requested it.
  mask = (dside EQ -1 OR (tnconfig.ulimited AND (x GT tnconfig.ulimit-h)))
  wh = where(mask, ct)
  if ct GT 0 then h(wh) = -h(wh)

  dx = x * 0.
  f = tnmin_call(fcn, x)
  for i = 0L, n_elements(x)-1 do begin
      if tnconfig.pfixed(i) EQ 1 then goto, NEXT_PAR
      hh = h(i)

      RESTART_PAR:
      xp = x
      xp(i) = xp(i) + hh

      fp = tnmin_call(fcn, xp)

      if abs(dside(i)) LE 1 then begin
          ;; COMPUTE THE ONE-SIDED DERIVATIVE
          dx(i) = (fp-f)/hh
      endif else begin
          ;; COMPUTE THE TWO-SIDED DERIVATIVE
          xp(i) = x(i) - hh

          fm = tnmin_call(fcn, xp)
          dx(i) = (fp-fm)/(2*hh)
      endelse
      NEXT_PAR:
  endfor

  return, f
end

;; Call user function or procedure, with _EXTRA or not, with
;; derivatives or not.
function tnmin_call, fcn, x1, dx, fullparam_=xall

;  on_error, 2
  common tnmin_config, tnconfig
  common tnmin_fcnargs, fcnargs

  ifree = tnconfig.ifree
  ;; Following promotes the byte array to a floating point array so
  ;; that users who simply re-fill the array aren't surprised when
  ;; their gradient comes out as bytes. :-)
  dx = tnconfig.pfixed + x1(0)*0.

  if n_elements(xall) GT 0 then begin
      x = xall
      x(ifree) = x1
  endif else begin
      x = x1
  endelse

  ;; Enforce TIEd parameters
  if keyword_set(tnconfig.qanytied) then tnmin_tie, x, tnconfig.ptied

  ;; Decide whether we are calling a procedure or function
  if tnconfig.proc then proc = 1 else proc = 0
  tnconfig.nfev = tnconfig.nfev + 1

  if n_params() EQ 3 then begin
      if tnconfig.autoderiv then $
        f = tnmin_autoder(fcn, x, dx) $
      else if n_elements(fcnargs) GT 0 then $
        f = call_function(fcn, x, dx, _EXTRA=fcnargs) $
      else $
        f = call_function(fcn, x, dx)

      dx = dx(ifree)
      if tnconfig.max then begin
          dx = -dx
          f = -f
      endif
  endif else begin
      if n_elements(fcnargs) GT 0 then $
        f = call_function(fcn, x, _EXTRA=fcnargs) $
      else $
        f = call_function(fcn, x)

      if n_elements(f) NE 1 then begin
          message, 'ERROR: function "'+fcn+'" returned a vector when '+$
            'a scalar was expected.'
      endif
  endelse

  if n_elements(f) GT 1 then return, f $
  else                       return, f(0)
end

function tnmin_enorm, vec

  common tnmin_config, tnconfig
  ;; Very simple-minded sum-of-squares
  if n_elements(tnconfig) GT 0 then if tnconfig.fastnorm then begin
      ans = sqrt(total(vec^2,1))
      goto, TERMINATE
  endif

  common tnmin_machar, machvals

  agiant = machvals.rgiant / n_elements(vec)
  adwarf = machvals.rdwarf * n_elements(vec)

  ;; This is hopefully a compromise between speed and robustness.
  ;; Need to do this because of the possibility of over- or underflow.
  mx = max(vec, min=mn)
  mx = max(abs([mx,mn]))
  if mx EQ 0 then return, vec(0)* 0.

  if mx GT agiant OR mx LT adwarf then ans = mx * sqrt(total((vec/mx)^2)) $
  else                                 ans = sqrt( total(vec^2) )

  TERMINATE:
  return, ans
end

;
; ROUTINES TO INITIALIZE PRECONDITIONER
;
pro tnmin_initpc, diagb, emat, n, upd1, yksk, gsk, yrsr, lreset
  ;; Rename common variables as they appear in INITP3.  Those
  ;; indicated in all caps are not used or renamed here.
; common tnmin_work, lsk, lyk, ldiagb, lsr, lyr
  common tnmin_work,  sk,  yk, LDIAGB,  sr,  yr
;                    I    I            I    I

  ;; From INITP3
  if keyword_set(upd1) then begin
      EMAT = DIAGB
  endif else if keyword_set(lreset) then begin
      BSK  = DIAGB*SK
      SDS  = TOTAL(SK*BSK)
      EMAT = DIAGB - DIAGB*DIAGB*SK*SK/SDS + YK*YK/YKSK
  endif else begin
      BSK  = DIAGB * SR
      SDS  = TOTAL(SR*BSK)
      SRDS = TOTAL(SK*BSK)
      YRSK = TOTAL(YR*SK)
      BSK  = DIAGB*SK - BSK*SRDS/SDS+YR*YRSK/YRSR
      EMAT = DIAGB-DIAGB*DIAGB*SR*SR/SDS+YR*YR/YRSR
      SDS  = TOTAL(SK*BSK)
      EMAT = EMAT - BSK*BSK/SDS+YK*YK/YKSK
  endelse

  return
end

pro tnmin_ssbfgs, n, gamma, sj, yj, hjv, hjyj, yjsj, yjhyj, $
         vsj, vhyj, hjp1v
;
; SELF-SCALED BFGS
;
  DELTA = (1. + GAMMA*YJHYJ/YJSJ)*VSJ/YJSJ - GAMMA*VHYJ/YJSJ
  BETA = -GAMMA*VSJ/YJSJ
  HJP1V = GAMMA*HJV + DELTA*SJ + BETA*HJYJ
  RETURN
end


;
; THIS ROUTINE ACTS AS A PRECONDITIONING STEP FOR THE
; LINEAR CONJUGATE-GRADIENT ROUTINE.  IT IS ALSO THE
; METHOD OF COMPUTING THE SEARCH DIRECTION FROM THE
; GRADIENT FOR THE NON-LINEAR CONJUGATE-GRADIENT CODE.
; IT REPRESENTS A TWO-STEP SELF-SCALED BFGS FORMULA.
;
pro tnmin_msolve, g, y, n, upd1, yksk, gsk, yrsr, lreset, first, $
                  hyr, hyk, ykhyk, yrhyr
  ;; Rename common variables as they appear in MSLV
; common tnmin_work, lsk, lyk, ldiagb, lsr, lyr
  common tnmin_work,  sk,  yk,  diagb,  sr,  yr
;                    I    I      I     I    I

  ;; From MSLV
  if keyword_set(UPD1) then begin
      Y = G / DIAGB
      RETURN
  endif

  ONE = G(0)*0 + 1.
  GSK = TOTAL(G*SK)

  if keyword_set(lreset) then begin
;
; COMPUTE GH AND HY WHERE H IS THE INVERSE OF THE DIAGONALS
;
      HG = G/DIAGB
      if keyword_set(FIRST) then begin
          HYK = YK/DIAGB
          YKHYK = TOTAL(YK*HYK)
      endif
      GHYK = TOTAL(G*HYK)
      TNMIN_SSBFGS,N,ONE,SK,YK,HG,HYK,YKSK, YKHYK,GSK,GHYK,Y
      RETURN
  endif

;
; COMPUTE HG AND HY WHERE H IS THE INVERSE OF THE DIAGONALS
;
  HG = G/DIAGB
  if keyword_set(FIRST) then begin
      HYK = YK/DIAGB
      HYR = YR/DIAGB
      YKSR = TOTAL(YK*SR)
      YKHYR = TOTAL(YK*HYR)
  endif
  GSR = TOTAL(G*SR)
  GHYR = TOTAL(G*HYR)
  if keyword_set(FIRST) then begin
      YRHYR = TOTAL(YR*HYR)
  endif

  TNMIN_SSBFGS,N,ONE,SR,YR,HG,HYR,YRSR, YRHYR,GSR,GHYR,HG
  if keyword_set(FIRST) then begin
      TNMIN_SSBFGS,N,ONE,SR,YR,HYK,HYR,YRSR, YRHYR,YKSR,YKHYR,HYK
  endif
  YKHYK = TOTAL(HYK*YK)
  GHYK  = TOTAL(HYK*G)
  TNMIN_SSBFGS,N,ONE,SK,YK,HG,HYK,YKSK, YKHYK,GSK,GHYK,Y

  RETURN
end

;
; THIS ROUTINE COMPUTES THE PRODUCT OF THE MATRIX G TIMES THE VECTOR
; V AND STORES THE RESULT IN THE VECTOR GV (FINITE-DIFFERENCE VERSION)
;
pro tnmin_gtims, v, gv, n, x, g, fcn, first, delta, accrcy, xnorm, $
        xnew

  IF keyword_set(FIRST) THEN BEGIN
      ;; Extra factor of ten is to avoid clashing with the finite
      ;; difference scheme which computes the derivatives
      DELTA = SQRT(100*ACCRCY)*(1.+XNORM)  ;; XXX diff than TN.F
;      DELTA = SQRT(ACCRCY)*(1.+XNORM)
      FIRST = 0
  ENDIF
  DINV = 1.  /DELTA

  F = tnmin_call(FCN, X + DELTA*V, GV, fullparam_=xnew)
  GV = (GV-G)*DINV
  return
end

;
; UPDATE THE PRECONDITIOING MATRIX BASED ON A DIAGONAL VERSION
; OF THE BFGS QUASI-NEWTON UPDATE.
;
pro tnmin_ndia3, n, e, v, gv, r, vgv
  VR = TOTAL(V*R)
  E = E - R*R/VR + GV*GV/VGV
  wh = where(e LE 1D-6, ct)
  if ct GT 0 then e(wh) = 1
  return
end

pro tnmin_fix, whlpeg, whupeg, z
  if whlpeg(0) NE -1 then z(whlpeg) = 0
  if whupeg(0) NE -1 then z(whupeg) = 0
end

;
; THIS ROUTINE PERFORMS A PRECONDITIONED CONJUGATE-GRADIENT
; ITERATION IN ORDER TO SOLVE THE NEWTON EQUATIONS FOR A SEARCH
; DIRECTION FOR A TRUNCATED-NEWTON ALGORITHM.  WHEN THE VALUE OF THE
; QUADRATIC MODEL IS SUFFICIENTLY REDUCED,
; THE ITERATION IS TERMINATED.
;
; PARAMETERS
;
; ZSOL        - COMPUTED SEARCH DIRECTION
; G           - CURRENT GRADIENT
; GV,GZ1,V    - SCRATCH VECTORS
; R           - RESIDUAL
; DIAGB,EMAT  - DIAGONAL PRECONDITONING MATRIX
; NITER       - NONLINEAR ITERATION #
; FEVAL       - VALUE OF QUADRATIC FUNCTION
pro tnmin_modlnp, zsol, gv, r, v, diagb, emat, $
         x, g, zk, n, niter, maxit, nmodif, nlincg, $
         upd1, yksk, gsk, yrsr, lreset, fcn, whlpeg, whupeg, $
         accrcy, gtp, gnorm, xnorm, xnew

;
; GENERAL INITIALIZATION
;
  zero = x(0)* 0.
  one = zero + 1
  IF (MAXIT EQ 0) THEN RETURN
  FIRST = 1
  RHSNRM = GNORM
  TOL = zero + 1.E-12
  QOLD = zero

;
; INITIALIZATION FOR PRECONDITIONED CONJUGATE-GRADIENT ALGORITHM
;
  tnmin_initpc, diagb, emat, n, upd1, yksk, gsk, yrsr, lreset

  R = -G
  V = G*0.
  ZSOL = V

;
; ************************************************************
; MAIN ITERATION
; ************************************************************
;
  FOR K = 1L, MAXIT DO BEGIN
      NLINCG = NLINCG + 1

;
; CG ITERATION TO SOLVE SYSTEM OF EQUATIONS
;
      tnmin_fix, whlpeg, whupeg, r
      TNMIN_MSOLVE, R, ZK, N, UPD1, YKSK, GSK, YRSR, LRESET, FIRST, $
        HYR, HYK, YKHYK, YRHYR
      tnmin_fix, whlpeg, whupeg, zk
      RZ = TOTAL(R*ZK)
      IF (RZ/RHSNRM LT TOL) THEN GOTO, MODLNP_80
      IF (K EQ 1) THEN BETA = ZERO $
      ELSE             BETA = RZ/RZOLD
      V = ZK + BETA*V
      tnmin_fix, whlpeg, whupeg, v
      TNMIN_GTIMS, V, GV, N, X, G, FCN, FIRST, DELTA, ACCRCY, XNORM, XNEW
      tnmin_fix, whlpeg, whupeg, gv
      VGV = TOTAL(V*GV)
      IF (VGV/RHSNRM LT TOL) THEN GOTO, MODLNP_50
      TNMIN_NDIA3, N,EMAT,V,GV,R,VGV
;
; COMPUTE LINEAR STEP LENGTH
;
      ALPHA = RZ / VGV
;
; COMPUTE CURRENT SOLUTION AND RELATED VECTORS
;
      ZSOL = ZSOL + ALPHA*V
      R = R - ALPHA*GV
;
; TEST FOR CONVERGENCE
;
      GTP = TOTAL(ZSOL*G)
      PR = TOTAL(R*ZSOL)
      QNEW = 5.E-1 * (GTP + PR)
      QTEST = K * (1.E0 - QOLD/QNEW)
      IF (QTEST LT 0.D0) THEN GOTO, MODLNP_70
      QOLD = QNEW
      IF (QTEST LE 5.D-1) THEN GOTO, MODLNP_70
;
; PERFORM CAUTIONARY TEST
;
      IF (GTP GT 0) THEN GOTO, MODLNP_40
      RZOLD = RZ
  ENDFOR

;
; TERMINATE ALGORITHM
;
  K = K-1
  goto, MODLNP_70

MODLNP_40:
  ZSOL = ZSOL - ALPHA*V
  GTP = TOTAL(ZSOL*G)
  goto, MODLNP_90

MODLNP_50:
  ;; printed output
MODLNP_60:
  IF (K GT 1) THEN GOTO, MODLNP_70
  TNMIN_MSOLVE,G,ZSOL,N,UPD1,YKSK,GSK,YRSR,LRESET,FIRST, $
    HYR, HYK, YKHYK, YRHYR
  ZSOL = -ZSOL
  tnmin_fix, whlpeg, whupeg, zsol
  GTP = TOTAL(ZSOL*G)
MODLNP_70:
  goto, MODLNP_90
MODLNP_80:
  IF (K  GT  1) THEN GOTO, MODLNP_70
  ZSOL = -G
  tnmin_fix, whlpeg, whupeg, zsol
  GTP = TOTAL(ZSOL*G)
  goto, MODLNP_70

;
; STORE (OR RESTORE) DIAGONAL PRECONDITIONING
;
MODLNP_90:
  diagb = emat
  return
end

function tnmin_step1, fnew, fm, gtp, smax, epsmch

; ********************************************************
; STEP1 RETURNS THE LENGTH OF THE INITIAL STEP TO BE TAKEN ALONG THE
; VECTOR P IN THE NEXT LINEAR SEARCH.
; ********************************************************

  D = ABS(FNEW-FM)
  ALPHA = FNEW(0)*0  + 1.
  IF (2.D0*D LE (-GTP) AND D GE EPSMCH) THEN $
    ALPHA = -2.*D/GTP
  IF (ALPHA GE SMAX) THEN ALPHA = SMAX
  return, alpha
end

;
; ************************************************************
; GETPTC, AN ALGORITHM FOR FINDING A STEPLENGTH, CALLED REPEATEDLY BY
; ROUTINES WHICH REQUIRE A STEP LENGTH TO BE COMPUTED USING CUBIC
; INTERPOLATION. THE PARAMETERS CONTAIN INFORMATION ABOUT THE INTERVAL
; IN WHICH A LOWER POINT IS TO BE FOUND AND FROM THIS GETPTC COMPUTES A
; POINT AT WHICH THE FUNCTION CAN BE EVALUATED BY THE CALLING PROGRAM.
; THE VALUE OF THE INTEGER PARAMETERS IENTRY DETERMINES THE PATH TAKEN
; THROUGH THE CODE.
; ************************************************************
pro tnmin_getptc, big, small, rtsmll, reltol, abstol, tnytol, $
         fpresn, eta, rmu, xbnd, u, fu, gu, xmin, fmin, gmin, $
         xw, fw, gw, a, b, oldf, b1, scxbnd, e, step, factor, $
         braktd, gtest1, gtest2, tol, ientry, itest

  ;; This chicanery is so that we get the data types right
  ZERO = fu(0)* 0.
; a1 = zero & scale = zero & chordm = zero
; chordu = zero & d1 = zero & d2 = zero
; denom = zero
  POINT1 = ZERO + 0.1
  HALF = ZERO + 0.5
  ONE = ZERO + 1
  THREE = ZERO + 3
  FIVE = ZERO + 5
  ELEVEN = ZERO + 11

  if ientry EQ 1 then begin ;; else clause = 20 (OK)
;
;      IENTRY=1
;      CHECK INPUT PARAMETERS
;
      ;; GETPTC_10:
      ITEST = 2
      IF (U LE ZERO OR XBND LE TNYTOL OR GU GT ZERO) THEN RETURN
      ITEST = 1
      IF (XBND LT ABSTOL) THEN ABSTOL = XBND
      TOL = ABSTOL
      TWOTOL = TOL + TOL
;
; A AND B DEFINE THE INTERVAL OF UNCERTAINTY, X AND XW ARE POINTS
; WITH LOWEST AND SECOND LOWEST FUNCTION VALUES SO FAR OBTAINED.
; INITIALIZE A,SMIN,XW AT ORIGIN AND CORRESPONDING VALUES OF
; FUNCTION AND PROJECTION OF THE GRADIENT ALONG DIRECTION OF SEARCH
; AT VALUES FOR LATEST ESTIMATE AT MINIMUM.
;
      A = ZERO
      XW = ZERO
      XMIN = ZERO
      OLDF = FU
      FMIN = FU
      FW = FU
      GW = GU
      GMIN = GU
      STEP = U
      FACTOR = FIVE
;
;      THE MINIMUM HAS NOT YET BEEN BRACKETED.
;
      BRAKTD = 0
;
; SET UP XBND AS A BOUND ON THE STEP TO BE TAKEN. (XBND IS NOT COMPUTED
; EXPLICITLY BUT SCXBND IS ITS SCALED VALUE.)  SET THE UPPER BOUND
; ON THE INTERVAL OF UNCERTAINTY INITIALLY TO XBND + TOL(XBND).
;
      SCXBND = XBND
      B = SCXBND + RELTOL*ABS(SCXBND) + ABSTOL
      E = B + B
      B1 = B
;
; COMPUTE THE CONSTANTS REQUIRED FOR THE TWO CONVERGENCE CRITERIA.
;
      GTEST1 = -RMU*GU
      GTEST2 = -ETA*GU
;
; SET IENTRY TO INDICATE THAT THIS IS THE FIRST ITERATION
;
      IENTRY = 2
      goto, GETPTC_210
  endif

;
; IENTRY = 2
;
; UPDATE A,B,XW, AND XMIN
;
  ;; GETPTC_20:
  IF (FU GT FMIN) THEN GOTO, GETPTC_60
;
; IF FUNCTION VALUE NOT INCREASED, NEW POINT BECOMES NEXT
; ORIGIN AND OTHER POINTS ARE SCALED ACCORDINGLY.
;
  CHORDU = OLDF - (XMIN + U)*GTEST1
  if NOT (FU LE CHORDU) then begin
;
; THE NEW FUNCTION VALUE DOES NOT SATISFY THE SUFFICIENT DECREASE
; CRITERION. PREPARE TO MOVE THE UPPER BOUND TO THIS POINT AND
; FORCE THE INTERPOLATION SCHEME TO EITHER BISECT THE INTERVAL OF
; UNCERTAINTY OR TAKE THE LINEAR INTERPOLATION STEP WHICH ESTIMATES
; THE ROOT OF F(ALPHA)=CHORD(ALPHA).
;
      CHORDM = OLDF - XMIN*GTEST1
      GU = -GMIN
      DENOM = CHORDM-FMIN
      IF (ABS(DENOM) LT 1.D-15) THEN BEGIN
          DENOM = ZERO + 1.E-15
          IF (CHORDM-FMIN LT 0.D0) THEN DENOM = -DENOM
      ENDIF
      IF (XMIN NE ZERO) THEN GU = GMIN*(CHORDU-FU)/DENOM
      FU = (HALF*U*(GMIN+GU) + FMIN) > FMIN
;
; IF FUNCTION VALUE INCREASED, ORIGIN REMAINS UNCHANGED
; BUT NEW POINT MAY NOW QUALIFY AS W.
;
      GETPTC_60:
      IF (U GE ZERO) THEN BEGIN
          B = U
          BRAKTD = 1
      ENDIF ELSE BEGIN
          A = U
      ENDELSE
      XW = U
      FW = FU
      GW = GU
  endif else begin
      ;; GETPTC_30:
      FW = FMIN
      FMIN = FU
      GW = GMIN
      GMIN = GU
      XMIN = XMIN + U
      A = A-U
      B = B-U
      XW = -U
      SCXBND = SCXBND - U
      IF (GU GT ZERO) THEN BEGIN
          B = ZERO
          BRAKTD = 1
      ENDIF ELSE BEGIN
          A = ZERO
      ENDELSE
      TOL = ABS(XMIN)*RELTOL + ABSTOL
  endelse

  TWOTOL = TOL + TOL
  XMIDPT = HALF*(A + B)

;
; CHECK TERMINATION CRITERIA
;
  CONVRG = ABS(XMIDPT) LE TWOTOL - HALF*(B-A) OR $
    ABS(GMIN) LE GTEST2 AND FMIN LT OLDF AND $
    (ABS(XMIN - XBND) GT TOL OR NOT BRAKTD)
  IF CONVRG THEN BEGIN
      ITEST = 0
      IF (XMIN NE ZERO) THEN RETURN
;
; IF THE FUNCTION HAS NOT BEEN REDUCED, CHECK TO SEE THAT THE RELATIVE
; CHANGE IN F(X) IS CONSISTENT WITH THE ESTIMATE OF THE DELTA-
; UNIMODALITY CONSTANT, TOL.  IF THE CHANGE IN F(X) IS LARGER THAN
; EXPECTED, REDUCE THE VALUE OF TOL.
;
      ITEST = 3
      IF (ABS(OLDF-FW) LE FPRESN*(ONE + ABS(OLDF))) THEN RETURN
      TOL = POINT1*TOL
      IF (TOL LT TNYTOL) THEN RETURN
      RELTOL = POINT1*RELTOL
      ABSTOL = POINT1*ABSTOL
      TWOTOL = POINT1*TWOTOL
  endif

;
; CONTINUE WITH THE COMPUTATION OF A TRIAL STEP LENGTH
;
  ;; GETPTC_100:
  R = ZERO
  Q = ZERO
  S = ZERO
  IF (ABS(E) GT TOL) THEN BEGIN
;
; FIT CUBIC THROUGH XMIN AND XW
;
      R = THREE*(FMIN-FW)/XW + GMIN + GW
      ABSR = ABS(R)
      Q = ABSR
      IF (GW EQ ZERO OR GMIN EQ ZERO) EQ 0 THEN BEGIN ;; else clause = 140 (OK)
;
; COMPUTE THE SQUARE ROOT OF (R*R - GMIN*GW) IN A WAY
; WHICH AVOIDS UNDERFLOW AND OVERFLOW.
;
          ABGW = ABS(GW)
          ABGMIN = ABS(GMIN)
          S = SQRT(ABGMIN)*SQRT(ABGW)
          IF ((GW/ABGW)*GMIN LE ZERO) THEN BEGIN
;
; COMPUTE THE SQUARE ROOT OF R*R + S*S.
;
              SUMSQ = ONE
              P = ZERO
              IF (ABSR LT S) THEN BEGIN ;; else clause = 110 (OK)
;
; THERE IS A POSSIBILITY OF OVERFLOW.
;
                  IF (S GT RTSMLL) THEN P = S*RTSMLL
                  IF (ABSR GE P) THEN SUMSQ = ONE +(ABSR/S)^2
                  SCALE = S
              endif else begin ;; flow to 120 (OK)
;
; THERE IS A POSSIBILITY OF UNDERFLOW.
;
                  ;; GETPTC_110:
                  IF (ABSR GT RTSMLL) THEN P = ABSR*RTSMLL
                  IF (S GE P) THEN SUMSQ = ONE + (S/ABSR)^2
                  SCALE = ABSR
              ENDELSE ;; flow to 120 (OK)
              ;; GETPTC_120:
              SUMSQ = SQRT(SUMSQ)
              Q = BIG
              IF (SCALE LT BIG/SUMSQ) THEN Q = SCALE*SUMSQ
          endif else begin ;; flow to 140
;
; COMPUTE THE SQUARE ROOT OF R*R - S*S
;
              ;; GETPTC_130:
              Q = SQRT(ABS(R+S))*SQRT(ABS(R-S))
              IF (R GE S OR R LE (-S)) EQ 0 THEN BEGIN
                  R = ZERO
                  Q = ZERO
                  goto, GETPTC_150
              endif
          endelse
      endif
;
; COMPUTE THE MINIMUM OF FITTED CUBIC
;
      ;; GETPTC_140:
      IF (XW LT ZERO) THEN Q = -Q
      S = XW*(GMIN - R - Q)
      Q = GW - GMIN + Q + Q
      IF (Q GT ZERO) THEN S = -S
      IF (Q LE ZERO) THEN Q = -Q
      R = E
      IF (B1 NE STEP OR BRAKTD) THEN E = STEP
  endif

;
; CONSTRUCT AN ARTIFICIAL BOUND ON THE ESTIMATED STEPLENGTH
;
GETPTC_150:
  A1 = A
  B1 = B
  STEP = XMIDPT
  IF (BRAKTD) EQ 0 THEN BEGIN ;; else flow to 160 (OK)
      STEP = -FACTOR*XW
      IF (STEP GT SCXBND) THEN STEP = SCXBND
      IF (STEP NE SCXBND) THEN FACTOR = FIVE*FACTOR
      ;; flow to 170 (OK)
  endif else begin
;
; IF THE MINIMUM IS BRACKETED BY 0 AND XW THE STEP MUST LIE
; WITHIN (A,B).
;
      ;; GETPTC_160:
      if (a NE zero OR xw GE zero) AND (b NE zero OR xw LE zero) then $
        goto, GETPTC_180
;
; IF THE MINIMUM IS NOT BRACKETED BY 0 AND XW THE STEP MUST LIE
; WITHIN (A1,B1).
;
      D1 = XW
      D2 = A
      IF (A EQ ZERO) THEN D2 = B
; THIS LINE MIGHT BE
;     IF (A EQ ZERO) THEN D2 = E
      U = - D1/D2
      STEP = FIVE*D2*(POINT1 + ONE/U)/ELEVEN
      IF (U LT ONE) THEN STEP = HALF*D2*SQRT(U)
  endelse
  ;; GETPTC_170:
  IF (STEP LE ZERO) THEN A1 = STEP
  IF (STEP GT ZERO) THEN B1 = STEP
;
; REJECT THE STEP OBTAINED BY INTERPOLATION IF IT LIES OUTSIDE THE
; REQUIRED INTERVAL OR IT IS GREATER THAN HALF THE STEP OBTAINED
; DURING THE LAST-BUT-ONE ITERATION.
;
GETPTC_180:
  if NOT (abs(s) LE abs(half*q*r) OR s LE q*a1 OR s GE q*b1) then begin
      ;; else clause = 200 (OK)
;
; A CUBIC INTERPOLATION STEP
;
      STEP = S/Q
;
; THE FUNCTION MUST NOT BE EVALUTATED TOO CLOSE TO A OR B.
;
      if NOT (step - a GE twotol AND b - step GE twotol) then begin
          ;; else clause = 210 (OK)
          IF (XMIDPT LE ZERO) THEN STEP = -TOL ELSE STEP = TOL
      endif ;; flow to 210 (OK)
  endif else begin
      ;; GETPTC_200:
      E = B-A
  endelse

;
; IF THE STEP IS TOO LARGE, REPLACE BY THE SCALED BOUND (SO AS TO
; COMPUTE THE NEW POINT ON THE BOUNDARY).
;
GETPTC_210:
  if (step GE scxbnd) then begin ;; else clause = 220 (OK)
      STEP = SCXBND
;
; MOVE SXBD TO THE LEFT SO THAT SBND + TOL(XBND) = XBND.
;
      SCXBND = SCXBND - (RELTOL*ABS(XBND)+ABSTOL)/(ONE + RELTOL)
  endif
  ;; GETPTC_220:
  U = STEP
  IF (ABS(STEP) LT TOL AND STEP LT ZERO) THEN U = -TOL
  IF (ABS(STEP) LT TOL AND STEP GE ZERO) THEN U = TOL
  ITEST = 1
  RETURN
end

;
;      LINE SEARCH ALGORITHMS OF GILL AND MURRAY
;
pro tnmin_linder, n, fcn, small, epsmch, reltol, abstol, $
         tnytol, eta, sftbnd, xbnd, p, gtp, x, f, alpha, g, $
         iflag, xnew

  zero = f(0) * 0.
  one = zero + 1.

  LSPRNT = 0L
  NPRNT  = 10000L
  RTSMLL = SQRT(SMALL)
  BIG = 1./SMALL
  ITCNT = 0L

;
;     SET THE ESTIMATED RELATIVE PRECISION IN F(X).
;
  FPRESN = 10.*EPSMCH
  U = ALPHA
  FU = F
  FMIN = F
  GU = GTP
  RMU = zero + 1E-4

;
;      FIRST ENTRY SETS UP THE INITIAL INTERVAL OF UNCERTAINTY.
;
  IENTRY = 1L

LINDER_10:
;
; TEST FOR TOO MANY ITERATIONS
;
  ITCNT = ITCNT + 1
  IF (ITCNT GT 30) THEN BEGIN
      ;; deviation from Nash: allow optimization to continue in outer
      ;; loop even if we fail to converge, if IFLAG EQ 0.  A value of
      ;; 1 indicates failure.  I believe that I tried IFLAG=0 once and
      ;; there was some problem, but I forget what it was.
      IFLAG = 1
      F = FMIN
      ALPHA = XMIN
      X = X + ALPHA*P
      RETURN
  ENDIF

  IFLAG = 0
  TNMIN_GETPTC,BIG,SMALL,RTSMLL,RELTOL,ABSTOL,TNYTOL, $
    FPRESN,ETA,RMU,XBND,U,FU,GU,XMIN,FMIN,GMIN, $
    XW,FW,GW,A,B,OLDF,B1,SCXBND,E,STEP,FACTOR, $
    BRAKTD,GTEST1,GTEST2,TOL,IENTRY,ITEST

;
;      IF ITEST=1, THE ALGORITHM REQUIRES THE FUNCTION VALUE TO BE
;      CALCULATED.
;
  IF (ITEST EQ 1) THEN BEGIN
      UALPHA = XMIN + U
      FU = TNMIN_CALL(FCN, X + UALPHA*P, LG, fullparam_=xnew)
      GU = TOTAL(LG*P)
;
;      THE GRADIENT VECTOR CORRESPONDING TO THE BEST POINT IS
;      OVERWRITTEN IF FU IS LESS THAN FMIN AND FU IS SUFFICIENTLY
;      LOWER THAN F AT THE ORIGIN.
;
      IF (FU LE FMIN AND FU LE OLDF-UALPHA*GTEST1) THEN $
        G = LG
;      print, 'fu = ', fu
      GOTO, LINDER_10
  ENDIF
;
;      IF ITEST=2 OR 3 A LOWER POINT COULD NOT BE FOUND
;
  IFLAG = 1
  IF (ITEST NE 0) THEN RETURN

;
;      IF ITEST=0 A SUCCESSFUL SEARCH HAS BEEN MADE
;
;  print, 'itcnt = ', itcnt
  IFLAG = 0
  F = FMIN
  ALPHA = XMIN
  X = X + ALPHA*P
  RETURN
END

pro tnmin_defiter, fcn, x, iter, fnorm, fmt=fmt, FUNCTARGS=fcnargs, $
                   quiet=quiet, deriv=df, dprint=dprint, pfixed=pfixed, $
                   maximize=maximize, _EXTRA=iterargs

  if keyword_set(quiet) then return
  if n_params() EQ 3 then begin
      fnorm = tnmin_call(fcn, x, df)
  endif

  if keyword_set(maximize) then f = -fnorm else f = fnorm
  print, iter, f, format='("Iter ",I6,"   FUNCTION = ",G20.8)'
  if n_elements(fmt) GT 0 then begin
      print, x, format=fmt
  endif else begin
      n = n_elements(x)
      ii = lindgen(n)
      p = '     P('+strtrim(ii,2)+') = '+string(x,format='(G)')
      if keyword_set(dprint) then begin
          p1 = strarr(n)
          wh = where(pfixed EQ 0, ct)
          if ct GT 0 AND n_elements(df) GE ct then begin
              if keyword_set(maximize) then df1 = -df else df1 = df
              p1(wh) = string(df1, format='(G)')
          endif
          wh = where(pfixed EQ 1, ct)
          if ct GT 0 then $
            p1(wh) = '          (FIXED)'
          p = p + '  :  DF/DP('+strtrim(ii,2)+') = '+p1
      endif
      print, p, format='(A)'
  endelse

  return
end

;      SUBROUTINE TNBC (IERROR, N, X, F, G, W, LW, SFUN, LOW, UP, IPIVOT)
;      IMPLICIT          DOUBLE PRECISION (A-H,O-Z)
;      INTEGER           IERROR, N, LW, IPIVOT(N)
;      DOUBLE PRECISION  X(N), G(N), F, W(LW), LOW(N), UP(N)
;
; THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM
;
;   MINIMIZE     F(X)
;      X
;   SUBJECT TO   LOW <= X <= UP
;
; WHERE X IS A VECTOR OF N REAL VARIABLES.  THE METHOD USED IS
; A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA
; THE LANCZOS ALGORITHM" BY S.G. NASH (TECHNICAL REPORT 378, MATH.
; THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984),
; PP. 770-778).  THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X).  IT DOES
; NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A
; GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW.
; IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS
; ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE.
;
; SUBROUTINE PARAMETERS:
;
; IERROR  - (INTEGER) ERROR CODE
;           ( 0 => NORMAL RETURN
;           ( 2 => MORE THAN MAXFUN EVALUATIONS
;           ( 3 => LINE SEARCH FAILED TO FIND LOWER
;           (          POINT (MAY NOT BE SERIOUS)
;           (-1 => ERROR IN INPUT PARAMETERS
; N       - (INTEGER) NUMBER OF VARIABLES
; X       - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL
;           ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION.
; G       - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL
;           VALUE OF THE GRADIENT
; F       - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE
;           OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE
;           OF THE OBJECTIVE FUNCTION AT THE SOLUTION
; W       - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N
; LW      - (INTEGER) THE DECLARED DIMENSION OF W
; SFUN    - A USER-SPECIFIED SUBROUTINE THAT COMPUTES THE FUNCTION
;           AND GRADIENT OF THE OBJECTIVE FUNCTION.  IT MUST HAVE
;           THE CALLING SEQUENCE
;             SUBROUTINE SFUN (N, X, F, G)
;             INTEGER           N
;             DOUBLE PRECISION  X(N), G(N), F
; LOW, UP - (REAL*8) VECTORS OF LENGTH AT LEAST N CONTAINING
;           THE LOWER AND UPPER BOUNDS ON THE VARIABLES.  IF
;           THERE ARE NO BOUNDS ON A PARTICULAR VARIABLE, SET
;           THE BOUNDS TO -1.D38 AND 1.D38, RESPECTIVELY.
; IPIVOT  - (INTEGER) WORK VECTOR OF LENGTH AT LEAST N, USED
;           TO RECORD WHICH VARIABLES ARE AT THEIR BOUNDS.
;
; THIS IS AN EASY-TO-USE DRIVER FOR THE MAIN OPTIMIZATION ROUTINE
; LMQNBC.  MORE EXPERIENCED USERS WHO WISH TO CUSTOMIZE PERFORMANCE
; OF THIS ALGORITHM SHOULD CALL LMQBC DIRECTLY.
;
;----------------------------------------------------------------------
; THIS ROUTINE SETS UP ALL THE PARAMETERS FOR THE TRUNCATED-NEWTON
; ALGORITHM.  THE PARAMETERS ARE:
;
; ETA    - SEVERITY OF THE LINESEARCH
; MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS
; XTOL   - DESIRED ACCURACY FOR THE SOLUTION X*
; STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH
; ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES
; MSGLVL - CONTROLS QUANTITY OF PRINTED OUTPUT
;          0 = NONE, 1 = ONE LINE PER MAJOR ITERATION.
; MAXIT  - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP
;
function tnmin, fcn, xall, fguess=fguess, functargs=fcnargs, parinfo=parinfo, $
                epsrel=epsrel0, epsabs=epsabs0, fastnorm=fastnorm, $
                nfev=nfev, maxiter=maxiter0, maxnfev=maxfun0, maximize=fmax, $
                errmsg=errmsg, nprint=nprint, status=status, nocatch=nocatch, $
                iterproc=iterproc, iterargs=iterargs, niter=niter,quiet=quiet,$
                autoderivative=autoderiv, iterderiv=iterderiv, bestmin=f

  if n_elements(nprint) EQ 0 then nprint = 1
  if n_elements(iterproc) EQ 0 then iterproc = 'TNMIN_DEFITER'
  if n_elements(autoderiv) EQ 0 then autoderiv = 0
  if n_elements(msglvl) EQ 0 then msglvl = 0
  if n_params() EQ 0 then begin
      message, "USAGE: PARMS = TNMIN('MYFUNCT', START_PARAMS, ... )", /info
      return, !values.d_nan
  endif
  iterd = keyword_set(iterderiv)
  maximize = keyword_set(fmax)
  status = 0L
  nfev = 0L
  errmsg = ''
  catch_msg = 'in TNMIN'

  common tnmin_config, tnconfig
  tnconfig = {fastnorm: keyword_set(fastnorm), proc: 0, nfev: 0L, $
              autoderiv: keyword_set(autoderiv), max: maximize}

  ;; Handle error conditions gracefully
  if NOT keyword_set(nocatch) then begin
      catch, catcherror
      if catcherror NE 0 then begin
          catch, /cancel
          err_string = ''+!error_state.msg
          message, /cont, 'Error detected while '+catch_msg+':'
          message, /cont,    err_string
          message, /cont, 'Error condition detected. Returning to MAIN level.'
          if errmsg EQ '' then $
            errmsg = 'Error detected while '+catch_msg+': '+err_string
          if status EQ 0 then status = -18
          return, !values.d_nan
      endif
  endif

  ;; Parinfo:
  ;; --------------- STARTING/CONFIG INFO (passed in to routine, not changed)
  ;; .value   - starting value for parameter
  ;; .fixed   - parameter is fixed
  ;; .limited - a two-element array, if parameter is bounded on
  ;;            lower/upper side
  ;; .limits - a two-element array, lower/upper parameter bounds, if
  ;;           limited vale is set
  ;; .step   - step size in Jacobian calc, if greater than zero

  catch_msg = 'parsing input parameters'
  ;; Parameters can either be stored in parinfo, or x.  Parinfo takes
  ;; precedence if it exists.
  if n_elements(xall) EQ 0 AND n_elements(parinfo) EQ 0 then begin
      errmsg = 'ERROR: must pass parameters in X or PARINFO'
      goto, TERMINATE
  endif

  ;; Be sure that PARINFO is of the right type
  if n_elements(parinfo) GT 0 then begin
      parinfo_size = size(parinfo)
      if parinfo_size(parinfo_size(0)+1) NE 8 then begin
          errmsg = 'ERROR: PARINFO must be a structure.'
          goto, TERMINATE
      endif
      if n_elements(xall) GT 0 AND n_elements(xall) NE n_elements(parinfo) $
        then begin
          errmsg = 'ERROR: number of elements in PARINFO and X must agree'
          goto, TERMINATE
      endif
  endif

  ;; If the parameters were not specified at the command line, then
  ;; extract them from PARINFO
  if n_elements(xall) EQ 0 then begin
      tnmin_parinfo, parinfo, tagnames, 'VALUE', xall, status=stx
      if stx EQ 0 then begin
          errmsg = 'ERROR: either X or PARINFO(*).VALUE must be supplied.'
          goto, TERMINATE
      endif

      sz = size(xall)
      ;; Convert to double if parameters are not float or double
      if sz(sz(0)+1) NE 4 AND sz(sz(0)+1) NE 5 then $
        xall = double(xall)
  endif
  npar = n_elements(xall)
  zero = xall(0) * 0.
  one  = zero + 1
  ten  = zero + 10
  twothird = (zero+2)/(zero+3)
  quarter = zero + 0.25
  half = zero + 0.5

  ;; Extract machine parameters
  sz = size(xall)
  tp = sz(sz(0)+1)
  if tp NE 4 AND tp NE 5 then begin
      if NOT keyword_set(quiet) then begin
          message, 'WARNING: input parameters must be at least FLOAT', /info
          message, '         (converting parameters to FLOAT)', /info
      endif
      xall = float(xall)
      sz = size(xall)
  endif
  isdouble = (sz(sz(0)+1) EQ 5)

  common tnmin_machar, machvals
  tnmin_setmachar, double=isdouble
  MCHPR1 = machvals.machep

  ;; TIED parameters?
  tnmin_parinfo, parinfo, tagnames, 'TIED', ptied, default='', n=npar
  ptied = strtrim(ptied, 2)
  wh = where(ptied NE '', qanytied)
  qanytied = qanytied GT 0
  tnconfig = create_struct(tnconfig, 'QANYTIED', qanytied, 'PTIED', ptied)

  ;; FIXED parameters ?
  tnmin_parinfo, parinfo, tagnames, 'FIXED', pfixed, default=0, n=npar
  pfixed = pfixed EQ 1
  pfixed = pfixed OR (ptied NE '')   ;; Tied parameters are also effectively fixed

  ;; Finite differencing step, absolute and relative, and sidedness of derivative
  tnmin_parinfo, parinfo, tagnames, 'STEP',     step, default=zero, n=npar
  tnmin_parinfo, parinfo, tagnames, 'RELSTEP', dstep, default=zero, n=npar
  tnmin_parinfo, parinfo, tagnames, 'TNSIDE',  dside, default=2,    n=npar

  ;; Maximum and minimum steps allowed to be taken in one iteration
  tnmin_parinfo, parinfo, tagnames, 'TNMAXSTEP', maxstep, default=zero, n=npar
  tnmin_parinfo, parinfo, tagnames, 'TNMINSTEP', minstep, default=zero, n=npar
  qmin = minstep *  0 ;; Disable minstep for now
  qmax = maxstep NE 0
  wh = where(qmin AND qmax AND maxstep LT minstep, ct)
  if ct GT 0 then begin
      errmsg = 'ERROR: TNMINSTEP is greater than TNMAXSTEP'
      goto, TERMINATE
  endif
  wh = where(qmin AND qmax, ct)
  qminmax = ct GT 0

  ;; Finish up the free parameters
  ifree = where(pfixed NE 1, ct)
  if ct EQ 0 then begin
      errmsg = 'ERROR: no free parameters'
      goto, TERMINATE
  endif

  ;; Compose only VARYING parameters
  xnew = xall      ;; xnew is the set of parameters to be returned
  x = xnew(ifree)  ;; x is the set of free parameters

  ;; LIMITED parameters ?
  tnmin_parinfo, parinfo, tagnames, 'LIMITED', limited, status=st1
  tnmin_parinfo, parinfo, tagnames, 'LIMITS',  limits,  status=st2
  if st1 EQ 1 AND st2 EQ 1 then begin

      ;; Error checking on limits in parinfo
      wh = where((limited[0,*] AND xall LT limits[0,*]) OR $
                 (limited[1,*] AND xall GT limits[1,*]), ct)
      if ct GT 0 then begin
          errmsg = 'ERROR: parameters are not within PARINFO limits'
          goto, TERMINATE
      endif
      wh = where(limited[0,*] AND limited[1,*] AND $
                 limits[0,*] GE limits[1,*] AND pfixed EQ 0, ct)
      if ct GT 0 then begin
          errmsg = 'ERROR: PARINFO parameter limits are not consistent'
          goto, TERMINATE
      endif


      ;; Transfer structure values to local variables
      qulim = limited[1, ifree]
      ulim  = limits [1, ifree]
      qllim = limited[0, ifree]
      llim  = limits [0, ifree]

      wh = where(qulim OR qllim, ct)
      if ct GT 0 then qanylim = 1 else qanylim = 0

  endif else begin

      ;; Fill in local variables with dummy values
      qulim = lonarr(n_elements(ifree))
      ulim  = x * 0.
      qllim = qulim
      llim  = x * 0.
      qanylim = 0

  endelse

  tnconfig = create_struct(tnconfig, $
                           'PFIXED', pfixed, 'IFREE', ifree, $
                           'STEP', step, 'DSTEP', dstep, 'DSIDE', dside, $
                           'ULIMITED', qulim, 'ULIMIT', ulim)


  common tnmin_fcnargs, tnfcnargs
  tnfcnargs = 0 & dummy = temporary(tnfcnargs)
  if n_elements(fcnargs) GT 0 then tnfcnargs = fcnargs

  ;; SET UP CUSTOMIZING PARAMETERS
  ;; ETA    - SEVERITY OF THE LINESEARCH
  ;; MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS
  ;; XTOL   - DESIRED ACCURACY FOR THE SOLUTION X*
  ;; STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH
  ;; ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES
  ;; MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT
  ;;          0 = NONE, 1 = ONE LINE PER MAJOR ITERATION.
  ;; MAXIT  - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP

  n = n_elements(x)
  if n_elements(maxit) EQ 0 then begin
      maxit = (n/2) < 50 > 2    ;; XXX diff than TN.F
  endif

  if n_elements(maxfun0) EQ 0 then $
    maxfun = 0L $
  else $
    maxfun = floor(maxfun0(0)) > 1
;  maxfun = 150L*n
;  if keyword_set(autoderiv) then $
;    maxfun = maxfun*(1L + round(total(abs(dside)> 1 < 2)))
  eta = zero + 0.25
  stepmx = zero + 10

  if n_elements(maxiter0) EQ 0 then $
    maxiter = 200L $
  else $
    maxiter = floor(maxiter0(0)) > 1

  g = replicate(x(0)* 0., n)

;; call minimizer
;
; THIS ROUTINE IS A BOUNDS-CONSTRAINED TRUNCATED-NEWTON METHOD.
; THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY
; QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED
; IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3).
; FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TNBC.
;

;
; initialize variables
;
  common tnmin_work, lsk, lyk, ldiagb, lsr, lyr
;                    I/O  I/O     I/O  I/O  I/O
  lsk = 0 & lyk = 0 & ldiagb = 0 & lsr = 0 & lyr = 0

  zero   = x(0)* 0.
  one    = zero + 1

  if n_elements(fguess) EQ 0 then fguess = one
  if maximize then f = -fguess else f = fguess
  conv = 0 & lreset = 0 & upd1 = 0 & newcon = 0
  gsk = zero & yksk = zero & gtp = zero & gtpnew = zero & yrsr = zero

  upd1 = 1
  ireset = 0L
  nmodif = 0L
  nlincg = 0L
  fstop  = f
  conv   = 0
  nm1    = n - 1

;; From CHKUCP
;
; CHECKS PARAMETERS AND SETS CONSTANTS WHICH ARE COMMON TO BOTH
; DERIVATIVE AND NON-DERIVATIVE ALGORITHMS
;
  EPSMCH = MCHPR1
  SMALL = EPSMCH*EPSMCH
  TINY = SMALL
  NWHY = -1L
;
; SET CONSTANTS FOR LATER
;
  ;; Some of these constants have been moved around for clarity (!)
  if n_elements(epsrel0) EQ 0 then epsrel = 100*MCHPR1 $
  else                             epsrel = epsrel0(0)+0.
  if n_elements(epsabs0) EQ 0 then epsabs = zero $
  else                             epsabs = abs(epsabs0(0))+0.

  ACCRCY = epsrel
  XTOL   = sqrt(ACCRCY)

  RTEPS = SQRT(EPSMCH)
  RTOL = XTOL
  IF (ABS(RTOL) LT ACCRCY) THEN RTOL = 10. *RTEPS

  FTOL2 = 0
  FTOL1 = RTOL^2 + EPSMCH    ;; For func chg convergence test (U1a)
  if epsabs NE 0 then $
    FTOL2 = EPSABS + EPSMCH ;; For absolute func convergence test (U1b)
  PTOL = RTOL + RTEPS       ;; For parm chg convergence test (U2)
  GTOL1 = ACCRCY^TWOTHIRD   ;; For gradient convergence test (U3, squared)
  GTOL2 = (1D-2*XTOL)^2     ;; For gradient convergence test (U4, squared)

;
; CHECK FOR ERRORS IN THE INPUT PARAMETERS
;
  IF (ETA LT 0.D0 OR STEPMX LT RTOL) THEN BEGIN
      errmsg = 'ERROR: input keywords are inconsistent'
      goto, TERMINATE
  endif
  ;; Check input parameters for errors
  if (n LE 0) OR (xtol LE 0) OR (maxit LE 0) then begin
      errmsg = 'ERROR: input keywords are inconsistent'
      goto, TERMINATE
  endif
  NWHY = 0L

  XNORM = TNMIN_ENORM(X)
  ALPHA = zero
  TEST = zero

; From SETUCR
;
; CHECK INPUT PARAMETERS, COMPUTE THE INITIAL FUNCTION VALUE, SET
; CONSTANTS FOR THE SUBSEQUENT MINIMIZATION
;
  fm = f
;
; COMPUTE THE INITIAL FUNCTION VALUE
;
  catch_msg = 'calling TNMIN_CALL'
  fnew = tnmin_call(fcn, x, g, fullparam_=xnew)

;
; SET CONSTANTS FOR LATER
;
  NITER = 0L
  OLDF = FNEW
  GTG = TOTAL(G*G)

  common tnmin_error, tnerr

  if nprint GT 0 AND iterproc NE '' then begin
      iflag = 0L
      if (niter-1) MOD nprint EQ 0 then begin
          catch_msg = 'calling '+iterproc
          tnerr = 0
          call_procedure, iterproc, fcn, xnew, niter, fnew, $
            FUNCTARGS=fcnargs, parinfo=parinfo, quiet=quiet, $
            dprint=iterd, deriv=g, pfixed=pfixed, maximize=maximize, $
            _EXTRA=iterargs
          iflag = tnerr
          if iflag LT 0 then begin
              errmsg = 'WARNING: premature termination by "'+iterproc+'"'
              nwhy = 4L
              goto, CLEANUP
          endif
      endif
  endif


  fold = fnew
  flast = fnew

; From LMQNBC
;
; TEST THE LAGRANGE MULTIPLIERS TO SEE IF THEY ARE NON-NEGATIVE.
; BECAUSE THE CONSTRAINTS ARE ONLY LOWER BOUNDS, THE COMPONENTS
; OF THE GRADIENT CORRESPONDING TO THE ACTIVE CONSTRAINTS ARE THE
; LAGRANGE MULTIPLIERS.  AFTERWORDS, THE PROJECTED GRADIENT IS FORMED.
;
  catch_msg = 'zeroing derivatives of pegged parameters'
  lmask = qllim AND (x EQ llim) AND (g GE 0)
  umask = qulim AND (x EQ ulim) AND (g LE 0)
  whlpeg = where(lmask, nlpeg)
  whupeg = where(umask, nupeg)
  tnmin_fix, whlpeg, whupeg, g
  GTG = TOTAL(G*G)

;
; CHECK IF THE INITIAL POINT IS A LOCAL MINIMUM.
;
  FTEST = ONE + ABS(FNEW)
  IF (GTG LT GTOL2*FTEST*FTEST) THEN GOTO, CLEANUP

;
; SET INITIAL VALUES TO OTHER PARAMETERS
;
  ICYCLE = NM1
  GNORM  = SQRT(GTG)
  DIFNEW = ZERO
  EPSRED = HALF/TEN
  FKEEP  = FNEW

;
; SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY.
;
  LDIAGB = replicate(one, n)


;
; ..................START OF MAIN ITERATIVE LOOP..........
;
; COMPUTE THE NEW SEARCH DIRECTION
;
  catch_msg = 'calling TNMIN_MODLNP'
  tnmin_modlnp, lpk, lgv, lz1, lv, ldiagb, lemat, $
    x, g, lzk, n, niter, maxit, nmodif, nlincg, upd1, yksk, $
    gsk, yrsr, lreset, fcn, whlpeg, whupeg, accrcy, gtpnew, gnorm, xnorm, $
    xnew

ITER_LOOP:
  catch_msg = 'computing step length'
  LOLDG = G
  PNORM = tnmin_enorm(LPK)
  OLDF = FNEW
  OLDGTP = GTPNEW

;
; PREPARE TO COMPUTE THE STEP LENGTH
;
  PE = PNORM + EPSMCH

;
; COMPUTE THE ABSOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH
;
  RELTOL = RTEPS*(XNORM + ONE)/PE
  ABSTOL = - EPSMCH*FTEST/(OLDGTP - EPSMCH)

;
; COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN
; THE LINEAR SEARCH
;
  TNYTOL = EPSMCH*(XNORM + ONE)/PE

  ;; From STPMAX
  SPE = STEPMX/PE
  mmask = (NOT lmask AND NOT umask)
  wh = where(mmask AND (lpk GT 0) AND qulim AND (ulim - x LT spe*lpk), ct)
  if ct GT 0 then begin
      spe = min( (ulim(wh)-x(wh)) / lpk(wh))
  endif
  wh = where(mmask AND (lpk LT 0) AND qllim AND (llim - x GT spe*lpk), ct)
  if ct GT 0 then begin
      spe = min( (llim(wh)-x(wh)) / lpk(wh))
  endif

  ;; From LMQNBC
;
; SET THE INITIAL STEP LENGTH.
;
  ALPHA = TNMIN_STEP1(FNEW,FM,OLDGTP,SPE, epsmch)

;
; PERFORM THE LINEAR SEARCH
;
  catch_msg = 'performing linear search'
  tnmin_linder, n, fcn, small, epsmch, reltol, abstol, tnytol, $
    eta, zero, spe, lpk, oldgtp, x, fnew, alpha, g, nwhy, xnew

  NEWCON = 0
  IF (ABS(ALPHA-SPE) GT 1.D1*EPSMCH) EQ 0 THEN BEGIN
      NEWCON = 1
      NWHY   = 0L

      ;; From MODZ
      mmask = (NOT lmask AND NOT umask)
      wh = where(mmask AND (lpk LT 0) AND qllim $
                 AND (x-llim LE 10*epsmch*(abs(llim)+one)),ct)
      if ct GT 0 then begin
          flast = fnew
          lmask(wh) = 1
          x(wh) = llim(wh)
          whlpeg = where(lmask, nlpeg)
      endif
      wh = where(mmask AND (lpk GT 0) AND qulim $
                 AND (ulim-x LE 10*epsmch*(abs(ulim)+one)),ct)
      if ct GT 0 then begin
          flast = fnew
          umask(wh) = 1
          x(wh) = ulim(wh)
          whupeg = where(umask, nupeg)
      endif
      xnew(ifree) = x

      ;; From LMQNBC
      FLAST = FNEW
  endif

  FOLD = FNEW
  NITER = NITER + 1

;
; IF REQUIRED, PRINT THE DETAILS OF THIS ITERATION
;
  if nprint GT 0 AND iterproc NE '' then begin
      iflag = 0L
      xx = xnew
      xx(ifree) = x
      if (niter-1) MOD nprint EQ 0 then begin
          catch_msg = 'calling '+iterproc
          tnerr = 0
          call_procedure, iterproc, fcn, xx, niter, fnew, $
            FUNCTARGS=fcnargs, parinfo=parinfo, quiet=quiet, $
            dprint=iterd, deriv=g, pfixed=pfixed, maximize=maximize, $
            _EXTRA=iterargs
          iflag = tnerr
          if iflag LT 0 then begin
              errmsg = 'WARNING: premature termination by "'+iterproc+'"'
              nwhy = 4L
              goto, CLEANUP
          endif
      endif
  endif

  catch_msg = 'testing for convergence'
  IF (NWHY LT 0) THEN BEGIN
      NWHY = -2L
      goto, CLEANUP
  ENDIF
  IF (NWHY NE 0 AND NWHY NE 2) THEN BEGIN
      ;; THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT
      NWHY = 3L
      goto, CLEANUP
  endif
  if nwhy GT 1 then begin
      fnew = tnmin_call(fcn, x, g, fullparam_=xnew)
  endif
  wh = where(finite(x) EQ 0, ct)
  if ct GT 0 OR finite(fnew) EQ 0 then begin
      nwhy = -3L
      goto, CLEANUP
  endif

;
; TERMINATE IF MORE THAN MAXFUN EVALUATIONS HAVE BEEN MADE
;
  NWHY = 2L
  if maxfun GT 0 AND tnconfig.nfev GT maxfun then goto, CLEANUP
  if niter GT maxiter then goto, CLEANUP
  NWHY = 0L

;
; SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS
;
  DIFOLD = DIFNEW
  DIFNEW = OLDF - FNEW

;
; IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE
; PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST.
;
  IF (ICYCLE EQ 1) THEN BEGIN
      IF (DIFNEW GT 2.D0*DIFOLD) THEN EPSRED = EPSRED + EPSRED
      IF (DIFNEW LT 5.0D-1*DIFOLD) THEN EPSRED = HALF*EPSRED
  ENDIF
  LGV = G
  tnmin_fix, whlpeg, whupeg, lgv
  GTG = TOTAL(LGV*LGV)
  GNORM = SQRT(GTG)
  FTEST = ONE + ABS(FNEW)
  XNORM = tnmin_enorm(X)

  ;; From CNVTST
  LTEST = (FLAST - FNEW) LE (-5.D-1*GTPNEW)
  wh = where((lmask AND g LT 0) OR (umask AND g GT 0), ct)
  if ct GT 0 then begin
      conv = 0
      if NOT ltest then begin
          mx = max(abs(g(wh)), wh2)
          lmask(wh(wh2)) = 0 & umask(wh(wh2)) = 0
          FLAST = FNEW
          goto, CNVTST_DONE
      endif
  endif
  ;; Gill Murray and Wright tests are listed to the right.
  ;; Modifications due to absolute function value test are done here.

  fconv = abs(DIFNEW) LT FTOL1*FTEST              ;; U1a
  if ftol2 EQ 0 then begin
      pconv = ALPHA*PNORM LT PTOL*(ONE + XNORM)   ;; U2
      gconv = GTG LT GTOL1*FTEST*FTEST            ;; U3
  endif else begin
      ;; Absolute tolerance implies a loser constraint on parameters
      fconv = fconv OR (abs(difnew) LT ftol2)     ;; U1b
      acc2  = (FTOL2/FTEST + EPSMCH)
      pconv = ALPHA*PNORM LT SQRT(acc2)*(ONE + XNORM) ;; U2
      gconv = GTG LT (acc2^twothird)*FTEST*FTEST  ;; U3
  endelse
  IF ((PCONV AND FCONV AND GCONV) $               ;; U1 + U2 + U3
      OR (GTG LT GTOL2*FTEST*FTEST)) THEN BEGIN   ;; U4
      CONV = 1
  ENDIF ELSE BEGIN
      ;; Convergence failed
      CONV = 0
  ENDELSE

;
; FOR DETAILS, SEE GILL, MURRAY, AND WRIGHT (1981, P. 308) AND
; FLETCHER (1981, P. 116).  THE MULTIPLIER TESTS (HERE, TESTING
; THE SIGN OF THE COMPONENTS OF THE GRADIENT) MAY STILL NEED TO
; MODIFIED TO INCORPORATE TOLERANCES FOR ZERO.
;

CNVTST_DONE:
  ;; From LMQNBC

  IF (CONV) THEN GOTO, CLEANUP
  tnmin_fix, whlpeg, whupeg, g

;
; COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE
; IN THE GRADIENTS
;
  IF NEWCON EQ 0 THEN BEGIN
      LYK = G - LOLDG
      LSK = ALPHA*LPK
;
; SET UP PARAMETERS USED IN UPDATING THE PRECONDITIONING STRATEGY.
;
      YKSK = TOTAL(LYK*LSK)
      LRESET = 0

      IF (ICYCLE EQ NM1 OR DIFNEW LT EPSRED*(FKEEP-FNEW)) THEN LRESET = 1
      IF (LRESET EQ 0) THEN BEGIN
          YRSR = TOTAL(LYR*LSR)
          IF (YRSR LE ZERO) THEN LRESET = 1
      ENDIF
      UPD1 = 0
  ENDIF

;
;      COMPUTE THE NEW SEARCH DIRECTION
;
  ;; TNMIN_90:
  catch_msg = 'calling TNMIN_MODLNP'

  tnmin_modlnp, lpk, lgv, lz1, lv, ldiagb, lemat, $
    x, g, lzk, n, niter, maxit, nmodif, nlincg, upd1, yksk, $
    gsk, yrsr, lreset, fcn, whlpeg, whupeg, accrcy, gtpnew, gnorm, xnorm, $
    xnew

  IF (NEWCON) THEN GOTO, ITER_LOOP
;  IF (NOT LRESET) OR ICYCLE EQ 1 AND n_elements(LSR) GT 0 THEN BEGIN   ;; For testing
  IF (LRESET EQ 0) THEN BEGIN
;
; COMPUTE THE ACCUMULATED STEP AND ITS CORRESPONDING
; GRADIENT DIFFERENCE.
;
      LSR = LSR + LSK
      LYR = LYR + LYK
      ICYCLE = ICYCLE + 1
      goto, ITER_LOOP
  ENDIF

;
; RESET
;
  ;; TNMIN_110:
  IRESET = IRESET + 1
;
; INITIALIZE THE SUM OF ALL THE CHANGES IN X.
;
  LSR = LSK
  LYR = LYK
  FKEEP = FNEW
  ICYCLE = 1L
  goto, ITER_LOOP

;
; ...............END OF MAIN ITERATION.......................
;
  CLEANUP:

  nfev = tnconfig.nfev
  tnfcnargs = 0
  catch, /cancel
  case NWHY of
      -3: begin
          ;; INDEFINITE VALUE
          status = -16L
          if errmsg EQ '' then $
            errmsg = ('ERROR: parameter or function value(s) have become '+$
                      'infinite; check model function for over- '+$
                      'and underflow')
          return, !values.d_nan
      end
      -2: begin
          ;; INTERNAL ERROR IN LINE SEARCH
          status = -18L
          if errmsg EQ '' then $
            errmsg = 'ERROR: Fatal error during line search'
          return, !values.d_nan
      end
      -1: begin
          TERMINATE:
          ;; FATAL TERMINATION
          status = 0L
          if errmsg EQ '' then errmsg = 'ERROR: Invalid inputs'
          return, !values.d_nan
      end
      0: begin
          CONVERGED:
          status = 1L
      end
      2: begin
          ;; MAXIMUM NUMBER of FUNC EVALS or ITERATIONS REACHED
          if maxfun GT 0 AND nfev GT maxfun then begin
              status = -17L
              if errmsg EQ '' then $
                errmsg = ('WARNING: no convergence within maximum '+$
                          'number of function calls')
          endif else begin
              status = 5L
              if errmsg EQ '' then $
                errmsg = ('WARNING: no convergence within maximum '+$
                          'number of iterations')
          endelse
          FNEW = OLDF
      end
      3: begin
          status = -18L
          if errmsg EQ '' then errmsg = 'ERROR: Line search failed to converge'
      end
      4: begin
          ;; ABNORMAL TERMINATION BY USER ROUTINE
          status = iflag
      end
  endcase


  ;; Successful return
  F = FNEW
  xnew(ifree) = x
  if keyword_set(tnconfig.qanytied) then tnmin_tie, xnew, tnconfig.ptied

  return, xnew
end