-
Notifications
You must be signed in to change notification settings - Fork 0
/
nlscon.m
4125 lines (4088 loc) · 155 KB
/
nlscon.m
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
function [x,info,wk] = nlscon(m,fcn,x,xscal,fi,fscal,rtol,iopt,par,wk,info)
% ------------------------------------------------------------
%
%* Title
%
% Numerical solution of nonlinear (NL) least squares (S)
% problems with nonlinear constraints (CON), especially
% designed for numerically sensitive problems.
%
%* Written by U. Nowak, L. Weimann
%* Purpose Solution of highly nonlinear, optionally
% constrained least squares problems
%* Method Damped affine invariant Gauss-Newton method
% (see references below)
%* Category K1b2b. - Nonlinear least squares approxi-
% mation with nonlinear constraints
%* Keywords Nonlinear least squares problems,
% Gauss-Newton methods
%* Version 2.3.2
%* Revision December 1993
%* Latest Change August 2001
%* Library CodeLib
%* Code Matlab 6.0
%* Environment Standard Matlab 6.0 environment on PC's,
% workstations and hosts.
%* Copyright (c) Konrad-Zuse-Zentrum fuer
% Informationstechnik Berlin (ZIB)
% Takustrasse 7, D-14195 Berlin-Dahlem
% phone : + 49/30/84185-0
% fax : + 49/30/84185-125
%* Contact Lutz Weimann
% ZIB, Division Scientific Computing,
% Department Numerical Analysis and Modelling
% phone : + 49/30/84185-185
% fax : + 49/30/84185-107
% e-mail: [email protected]
%
%* References:
%
% /1/ P. Deuflhard:
% Newton Methods for Nonlinear Problems. -
% Affine Invariance and Adaptive Algorithms.
% Series Computational Mathematics 35, Springer (2004)
%
% /2/ U. Nowak, L. Weimann:
% A Family of Newton Codes for Systems of Highly Nonlinear
% Equations - Algorithm, Implementation, Application.
% ZIB, Technical Report TR 90-10 (December 1990)
%
% ---------------------------------------------------------------
%
%* Licence
% You may use or modify this code for your own non commercial
% purposes for an unlimited time.
% In any case you should not deliver this code without a special
% permission of ZIB.
% In case you intend to use the code commercially, we oblige you
% to sign an according licence agreement with ZIB.
%
%* Warranty
% This code has been tested up to a certain level. Defects and
% weaknesses, which may be included in the code, do not establish
% any warranties by ZIB. ZIB does not take over any liabilities
% which may follow from acquisition or application of this code.
%
%* Software status
% This code is under care of ZIB and belongs to ZIB software class 1.
%
% ------------------------------------------------------------
%
%* Summary:
% ========
% Damped Gauss-Newton-algorithm with rank strategy for highly
% nonlinear least squares approximation problems (optionally
% constrained, over- and underdetermined problems) -
% due to Ref.(1).
%
% (The iteration is done by function NCINT currently. NLSCON
% itself does some house keeping and builds up workspace.)
%
% The problem solved by this program looks as follows:
%
% Denoting below the n-dimensional real space with IR(n),
% the number of parameters to be estimated with N,
% the number of measurements (data to be fitted) with MFIT, and
% the number of equality constraints with MCON,
% M := MCON + MFIT ,
% let F : IR(N) --> IR(MFIT) , G : IR(N) --> IR(MCON)
% be nonlinear functions,
% FI in IR(MFIT) the vector of measurement data and
% FSCAL in IR(MFIT) the vector of measurement weights.
%
% For M >= N, find a parameter vector X in IR(N), which
% minimizes Sum (j=1,...,MFIT) ((Fj(X)-FIj)/FSCALj)**2 and
% satisfies G(X) = 0.
% For M < N, find the parameter vector X in IR(N) with the
% smallest possible euclidian norm,
% which satisfies F(X) = 0 and G(X) = 0 .
%
% Jacobian approximation by numerical differences or user
% supplied function FCN, called with 'jacobian'-flag.
%
% The numerical solution of the arising linear least squares
% problem is done by means of the functions DECCON and SOLCON
% (QR decomposition with subcondition estimation, rank decision
% and computation of the rank-deficient pseudoinverse) .
% For special purposes these routines may be substituted.
%
% A statistical a posteriori analysis of the parameter estimate
% is optionally available.
%
% This is a driver routine for the core solver NCINT.
%
% ------------------------------------------------------------
%
%* Parameters list description (* marks inout parameters)
% ======================================================
%
%* External functions (to be supplied by the user)
% =================================================
%
% [FOUT,FAIL] =FCN(X,FLAG,PAR) String Name of problem- and Jacobian function
% X() Float Vector of unknowns (input)
% FLAG String Operation flag - the following values must be
% supported:
% '' (empty string) : FOUT must return the problem-
% function value f(X) ;
% 'jacobian' : FOUT must return the associated
% Jacobian Jac(x)
% PAR AnyType A (required!) user parameter
%
% FOUT() Float Vector of returned function values or Jacobian matrix
% FAIL Int evaluation-failure indicator. (output)
% On output: Indicates failure of FCN eval-
% uation, if having a value <= 2.
% If <0: NLEQ1 will be terminated with
% error code = 82, and FAIL stored
% to wk.ifail.
% If =1: A new trial Newton iterate will
% computed, with the damping factor
% reduced to its half.
% If =2: A new trial Newton iterate will
% computed, with the damping factor
% reduced by a reduct. factor, which
% must be output through F(1) by FCN,
% and it is value must be >0 and < 1.
% Note, that if FAIL = 1 or 2, additional
% conditions concerning the damping factor,
% e.g. the minimum damping factor or the
% bounded damping strategy may also influ-
% ence the value of the reduced damping
% factor.
%
%
%* Input parameters of NLSCON
% ==========================
%
% M Int Sum of number of measurement data and
% equality constraints
% X(N) Dble Initial estimate of parameters
% XSCAL(N) Dble User scaling (lower threshold) of the
% iteration vector X(N)
% FI(MFIT) Dble Data obtained by measurements
% FSCAL(MFIT) Dble User weighting vector of measurements
% RTOL Dble Required relative precision of
% solution components -
% RTOL >= EPMACH*TEN*N
% IOPT(50) Int Structure of run-time options. Set the
% fields to zero or omit them to get
% default values (details see below)
% PAR AnyTyp This parameter will be passed on every
% call of the user function FCN as the
% third input argument.
% WK Struct A structured workspace. Normally an empty
% array should be passed to this argument.
% However, on successor steps in onestep mode,
% the output parameter WK from the previous step
% call must be passed to this parameter.
%
%* Output parameters of NLSCON
% ===========================
%
% X(N) Dble Solution parameters ( or final values,
% respectively )
% INFO Struct A structure, holding additional info
% obtained from the iteration run. The
% fields are as follows:
% XSCAL(N) Dble After return with IERR >= 0, it contains
% the latest internal scaling vector used
% After return with IERR == -1 in onestep-
% mode it contains a possibly adapted
% (as described below) user scaling vector:
% if XSCAL(I) < SMALL) XSCAL(I) = SMALL ,
% if XSCAL(I) > GREAT) XSCAL(I) = GREAT .
% For SMALL and GREAT, see section machine
% constants below and regard note 1.
% RTOL Dble Finally achieved (relative) accuracy
% The estimated absolute error of component i
% of x_out is approximately given by
% abs_err(i) = RTOL * XSCAL_out(i) ,
% where (approximately)
% XSCAL_out(i) =
% max(abs(X_out(i)),XSCAL_in(i)).
% Note that RTOL_out may be greater than
% RTOL_in, but NLSCON claims 'solution found'
% - see IOPT(36).
% IERR Int Return value parameter
% =-1 sucessfull completion of one iteration
% step, subsequent iterations are needed
% to get a solution. (stepwise mode only)
% = 0 successfull completion of iteration
% > 0 see list of error messages below
% XITER(N,:) Dble An array holding all iterates of the
% Newton iteration run, stored as column
% vectors.
% NATLEVEL(:) Dble The sequence of natural levels of the
% newton corrections over the iteration steps
% SIMLEVEL(:) Dble The sequence of natural levels of the simplified
% newton corrections over the iteration steps
% STDLEVEL(:) Dble The sequence of standard levels over the
% iteration steps
% PRECISION(:) Dble The sequence of acheived precisions over
% the iteration steps.
% DAMPINGFC(:) Dble The sequence of accepted damping factors
% over the iteration steps.
% NITER Int Number of Newton-iterations done
% NCORR Int Number of corrector steps done
% NREJR1 Int Number of rejected Newton iteration steps
% done with a rank-1 approximated Jacobian
% NJAC Int Number of Jacobian generations or JAC-calls done
% NFCN Int Number of FCN-evaluations done
% NFCNJ Int Number of FCN-evaluations for Jacobian
% approximation done
%
% Note 1.
% The machine dependent values SMALL, GREAT and EPMACH are
% gained from calls of the machine constants function D1MACH.
% As delivered, this function is adapted to use constants
% suitable for all machines with IEEE arithmetic. If you use
% another type of machine, you have to change the DATA state-
% ments for IEEE arithmetic in D1MACH into comments and to
% uncomment the set of DATA statements suitable for your machine.
%
%* Options IOPT.name:
% ==================
%
% Name Default Meaning
%
% mode 0 =0 Standard mode initial call:
% Return when the required accuracy for the
% iteration vector is reached. User defined
% parameters are evaluated and checked.
% Standard mode successive call:
% If NLSCON was called previously with
% MODE=1, it performs all remaining
% iteration steps.
% =1 Stepwise mode:
% Return after one Gauss-Newton
% iteration step.
% jacgen 0 Method of Jacobian generation
% =0 Standard method is JACGEN=2
% =1 User supplied function FCN will be
% called to generate Jacobian matrix
% =2 Jacobian approximation by numerical
% differentation (see function NCJAC)
% =3 Jacobian approximation by numerical
% differentation with feedback control
% (see function NCJCF)
% iscal 0 Determines how to scale the iterate-vector:
% =0 The user supplied scaling vector XSCAL is
% used as a (componentwise) lower threshold
% of the current scaling vector
% =1 The vector XSCAL is always used as the
% current scaling vector
% mprerr 0 Print error messages
% =0 No output
% =1 Error messages
% =2 Warnings additionally
% =3 Informal messages additionally
% luerr 1 Logical unit number for error messages
% mprmon 0 Print iteration Monitor
% =0 No output
% =1 Standard output
% =2 Summary iteration monitor additionally
% =3 Detailed iteration monitor additionally
% =4,5,6 Outputs with increasing level addi-
% tional increasing information for code
% testing purposes. Level 6 produces
% in general extremely large output!
% lumon 1 Logical unit number for iteration monitor
% mprsol 0 Print solutions
% =0 No output
% =1 Initial values and solution values
% =2 Intermediate iterates additionally
% lusol 1 Logical unit number for solutions
% mprtim 0 Output level for the time monitor
% = 0 : no time measurement and no output
% = 1 : time measurement will be done and
% summary output will be written -
% regard note 4a.
% lutim 1 Logical output unit for time monitor
% qstat 0 Statistical Analysis of the final
% least squares estimate:
% = 0 : Analysis will not be done
% = 1 : Analysis will be done,
% and certain results are stored
% to the RWK array (for details, see
% RWK description below)
% mprsta 0 Printing of statistical Analysis for
% the final least squares estimate:
% = 0 : Printing will not be done
% = 1 : Printing will be done (and this
% implies QSTAT to be set to 1)
% nonlin 3 Problem type specification
% =1 Linear problem
% Warning: If specified, no check will be
% done, if the problem is really linear, and
% NLSCON terminates unconditionally after
% one Gauss-Newton-iteration step.
% =2 Mildly nonlinear problem
% =3 Highly nonlinear problem
% =4 Extremely nonlinear problem
% qrank1 0 =0 (0) Rank-1 updates by Broyden-
% approximation are inhibited.
% =1 (1) Rank-1 updates by Broyden-
% approximation are allowed.
% qnscal 0 Inhibit automatic row scaling:
% =0 (0) Automatic row scaling of
% the linear system is activ:
% Rows i=1,...,N will be divided by
% max j=1,...,N (abs(a(i,j)))
% =1 (1) No row scaling of the linear
% system. Recommended only for well row-
% scaled nonlinear least squares problems.
% iterm 0 Determines the iteration termination cri-
% terium to be chosen:
% =0 Iteration is terminated, if one of the
% stopping criteria used for ITERM=1 and
% ITERM=2 is satisfied (see below).
% =1 Iteration is terminated, if, from a
% statistical point of view, a resonable
% precision is achieved, i.e.
% simplified GN-correction < RTOL, and an
% estimate of the accuracy is available.
% Recommended to be used for incompatible
% problems.
% =2 Iteration is terminated, if
% GN-correction < RTOL . Using this option
% may force 'to much' precision from
% the statistical point of view.
% ibdamp Bounded damping strategy switch:
% =0 means currently always IBDAMP = off
% (but may depend on the settings of other
% options in future versions)
% =1 means always IBDAMP = on
% =2 means always IBDAMP = off
% Note 3:
% If NLSCON terminates with IERR=2 (maximum iterations)
% or IERR=3 (small damping factor), you may try to continue
% the iteration by increasing NITMAX or decreasing FCMIN
% (see WK) and setting WK.QSUCC to 1.
%
%* Optional input/output in WK:
% ============================
%
% Name Meaning
%
% niter IN/OUT Number of Gauss-Newton-iterations
% ncorr IN/OUT Number of corrector steps
% nfcn IN/OUT Number of FCN-evaluations
% njac IN/OUT Number of Jacobian generations or
% JAC-calls
% nfcnj IN/OUT Number of FCN-evaluations for Jacobian
% approximation
% nrejr1 IN/OUT Number of rejected Gauss-Newton iteration
% steps done with a rank-1 computed Jacobian
% idcode IN/OUT Output: The 8 decimal digits program identi-
% fication number ppppvvvv, consisting of the
% program code pppp and the version code vvvv.
% Input: If containing a negative number,
% it will only be overwritten by the identi-
% fication number, immediately followed by
% a return to the calling program.
% ifail OUT Set in case of failure of NCFACT (IERR=80),
% N2SOLV (IERR=81), FCN (IERR=82) or JAC(IERR=83)
% to the nonzero IFAIL value returned by the
% routine indicating the failure .
% nitmax IN Maximum number of permitted iteration
% steps (default: 50)
% irank IN Initial rank
% =0 : full rank N
% =k with 0 < k < N : deficient rank assumed
% for the Jacobian in the starting point
% new IN/OUT Count of consecutive rank-1 updates
% ifccnt INTERN Count of consecutive special steps before
% convergence stop of iteration
% conv OUT The maximum norm of the latest ordinary
% respective simplified (scaled) Gauss-Newton
% correction.
% sumx OUT Natural level (not Normx of printouts)
% of the current iterate, i.e. Sum(DX(i)**2),
% where DX = scaled Gauss-Newton correction.
% dlevf OUT Standard level (not Normf of printouts)
% of the current iterate, i.e. Norm2(F(X)),
% where F = nonlinear model function.
% fcbnd IN Bounded damping strategy restriction factor
% (Default is 10)
% fcstrt IN Damping factor for first Gauss-Newton iteration
% - overrides option NONLIN, if set (see note 5)
% fcmin IN Minimal allowed damping factor (see note 5)
% sigma IN Broyden-approximation decision parameter
% Required choice: SIGMA >= 1. Increasing this
% parameter make it less probable that the algo-
% rith performs Broyden steps.
% Rank1 updates are inhibited, if
% SIGMA > 1/FCMIN is set. (see note 5)
% cond IN Maximum permitted subcondition for rank-
% decision by linear solver
% (Default: 1/epmach, epmach: relative
% machine precision)
% ajdel IN Jacobian approximation without feedback:
% Relative pertubation for components
% (Default: sqrt(epmach*10), epmach: relative
% machine precision)
% ajmin IN Jacobian approximation without feedback:
% Threshold value (Default: 0.0)
% The absolute pertubation for component k is
% computed by
% DELX := AJDEL*max(abs(Xk),AJMIN)
% etadif IN Jacobian approximation with feedback:
% Target value for relative pertubation ETA of X
% (Default: 1.0d-6)
% etaini IN Jacobian approximation with feedback:
% Initial value for denominator differences
% (Default: 1.0d-6)
% prec OUT An estimate for the achieved relative accuracy.
% This number is only available, if IERR=0 or
% IERR=1 and an estimate for the incompatibility
% factor kappa (SKAP=RWK(32)) can be made. If no
% meaningful information is at hand, a -1.0 is
% stored.
% skap OUT An estimate of the incompatibility factor of
% the given nonlinear least squares problem.
% This number is only available, if IERR=0 or
% IERR=1 and a certain further condition holds.
% If no meaningful information is at hand, a -1.0
% is stored.
% A small value of SKAP indicates that the given
% measurement data can be well fitted by the
% model, whereas a value SKAP >= 0.5 may give a
% hint to an insufficient modelling of the
% experiment from which the measurements
% originate.
% sigma2 OUT Holds the estimated variance of the residual
% on final exit of NLSCON, if IOPT(21)=1 is set.
% xl(1:n) OUT Holds the left bounds of the final parameters
% confidence intervals;
% on final exit of NLSCON, if IOPT(21)=1 is set.
% xr(1:n) OUT Holds the right bounds of the final parameters
% confidence intervals;
% on final exit of NLSCON, if IOPT(21)=1 is set.
% vcv(1:n,1:n)
% OUT Holds the columnwise stored correlation-matrix
% on final exit of NLSCON, if IOPT(21)=1 is set.
%
% Note 5:
% The default values of the internal parameters may be obtained
% from the monitor output with at least IOPT field MPRMON set to 2.
%
%* Error messages:
% ===============
%
% 1 Termination at stationary point (rank deficient Jacobian
% and termination criterion fulfilled)
% 2 Termination after NITMAX iterations ( as indicated by
% input parameter NITMAX=IWK(31) )
% 3 Termination, since damping factor became to small and
% Jacobian rank was already reduced to zero
% 20 Bad or inconsistent input to the dimensional parameter M
% 21 Nonpositive value for RTOL supplied
% 22 Negative scaling value via vector XSCAL supplied
% 30 One or more fields specified in IOPT are invalid
% (for more information, see error-printout)
% 80 Error signalled by linear solver routine NCFACT,
% for more detailed information see IFAIL-value
% stored to IWK(23)
% (not used with standard routine NCFACT)
% 81 Error signalled by linear solver routine NCFIT,
% for more detailed information see IFAIL-value
% stored to IWK(23)
% (not used with standard routine NCFIT)
% 82 Error signalled by user routine FCN when called
% with ''-flag (blank-flag) (Nonzero value
% returned via IFAIL-flag; stored to WK.ifail )
% 83 Error signalled by user routine FCN when called
% with 'jacobian'-flag (Nonzero value
% returned via IFAIL-flag; stored to WK.ifail )
% 180,182,183
% see error codes 80,82,83, but the failure of NCFACT or
% FCN occured when preparing the call of the statistics
% function STACON for the final iterate of NLSCON.
%
% Note 6 : in case of failure:
% - use non-standard options
% - use another initial guess
% - or reformulate model
% - or turn to general optimization routine
%
%* Machine dependent constants used:
% =================================
%
% DOUBLE PRECISION EPMACH in NLSCON, NCPCHK, NCINT
% DOUBLE PRECISION GREAT in NLSCON, NCPCHK
% DOUBLE PRECISION SMALL in NLSCON, NCPCHK, NCINT, NCSCAL
%
%* Functions called: NCPCHK, NCINT, D1MACH
%
% ------------------------------------------------------------
%* End Prologue
%
%* Summary of changes of the underlying Fortran code:
% ==================================================
%
% 2.2.1 91, June 3 Time monitor included
% 2.2.2 91, June 3 Bounded damping strategy implemented
% 2.2.3 91, July 26 AJDEL, AJMIN as RWK-options for JACGEN == 2,
% ETADIF, ETAINI as RWK-opts. for JACGEN == 3
% FCN-count changed for anal. Jacobian
% 91, Sept. DECCON with new fail exit, for the case that
% the square root of a negative number will
% appear during pseudo inverse computation.
% (Occured, although theoretical impossible!)
% 2.2.6 91, Sept. 17 Damping factor reduction by FCN-fail imple-
% mented
% Sept. 24 Meaning of option ITERM modified.
% 2.3 91, Dec. 20 New Release for CodeLib
% 92, June 2 Level of accepted simplified correction
% stored to RWK(IRWKI+4)
% 2.3.1 92, Oct. 13 Corrected errors in subr. NCJAC, NCJCF:
% dimensions of FX(M), FU(M),
% description of FCN in these functions
% 2.3.2 93, Nov. 24 Optional call of statistical analysis
% routine STACON to get statistics about
% the quality of the parameter estimate
% which has been computed by NLSCON.
% 00, July 12 RTOL output-value bug fixed
% 05, Sept. 12 FCREDU setting bug near 3.6.1 fixed
% 2.3.3 06, Nov. 23 Spurious number output fixed, and made
% several changes recommended by M-Lint
% tool of MATLAB version 7.2.0.
%
% ------------------------------------------------------------
%
% Further WK positions (only internally used)
%
% Name Meaning
%
% FCA Previous damping factor
% SUMXS natural level of accepted simplified correction
%
% Internal arrays stored in WK (see routine NCINT for descriptions)
%
% Array Type Remarks
%
% AA(M,N) Perm
% DX(N) Perm
% DXQ(N) Perm
% XA(N) Perm
% FMODEL(M) Perm
% F(M) Perm
% FA(M) Perm
% FW(M) Perm
% ETA(N) Perm Only used for JACGEN=3
% A(M,N) Temp
% QA(N,N) Temp
% XW(N) Temp
% DXQA(N) Temp
% QU(M) Temp
% RQ(M) Temp
% RQKEEP(M) Temp
% DELXQ(N) Temp
%
%
persistent xiter sumxall dlevfall sumxqall tolall fcall ;
one = 1.0 ;
ten = 10.0 ;
zero = 0.0 ;
iver = 22112323 ;
%
% Version: 2.3.3 Latest change:
% -----------------------------------------
%
chgdat = 'November 23, 2006 ' ;
prodct = 'NLSCON ' ;
%* Begin
n = length(x) ;
mfit = length(fi) ;
epmach = d1mach(3) ;
great = sqrt(d1mach(2)/ten) ;
small = d1mach(6) ;
tolmin = epmach*ten*n ;
ierr = 0 ;
if isfield(wk,'iver')
qvchk = wk.iver < 0 ;
wk.iver = iver ;
if qvchk
huch=qvchk
return ;
end ;
else
wk.iver = iver ;
end
% Print error messages?
[mprerr,iopt] = getopt(iopt,'mprerr',0) ;
[luerr, iopt] = getopt(iopt,'luerr', 1) ;
if luerr == 0
luerr = 1 ;
iopt.luerr = luerr ;
end
% print iteration monitor?
[mprmon,iopt] = getopt(iopt,'mprmon',0) ;
[lumon ,iopt] = getopt(iopt,'lumon', 1) ;
if lumon <= 0 || lumon > 99
lumon = 1 ;
iopt.lumon = lumon ;
end
% print intermediate solutions?
[mprsol,iopt] = getopt(iopt,'mprsol',0) ;
[lusol, iopt] = getopt(iopt,'lusol', 1) ;
if lusol == 0
lusol = 1 ;
iopt.lusol=lusol ;
end
% print time summary statistics?
[mprtim,iopt] = getopt(iopt,'mprtim',0) ;
[lutim, iopt] = getopt(iopt,'lutim', 1) ;
if lutim == 0
lutim = 1 ;
iopt.lutim=lutim ;
end
[qsucc,wk] = getopt(wk,'qsucc',0) ;
qinimo = mprmon >= 1 && ~ qsucc ;
% Print NLSCON heading lines
if qinimo
fprintf(lumon,'%s\n\n%s\n\n',' N L S C O N ***** V e r s i o n 2 . 3 . 2 ***', ...
' Gauss-Newton-Method for the solution of nonlinear least squares problems') ;
end
% Check input parameters and options
[ierr,iopt,rtol,xscal,fscal] = ncpchk(n,m,mfit,x,xscal,fi,fscal,rtol,iopt,wk) ;
% Exit, if any parameter error was detected till here
if ierr ~= 0
huch1=ierr
return ;
end
if ~ qsucc
xiter = [] ;
sumxall = [] ;
dlevfall = [] ;
sumxqall = [] ;
tolall = [] ;
fcall = [] ;
end
mcon = m-mfit ;
m1 = m ;
m2 = m ;
maxmn = max(m,n) ;
[jacgen,iopt] = getopt(iopt,'jacgen',0) ;
if jacgen == 0
jacgen=2 ;
end
iopt.jacgen = jacgen ;
% WorkSpace: WK
wk = iniopt(wk,'xl', zeros(n,1)) ;
wk = iniopt(wk,'xr', zeros(n,1)) ;
% wk = iniopt(wk,'vcv', zeros(n,n)) ;
% wk = iniopt(wk,'aa', zeros(m2,n)) ;
% wk = iniopt(wk,'a', zeros(m1,n)) ;
% wk = iniopt(wk,'qa', zeros(n,n)) ;
wk = iniopt(wk,'vcv', 0.0) ;
wk = iniopt(wk,'aa', 0.0) ;
wk = iniopt(wk,'a', 0.0) ;
wk = iniopt(wk,'qa', 0.0) ;
wk = iniopt(wk,'dx', zeros(n,1)) ;
wk = iniopt(wk,'dxq', zeros(n,1)) ;
wk = iniopt(wk,'xa', zeros(n,1)) ;
wk = iniopt(wk,'xwa', zeros(n,1)) ;
wk = iniopt(wk,'fmodel',zeros(m,1)) ;
wk = iniopt(wk,'f', zeros(m,1)) ;
wk = iniopt(wk,'fa', zeros(m,1)) ;
wk = iniopt(wk,'eta', zeros(n,1)) ;
wk = iniopt(wk,'xw', zeros(n,1)) ;
wk = iniopt(wk,'fw', zeros(m,1)) ;
wk = iniopt(wk,'dxqa', zeros(n,1)) ;
wk = iniopt(wk,'qu', zeros(m,1)) ;
wk = iniopt(wk,'rq', zeros(m,1)) ;
wk = iniopt(wk,'rqkeep',zeros(m,1)) ;
wk = iniopt(wk,'delxq', zeros(n,1)) ;
wk = iniopt(wk,'fcmin',0.0) ;
wk = iniopt(wk,'sigma',0.0) ;
wk = iniopt(wk,'fca',0.0) ;
wk = iniopt(wk,'conv',0.0) ;
wk = iniopt(wk,'sumx',0.0) ;
wk = iniopt(wk,'sumxs',0.0) ;
wk = iniopt(wk,'dlevf',0.0) ;
wk = iniopt(wk,'niter',0) ;
wk = iniopt(wk,'ncorr',0) ;
wk = iniopt(wk,'nfcn',0) ;
wk = iniopt(wk,'njac',0) ;
wk = iniopt(wk,'nfcnj',0) ;
wk = iniopt(wk,'nrejr1',0) ;
wk = iniopt(wk,'new',0) ;
wk = iniopt(wk,'iconv',0) ;
wk = iniopt(wk,'ifccnt',0) ;
%
iopt = iniopt(iopt,'norowscal',0) ;
if qinimo
fprintf(lumon,'\n %s : %4i\n %s : %4i\n %s : %4i\n\n %s : %10.2e\n', ...
'Number of parameters to be estimated (N)',n, ...
'Number of data to fitted, e.g. observations (MFIT)',mfit, ...
'Number of equality constraints (MCON) : ',mcon, ...
'Prescribed relative precision (RTOL) : ',rtol) ;
if jacgen == 1
jacg = 'a user function' ;
elseif jacgen == 2
jacg = 'numerical differentiation (without feedback strategy)' ;
elseif jacgen == 3
jacg = 'numerical differentiation (feedback strategy included)' ;
end
fprintf(lumon,'\n The Jacobian is supplied by %s\n',jacg) ;
if iopt.norowscal == 1
rsmode = 'inhibited' ;
else
rsmode = 'allowed' ;
end
fprintf(lumon,' Automatic row scaling of the jacobian is %s\n',rsmode) ;
end
[qrank1,iopt] = getopt(iopt,'qrank1',0) ;
[nonlin,iopt]=getopt(iopt,'nonlin',3) ;
iopt = iniopt(iopt,'boundeddamp',0) ;
if iopt.boundeddamp == 0
qbdamp = nonlin == 4 ;
elseif iopt.boundeddamp == 1
qbdamp = 1 ;
elseif iopt.boundeddamp == 2
qbdamp = 0 ;
end
wk = iniopt(wk,'fcbnd',0.0) ;
if qbdamp
if wk.fcbnd < 1.0
wk.fcbnd = 10.0 ;
end
end
if qrank1 && m > n
if mprerr >= 2
fprintf(luerr,'\n\n Warning: Broyden-steps set to inhibited for the overdetermined system') ;
end ;
qrank1 = 0 ;
iopt.qrank1 = 0 ;
end
if qinimo
if qrank1
rk1mode = 'allowed' ;
else
rk1mode = 'inhibited' ;
end
fprintf(lumon,'\n Rank-1 updates are %s\n',rk1mode) ;
if nonlin == 1
nonlinmode = 'linear' ;
elseif nonlin == 2
nonlinmode = 'mildly nonlinear' ;
elseif nonlin == 3
nonlinmode = 'highly nonlinear' ;
elseif nonlin == 4
nonlinmode = 'extremely nonlinear' ;
end
fprintf(lumon,' Problem is specified as being %s\n',nonlinmode) ;
if qbdamp
fprintf(lumon, ...
' Bounded damping strategy is active\n Bounding factor is %10.3e\n', wk.fcbnd) ;
else
fprintf(lumon, ' Bounded damping strategy is off\n') ;
end
end
% Maximum permitted number of iteration steps
[nitmax,wk]=getopt(wk,'nitmax',50) ;
if nitmax <= 0
nitmax=50 ;
end
wk.nitmax=nitmax ;
if qinimo
fprintf(lumon,' Maximum permitted number of iteration steps : %6i\n',nitmax) ;
end
% Initial damping factor for highly nonlinear problems
wk = iniopt(wk,'fcstart',0.0) ;
qfcstr = wk.fcstart > 0.0 ;
if ~ qfcstr
wk.fcstart=1.0e-2 ;
if nonlin == 4
wk.fcstart=1.0d-4 ;
end
end
% Minimal permitted damping factor
wk = iniopt(wk,'fcmin',0.0) ;
if wk.fcmin <= 0.0
wk.fcmin=1.0e-4 ;
if nonlin == 4
wk.fcmin=1.0e-8 ;
end
end
fcmin=wk.fcmin ;
% Broyden-update decision parameter SIGMA
wk = iniopt(wk,'sigma',0.0) ;
if wk.sigma < 1.0
wk.sigma=2.0 ;
end
if ~ qrank1
wk.sigma=10.0/fcmin ;
end
% Starting value of damping factor (FCMIN <= FC <= 1.0)
if nonlin <= 2 && ~ qfcstr
% for linear or mildly nonlinear problems
fc = 1.0 ;
else
% for highly or extremely nonlinear problems
fc = wk.fcstart ;
end
wk.fcstart = fc ;
% Initial rank
wk = iniopt(wk,'irank',0) ;
irank = wk.irank ;
minmn = min(m,n) ;
if irank <= 0 || irank > minmn
wk.irank = minmn ;
end
% Maximum permitted sub condition number of matrix A
[cond,wk] = getopt(wk,'cond',0.0) ;
if cond < one
cond = one/epmach ;
end
wk.cond = cond ;
if mprmon >= 2 && ~ qsucc
fprintf(lumon,'\n\n%s\n\n%s%9.2e\n%s%9.2e\n%s%9.2e\n%s%6i\n%s%9.2e\n', ...
' Internal parameters:', ...
' Starting value for damping factor FCSTART = ',wk.fcstart, ...
' Minimum allowed damping factor FCMIN = ',fcmin, ...
' Rank-1 updates decision parameter SIGMA = ',wk.sigma, ...
' Initial Jacobian pseudo-rank IRANK =',wk.irank, ...
' Maximum permitted subcondition COND = ',cond) ;
end
if ~ qsucc
for i=1:mfit
if fscal(i) >= small && fscal(i) <= great
wk.fw(mcon+i) = one/fscal(i) ;
else
wk.fw(mcon+i) = one ;
if fscal(i) ~= zero && mprerr >= 2
fprintf(luerr,'\n Warning: Bad scaling value fscal(%5i) = %10.2e replaced by 1.0\n', ...
i, fscal(i) ) ;
end
end
end
end
%
% Initialize and start time measurements monitor
%
dummy = 1 ;
if (~ qsucc) && mprtim ~= 0
monini(' NLSCON',lutim) ;
mondef(0,'NLSCON') ;
mondef(1,'FCN') ;
mondef(2,'Jacobi') ;
mondef(3,'Lin-Fact') ;
mondef(4,'Lin-Sol') ;
mondef(5,'Output') ;
monsrt(dummy) ;
end
%
ierr = -1 ;
% If IERR is unmodified on exit, successive steps are required
% to complete the Gauss-Newton iteration
[x,xscal,fi,rtol,wk.irank,ierr,wk.xl,wk.xr,wk.vcv, ...
wk.aa,wk.a,wk.qa,wk.dx,wk.dxq,wk.xa,wk.fmodel, ...
wk.f,wk.fa,wk.eta,wk.fw,wk.xw,wk.dxqa,wk.qu, ...
wk.rq,wk.rqkeep,wk.delxq,wk.fcstart,wk.fcmin,wk.sigma,wk.fca,cond, ...
wk.conv,wk.sumx,wk.sumxs,wk.dlevf, ...
xiter,sumxall,dlevfall,sumxqall,tolall,fcall, ...
wk.niter,wk.ncorr,wk.nfcn,wk.njac,wk.nfcnj,wk.nrejr1,wk.new,wk.ifccnt, ...
wk.qsucc] = ...
ncint(n,m,mfit,fcn,par,x,xscal,fi,rtol,nitmax,nonlin, ...
wk.irank,iopt,ierr,wk,m1,m2,maxmn, ...
wk.aa,wk.a,wk.qa,wk.dx,wk.dxq,wk.xa,wk.fmodel, ...
wk.f,wk.fa,wk.eta,wk.fw,wk.xw,wk.dxqa,wk.qu, ...
wk.rq,wk.rqkeep,wk.delxq,wk.fcstart,wk.fcmin,wk.sigma,wk.fca,cond, ...
wk.conv,wk.sumx,wk.sumxs,wk.dlevf, ...
xiter,sumxall,dlevfall,sumxqall,tolall,fcall, ...
tolmin,mprerr,mprmon,mprsol,luerr,lumon,lusol,wk.niter, ...
wk.ncorr,wk.nfcn,wk.njac,wk.nfcnj,wk.nrejr1,wk.new,wk.ifccnt,qbdamp) ;
%
if mprtim ~= 0 && ierr ~= -1 && ierr ~= 10
monend(dummy) ;
end
% Print statistics
if mprmon >= 1 && ierr ~= -1 && ierr ~= 10
fprintf(lumon,'\n%s%s%s\n%s%7i%s\n%s%7i%s\n%s%7i%s\n%s%7i%s\n%s%7i%s\n%s%7i%s\n%s\n\n', ...
' ****** Statistics * ',prodct, ' *******', ...
' *** Gauss-Newton iter.: ',wk.niter, ' ***', ...
' *** Corrector steps : ',wk.ncorr, ' ***', ...
' *** Rejected rk-1 st. : ',wk.nrejr1, ' ***', ...
' *** Jacobian eval. : ',wk.njac, ' ***', ...
' *** Function eval. : ',wk.nfcn, ' ***', ...
' *** ... for Jacobian : ',wk.nfcnj, ' ***', ...
' *************************************') ;
end
info.xscal = xscal ;
if ierr == -1
info.rtol = tolall(wk.niter) ;
else
info.rtol = rtol ;
end
info.ierr = ierr ;
info.xiter = xiter ;
info.natlevel = sumxall ;
info.simlevel = sumxqall ;
info.stdlevel = dlevfall ;
info.precision = tolall ;
info.dampingfc = fcall ;
info.niter = wk.niter ;
info.ncorr = wk.ncorr ;
info.nrejr1 = wk.nrejr1 ;
info.njac = wk.njac ;
info.nfcn = wk.nfcn ;
info.nfcnj = wk.nfcnj ;
% End of function NLSCON
function [ierr,iopt,rtol,xscal,fscal] = ncpchk(n,m,mfit,x,xscal,fi,fscal,rtol,iopt,wk)
% ------------------------------------------------------------
%
%* Summary :
%
% N C P C H K : Checking of input parameters and options
% for NLSCON.
%
%* Parameters:
% ===========
%
% See parameter descriptions in driver routine and NCINT.
%
%* functions called: D1MACH
%
%* Machine dependent constants used:
% =================================
%
% EPMACH = relative machine precision
% GREAT = squareroot of maxreal divided by 10
% SMALL = squareroot of "smallest positive machine number
% divided by relative machine precision"
%
% ------------------------------------------------------------
numopt=50 ;
%
% DATA IOPTL = [0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,0,0,1, ...
% 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, ...
% 0,0,0,0,0,0,0,0,0,0, ...
% -9999999,-9999999,-9999999,-9999999,-9999999]
% DATA IOPTU = [1,1,3,0,0,0,0,0,1,0,3,99,6,99,3,99,0,0,1,99, ...
% 1,1,0,0,0,0,0,0,0,0,4,1,0,0,1, ...
% 2,0,2,0,0,0,0,0,0,0, ...
% 9999999,9999999,9999999,9999999,9999999]
%
epmach = d1mach(3) ;
great = sqrt(d1mach(2)/10.0) ;
small = d1mach(6) ;
ierr = 0 ;
% print error messages?
mprerr = iopt.mprerr ;
luerr = iopt.luerr ;
if luerr <= 0 || luerr > 99
luerr = 1 ;
iopt.luerr=luerr ;
end
%
% Checking dimensional parameters N,M and MFIT
if n <= 0 || m <= 0 || mfit < 0 || mfit > m
if mprerr >= 1
fprintf(luerr,'\nError: %s\n %s\n %s%5i%s%5i%s%5i\n', ...
'Bad or inconsistent input to dimensional parameters supplied', ...
'choose N and M positive, and MFIT <= M, MFIT nonnegative', ...
'your input is: N = ',n,' M = ',m,' MFIT = ',mfit) ;
end
hier1=n
hier2=m
hier3=mfit
ierr = 20 ;