QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
lmdif.cpp
Go to the documentation of this file.
1/************************lmdif*************************/
2
3/*
4The original Fortran version is Copyright (C) 1999 University of Chicago.
5All rights reserved.
6
7Redistribution and use in source and binary forms, with or
8without modification, are permitted provided that the
9following conditions are met:
10
111. Redistributions of source code must retain the above
12copyright notice, this list of conditions and the following
13disclaimer.
14
152. Redistributions in binary form must reproduce the above
16copyright notice, this list of conditions and the following
17disclaimer in the documentation and/or other materials
18provided with the distribution.
19
203. The end-user documentation included with the
21redistribution, if any, must include the following
22acknowledgment:
23
24 "This product includes software developed by the
25 University of Chicago, as Operator of Argonne National
26 Laboratory.
27
28Alternately, this acknowledgment may appear in the software
29itself, if and wherever such third-party acknowledgments
30normally appear.
31
324. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
33WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
34UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
35THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
36IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
37OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
38OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
39OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
40USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
41THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
42DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
43UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
44BE CORRECTED.
45
465. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
47HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
48ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
49INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
50ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
51PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
52SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
53(INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
54EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
55POSSIBILITY OF SUCH LOSS OR DAMAGES.
56
57
58C translation Copyright (C) Steve Moshier
59
60What you see here may be used freely but it comes with no support
61or guarantee.
62*/
63
65#include <cmath>
66#include <cstdio>
67
68namespace QuantLib {
69 namespace MINPACK {
70#define BUG 0
71/* resolution of arithmetic */
72double MACHEP = 1.2e-16;
73/* smallest nonzero number */
74double DWARF = 1.0e-38;
75
76
77
78
79
81{
82/*
83* **********
84*
85* function enorm
86*
87* given an n-vector x, this function calculates the
88* euclidean norm of x.
89*
90* the euclidean norm is computed by accumulating the sum of
91* squares in three different sums. the sums of squares for the
92* small and large components are scaled so that no overflows
93* occur. non-destructive underflows are permitted. underflows
94* and overflows do not occur in the computation of the unscaled
95* sum of squares for the intermediate components.
96* the definitions of small, intermediate and large components
97* depend on two constants, rdwarf and rgiant. the main
98* restrictions on these constants are that rdwarf**2 not
99* underflow and rgiant**2 not overflow. the constants
100* given here are suitable for every known computer.
101*
102* the function statement is
103*
104* double precision function enorm(n,x)
105*
106* where
107*
108* n is a positive integer input variable.
109*
110* x is an input array of length n.
111*
112* subprograms called
113*
114* fortran-supplied ... dabs,dsqrt
115*
116* argonne national laboratory. minpack project. march 1980.
117* burton s. garbow, kenneth e. hillstrom, jorge j. more
118*
119* **********
120*/
121int i;
122Real agiant,floatn,s1,s2,s3,xabs,x1max,x3max;
123Real ans, temp;
124static double rdwarf = 3.834e-20;
125static double rgiant = 1.304e19;
126static double zero = 0.0;
127static double one = 1.0;
128
129s1 = zero;
130s2 = zero;
131s3 = zero;
132x1max = zero;
133x3max = zero;
134floatn = n;
135agiant = rgiant/floatn;
136
137for( i=0; i<n; i++ )
138{
139xabs = std::fabs(x[i]);
140if( (xabs > rdwarf) && (xabs < agiant) )
141 {
142/*
143* sum for intermediate components.
144*/
145 s2 += xabs*xabs;
146 continue;
147 }
148
149if(xabs > rdwarf)
150 {
151/*
152* sum for large components.
153*/
154 if(xabs > x1max)
155 {
156 temp = x1max/xabs;
157 s1 = one + s1*temp*temp;
158 x1max = xabs;
159 }
160 else
161 {
162 temp = xabs/x1max;
163 s1 += temp*temp;
164 }
165 continue;
166 }
167/*
168* sum for small components.
169*/
170if(xabs > x3max)
171 {
172 temp = x3max/xabs;
173 s3 = one + s3*temp*temp;
174 x3max = xabs;
175 }
176else
177 {
178 if(xabs != zero)
179 {
180 temp = xabs/x3max;
181 s3 += temp*temp;
182 }
183 }
184}
185/*
186* calculation of norm.
187*/
188if(s1 != zero)
189 {
190 temp = s1 + (s2/x1max)/x1max;
191 ans = x1max*std::sqrt(temp);
192 return(ans);
193 }
194if(s2 != zero)
195 {
196 if(s2 >= x3max)
197 temp = s2*(one+(x3max/s2)*(x3max*s3));
198 else
199 temp = x3max*((s2/x3max)+(x3max*s3));
200 ans = std::sqrt(temp);
201 }
202else
203 {
204 ans = x3max*std::sqrt(s3);
205 }
206return(ans);
207/*
208* last card of function enorm.
209*/
210}
211/************************lmmisc.c*************************/
212
214{
215if( a >= b )
216 return(a);
217else
218 return(b);
219}
220
222{
223if( a <= b )
224 return(a);
225else
226 return(b);
227}
228
229int min0(int a,int b)
230
231{
232if( a <= b )
233 return(a);
234else
235 return(b);
236}
237
238int mod( int k, int m )
239{
240return( k % m );
241}
242
243
244/***********Sample of user supplied function****************
245 * m = number of functions
246 * n = number of variables
247 * x = vector of function arguments
248 * fvec = vector of function values
249 * iflag = error return variable
250 */
251//void fcn(int m,int n, Real* x, Real* fvec,int *iflag)
252//{
253// QuantLib::LevenbergMarquardt::fcn(m, n, x, fvec, iflag);
254//}
255
256void fdjac2(int m,
257 int n,
258 Real* x,
259 const Real* fvec,
260 Real* fjac,
261 int,
262 int* iflag,
263 Real epsfcn,
264 Real* wa,
266 /*
267 * **********
268 *
269 * subroutine fdjac2
270 *
271 * this subroutine computes a forward-difference approximation
272 * to the m by n jacobian matrix associated with a specified
273 * problem of m functions in n variables.
274 *
275 * the subroutine statement is
276 *
277 * subroutine fdjac2(fcn,m,n,x,fvec,fjac,ldfjac,iflag,epsfcn,wa)
278 *
279 * where
280 *
281 * fcn is the name of the user-supplied subroutine which
282 * calculates the functions. fcn must be declared
283 * in an external statement in the user calling
284 * program, and should be written as follows.
285 *
286 * subroutine fcn(m,n,x,fvec,iflag)
287 * integer m,n,iflag
288 * double precision x(n),fvec(m)
289 * ----------
290 * calculate the functions at x and
291 * return this vector in fvec.
292 * ----------
293 * return
294 * end
295 *
296 * the value of iflag should not be changed by fcn unless
297 * the user wants to terminate execution of fdjac2.
298 * in this case set iflag to a negative integer.
299 *
300 * m is a positive integer input variable set to the number
301 * of functions.
302 *
303 * n is a positive integer input variable set to the number
304 * of variables. n must not exceed m.
305 *
306 * x is an input array of length n.
307 *
308 * fvec is an input array of length m which must contain the
309 * functions evaluated at x.
310 *
311 * fjac is an output m by n array which contains the
312 * approximation to the jacobian matrix evaluated at x.
313 *
314 * ldfjac is a positive integer input variable not less than m
315 * which specifies the leading dimension of the array fjac.
316 *
317 * iflag is an integer variable which can be used to terminate
318 * the execution of fdjac2. see description of fcn.
319 *
320 * epsfcn is an input variable used in determining a suitable
321 * step length for the forward-difference approximation. this
322 * approximation assumes that the relative errors in the
323 * functions are of the order of epsfcn. if epsfcn is less
324 * than the machine precision, it is assumed that the relative
325 * errors in the functions are of the order of the machine
326 * precision.
327 *
328 * wa is a work array of length m.
329 *
330 * subprograms called
331 *
332 * user-supplied ...... fcn
333 *
334 * minpack-supplied ... dpmpar
335 *
336 * fortran-supplied ... dabs,dmax1,dsqrt
337 *
338 * argonne national laboratory. minpack project. march 1980.
339 * burton s. garbow, kenneth e. hillstrom, jorge j. more
340 *
341 **********
342 */
343 int i, j, ij;
344 Real eps, h, temp;
345 static double zero = 0.0;
346
347
348 temp = dmax1(epsfcn, MACHEP);
349 eps = std::sqrt(temp);
350 ij = 0;
351 for (j = 0; j < n; j++) {
352 temp = x[j];
353 h = eps * std::fabs(temp);
354 if (h == zero)
355 h = eps;
356 x[j] = temp + h;
357 fcn(m, n, x, wa, iflag);
358 if (*iflag < 0)
359 return;
360 x[j] = temp;
361 for (i = 0; i < m; i++) {
362 fjac[ij] = (wa[i] - fvec[i]) / h;
363 ij += 1; /* fjac[i+m*j] */
364 }
365 }
366/*
367* last card of subroutine fdjac2.
368*/
369}
370/************************qrfac.c*************************/
371
372
373void
374qrfac(int m,int n,Real* a,int,int pivot,int* ipvt,
375 int,Real* rdiag,Real* acnorm,Real* wa)
376{
377/*
378* **********
379*
380* subroutine qrfac
381*
382* this subroutine uses householder transformations with column
383* pivoting (optional) to compute a qr factorization of the
384* m by n matrix a. that is, qrfac determines an orthogonal
385* matrix q, a permutation matrix p, and an upper trapezoidal
386* matrix r with diagonal elements of nonincreasing magnitude,
387* such that a*p = q*r. the householder transformation for
388* column k, k = 1,2,...,min(m,n), is of the form
389*
390* t
391* i - (1/u(k))*u*u
392*
393* where u has zeros in the first k-1 positions. the form of
394* this transformation and the method of pivoting first
395* appeared in the corresponding linpack subroutine.
396*
397* the subroutine statement is
398*
399* subroutine qrfac(m,n,a,lda,pivot,ipvt,lipvt,rdiag,acnorm,wa)
400*
401* where
402*
403* m is a positive integer input variable set to the number
404* of rows of a.
405*
406* n is a positive integer input variable set to the number
407* of columns of a.
408*
409* a is an m by n array. on input a contains the matrix for
410* which the qr factorization is to be computed. on output
411* the strict upper trapezoidal part of a contains the strict
412* upper trapezoidal part of r, and the lower trapezoidal
413* part of a contains a factored form of q (the non-trivial
414* elements of the u vectors described above).
415*
416* lda is a positive integer input variable not less than m
417* which specifies the leading dimension of the array a.
418*
419* pivot is a logical input variable. if pivot is set true,
420* then column pivoting is enforced. if pivot is set false,
421* then no column pivoting is done.
422*
423* ipvt is an integer output array of length lipvt. ipvt
424* defines the permutation matrix p such that a*p = q*r.
425* column j of p is column ipvt(j) of the identity matrix.
426* if pivot is false, ipvt is not referenced.
427*
428* lipvt is a positive integer input variable. if pivot is false,
429* then lipvt may be as small as 1. if pivot is true, then
430* lipvt must be at least n.
431*
432* rdiag is an output array of length n which contains the
433* diagonal elements of r.
434*
435* acnorm is an output array of length n which contains the
436* norms of the corresponding columns of the input matrix a.
437* if this information is not needed, then acnorm can coincide
438* with rdiag.
439*
440* wa is a work array of length n. if pivot is false, then wa
441* can coincide with rdiag.
442*
443* subprograms called
444*
445* minpack-supplied ... dpmpar,enorm
446*
447* fortran-supplied ... dmax1,dsqrt,min0
448*
449* argonne national laboratory. minpack project. march 1980.
450* burton s. garbow, kenneth e. hillstrom, jorge j. more
451*
452* **********
453*/
454int i,ij,jj,j,jp1,k,kmax,minmn;
455Real ajnorm,sum,temp;
456static double zero = 0.0;
457static double one = 1.0;
458static double p05 = 0.05;
459
460/*
461* compute the initial column norms and initialize several arrays.
462*/
463ij = 0;
464for( j=0; j<n; j++ )
465 {
466 acnorm[j] = enorm(m,&a[ij]);
467 rdiag[j] = acnorm[j];
468 wa[j] = rdiag[j];
469 if(pivot != 0)
470 ipvt[j] = j;
471 ij += m; /* m*j */
472 }
473/*
474* reduce a to r with householder transformations.
475*/
476minmn = min0(m,n);
477for( j=0; j<minmn; j++ )
478{
479if(pivot == 0)
480 goto L40;
481/*
482* bring the column of largest norm into the pivot position.
483*/
484kmax = j;
485for( k=j; k<n; k++ )
486 {
487 if(rdiag[k] > rdiag[kmax])
488 kmax = k;
489 }
490if(kmax == j)
491 goto L40;
492
493ij = m * j;
494jj = m * kmax;
495for( i=0; i<m; i++ )
496 {
497 temp = a[ij]; /* [i+m*j] */
498 a[ij] = a[jj]; /* [i+m*kmax] */
499 a[jj] = temp;
500 ij += 1;
501 jj += 1;
502 }
503rdiag[kmax] = rdiag[j];
504wa[kmax] = wa[j];
505k = ipvt[j];
506ipvt[j] = ipvt[kmax];
507ipvt[kmax] = k;
508
509L40:
510/*
511* compute the householder transformation to reduce the
512* j-th column of a to a multiple of the j-th unit vector.
513*/
514jj = j + m*j;
515ajnorm = enorm(m-j,&a[jj]);
516if(ajnorm == zero)
517 goto L100;
518if(a[jj] < zero)
519 ajnorm = -ajnorm;
520ij = jj;
521for( i=j; i<m; i++ )
522 {
523 a[ij] /= ajnorm;
524 ij += 1; /* [i+m*j] */
525 }
526a[jj] += one;
527/*
528* apply the transformation to the remaining columns
529* and update the norms.
530*/
531jp1 = j + 1;
532if(jp1 < n )
533{
534for( k=jp1; k<n; k++ )
535 {
536 sum = zero;
537 ij = j + m*k;
538 jj = j + m*j;
539 for( i=j; i<m; i++ )
540 {
541 sum += a[jj]*a[ij];
542 ij += 1; /* [i+m*k] */
543 jj += 1; /* [i+m*j] */
544 }
545 temp = sum/a[j+m*j];
546 ij = j + m*k;
547 jj = j + m*j;
548 for( i=j; i<m; i++ )
549 {
550 a[ij] -= temp*a[jj];
551 ij += 1; /* [i+m*k] */
552 jj += 1; /* [i+m*j] */
553 }
554 if( (pivot != 0) && (rdiag[k] != zero) )
555 {
556 temp = a[j+m*k]/rdiag[k];
557 temp = dmax1( zero, one-temp*temp );
558 rdiag[k] *= std::sqrt(temp);
559 temp = rdiag[k]/wa[k];
560 if( (p05*temp*temp) <= MACHEP)
561 {
562 rdiag[k] = enorm(m-j-1,&a[jp1+m*k]);
563 wa[k] = rdiag[k];
564 }
565 }
566 }
567}
568
569L100:
570 rdiag[j] = -ajnorm;
571}
572/*
573* last card of subroutine qrfac.
574*/
575}
576
577/************************qrsolv.c*************************/
578
579
580void qrsolv(int n,
581 Real* r,
582 int ldr,
583 const int* ipvt,
584 const Real* diag,
585 const Real* qtb,
586 Real* x,
587 Real* sdiag,
588 Real* wa) {
589 /*
590 * **********
591 *
592 * subroutine qrsolv
593 *
594 * given an m by n matrix a, an n by n diagonal matrix d,
595 * and an m-vector b, the problem is to determine an x which
596 * solves the system
597 *
598 * a*x = b , d*x = 0 ,
599 *
600 * in the least squares sense.
601 *
602 * this subroutine completes the solution of the problem
603 * if it is provided with the necessary information from the
604 * qr factorization, with column pivoting, of a. that is, if
605 * a*p = q*r, where p is a permutation matrix, q has orthogonal
606 * columns, and r is an upper triangular matrix with diagonal
607 * elements of nonincreasing magnitude, then qrsolv expects
608 * the full upper triangle of r, the permutation matrix p,
609 * and the first n components of (q transpose)*b. the system
610 * a*x = b, d*x = 0, is then equivalent to
611 *
612 * t t
613 * r*z = q *b , p *d*p*z = 0 ,
614 *
615 * where x = p*z. if this system does not have full rank,
616 * then a least squares solution is obtained. on output qrsolv
617 * also provides an upper triangular matrix s such that
618 *
619 * t t t
620 * p *(a *a + d*d)*p = s *s .
621 *
622 * s is computed within qrsolv and may be of separate interest.
623 *
624 * the subroutine statement is
625 *
626 * subroutine qrsolv(n,r,ldr,ipvt,diag,qtb,x,sdiag,wa)
627 *
628 * where
629 *
630 * n is a positive integer input variable set to the order of r.
631 *
632 * r is an n by n array. on input the full upper triangle
633 * must contain the full upper triangle of the matrix r.
634 * on output the full upper triangle is unaltered, and the
635 * strict lower triangle contains the strict upper triangle
636 * (transposed) of the upper triangular matrix s.
637 *
638 * ldr is a positive integer input variable not less than n
639 * which specifies the leading dimension of the array r.
640 *
641 * ipvt is an integer input array of length n which defines the
642 * permutation matrix p such that a*p = q*r. column j of p
643 * is column ipvt(j) of the identity matrix.
644 *
645 * diag is an input array of length n which must contain the
646 * diagonal elements of the matrix d.
647 *
648 * qtb is an input array of length n which must contain the first
649 * n elements of the vector (q transpose)*b.
650 *
651 * x is an output array of length n which contains the least
652 * squares solution of the system a*x = b, d*x = 0.
653 *
654 * sdiag is an output array of length n which contains the
655 * diagonal elements of the upper triangular matrix s.
656 *
657 * wa is a work array of length n.
658 *
659 * subprograms called
660 *
661 * fortran-supplied ... dabs,dsqrt
662 *
663 * argonne national laboratory. minpack project. march 1980.
664 * burton s. garbow, kenneth e. hillstrom, jorge j. more
665 *
666 * **********
667 */
668 int i, ij, ik, kk, j, jp1, k, kp1, l, nsing;
669 Real cos, cotan, qtbpj, sin, sum, tan, temp;
670 static double zero = 0.0;
671 static double p25 = 0.25;
672 static double p5 = 0.5;
673
674 /*
675 * copy r and (q transpose)*b to preserve input and initialize s.
676 * in particular, save the diagonal elements of r in x.
677 */
678 kk = 0;
679 for (j = 0; j < n; j++) {
680 ij = kk;
681 ik = kk;
682 for (i = j; i < n; i++) {
683 r[ij] = r[ik];
684 ij += 1; /* [i+ldr*j] */
685 ik += ldr; /* [j+ldr*i] */
686 }
687 x[j] = r[kk];
688 wa[j] = qtb[j];
689 kk += ldr+1; /* j+ldr*j */
690 }
691/*
692* eliminate the diagonal matrix d using a givens rotation.
693*/
694for( j=0; j<n; j++ )
695{
696/*
697* prepare the row of d to be eliminated, locating the
698* diagonal element using p from the qr factorization.
699*/
700l = ipvt[j];
701if(diag[l] == zero)
702 goto L90;
703for( k=j; k<n; k++ )
704 sdiag[k] = zero;
705sdiag[j] = diag[l];
706/*
707* the transformations to eliminate the row of d
708* modify only a single element of (q transpose)*b
709* beyond the first n, which is initially zero.
710*/
711qtbpj = zero;
712for( k=j; k<n; k++ )
713 {
714/*
715* determine a givens rotation which eliminates the
716* appropriate element in the current row of d.
717*/
718 if(sdiag[k] == zero)
719 continue;
720 kk = k + ldr * k;
721 if(std::fabs(r[kk]) < std::fabs(sdiag[k]))
722 {
723 cotan = r[kk]/sdiag[k];
724 sin = p5/std::sqrt(p25+p25*cotan*cotan);
725 cos = sin*cotan;
726 }
727 else
728 {
729 tan = sdiag[k]/r[kk];
730 cos = p5/std::sqrt(p25+p25*tan*tan);
731 sin = cos*tan;
732 }
733/*
734* compute the modified diagonal element of r and
735* the modified element of ((q transpose)*b,0).
736*/
737 r[kk] = cos*r[kk] + sin*sdiag[k];
738 temp = cos*wa[k] + sin*qtbpj;
739 qtbpj = -sin*wa[k] + cos*qtbpj;
740 wa[k] = temp;
741/*
742* accumulate the tranformation in the row of s.
743*/
744 kp1 = k + 1;
745 if( n > kp1 )
746 {
747 ik = kk + 1;
748 for( i=kp1; i<n; i++ )
749 {
750 temp = cos*r[ik] + sin*sdiag[i];
751 sdiag[i] = -sin*r[ik] + cos*sdiag[i];
752 r[ik] = temp;
753 ik += 1; /* [i+ldr*k] */
754 }
755 }
756 }
757L90:
758/*
759* store the diagonal element of s and restore
760* the corresponding diagonal element of r.
761*/
762 kk = j + ldr*j;
763 sdiag[j] = r[kk];
764 r[kk] = x[j];
765}
766/*
767* solve the triangular system for z. if the system is
768* singular, then obtain a least squares solution.
769*/
770nsing = n;
771for( j=0; j<n; j++ )
772 {
773 if( (sdiag[j] == zero) && (nsing == n) )
774 nsing = j;
775 if(nsing < n)
776 wa[j] = zero;
777 }
778if(nsing < 1)
779 goto L150;
780
781for( k=0; k<nsing; k++ )
782 {
783 j = nsing - k - 1;
784 sum = zero;
785 jp1 = j + 1;
786 if(nsing > jp1)
787 {
788 ij = jp1 + ldr * j;
789 for( i=jp1; i<nsing; i++ )
790 {
791 sum += r[ij]*wa[i];
792 ij += 1; /* [i+ldr*j] */
793 }
794 }
795 wa[j] = (wa[j] - sum)/sdiag[j];
796 }
797L150:
798/*
799* permute the components of z back to components of x.
800*/
801for( j=0; j<n; j++ )
802 {
803 l = ipvt[j];
804 x[l] = wa[j];
805 }
806/*
807* last card of subroutine qrsolv.
808*/
809}
810
811/************************lmpar.c*************************/
812
813
814void lmpar(int n,
815 Real* r,
816 int ldr,
817 int* ipvt,
818 const Real* diag,
819 Real* qtb,
820 Real delta,
821 Real* par,
822 Real* x,
823 Real* sdiag,
824 Real* wa1,
825 Real* wa2) {
826 /* **********
827 *
828 * subroutine lmpar
829 *
830 * given an m by n matrix a, an n by n nonsingular diagonal
831 * matrix d, an m-vector b, and a positive number delta,
832 * the problem is to determine a value for the parameter
833 * par such that if x solves the system
834 *
835 * a*x = b , sqrt(par)*d*x = 0 ,
836 *
837 * in the least squares sense, and dxnorm is the euclidean
838 * norm of d*x, then either par is zero and
839 *
840 * (dxnorm-delta) .le. 0.1*delta ,
841 *
842 * or par is positive and
843 *
844 * abs(dxnorm-delta) .le. 0.1*delta .
845 *
846 * this subroutine completes the solution of the problem
847 * if it is provided with the necessary information from the
848 * qr factorization, with column pivoting, of a. that is, if
849 * a*p = q*r, where p is a permutation matrix, q has orthogonal
850 * columns, and r is an upper triangular matrix with diagonal
851 * elements of nonincreasing magnitude, then lmpar expects
852 * the full upper triangle of r, the permutation matrix p,
853 * and the first n components of (q transpose)*b. on output
854 * lmpar also provides an upper triangular matrix s such that
855 *
856 * t t t
857 * p *(a *a + par*d*d)*p = s *s .
858 *
859 * s is employed within lmpar and may be of separate interest.
860 *
861 * only a few iterations are generally needed for convergence
862 * of the algorithm. if, however, the limit of 10 iterations
863 * is reached, then the output par will contain the best
864 * value obtained so far.
865 *
866 * the subroutine statement is
867 *
868 * subroutine lmpar(n,r,ldr,ipvt,diag,qtb,delta,par,x,sdiag,
869 * wa1,wa2)
870 *
871 * where
872 *
873 * n is a positive integer input variable set to the order of r.
874 *
875 * r is an n by n array. on input the full upper triangle
876 * must contain the full upper triangle of the matrix r.
877 * on output the full upper triangle is unaltered, and the
878 * strict lower triangle contains the strict upper triangle
879 * (transposed) of the upper triangular matrix s.
880 *
881 * ldr is a positive integer input variable not less than n
882 * which specifies the leading dimension of the array r.
883 *
884 * ipvt is an integer input array of length n which defines the
885 * permutation matrix p such that a*p = q*r. column j of p
886 * is column ipvt(j) of the identity matrix.
887 *
888 * diag is an input array of length n which must contain the
889 * diagonal elements of the matrix d.
890 *
891 * qtb is an input array of length n which must contain the first
892 * n elements of the vector (q transpose)*b.
893 *
894 * delta is a positive input variable which specifies an upper
895 * bound on the euclidean norm of d*x.
896 *
897 * par is a nonnegative variable. on input par contains an
898 * initial estimate of the levenberg-marquardt parameter.
899 * on output par contains the final estimate.
900 *
901 * x is an output array of length n which contains the least
902 * squares solution of the system a*x = b, sqrt(par)*d*x = 0,
903 * for the output par.
904 *
905 * sdiag is an output array of length n which contains the
906 * diagonal elements of the upper triangular matrix s.
907 *
908 * wa1 and wa2 are work arrays of length n.
909 *
910 * subprograms called
911 *
912 * minpack-supplied ... dpmpar,enorm,qrsolv
913 *
914 * fortran-supplied ... dabs,dmax1,dmin1,dsqrt
915 *
916 * argonne national laboratory. minpack project. march 1980.
917 * burton s. garbow, kenneth e. hillstrom, jorge j. more
918 *
919 * **********
920 */
921 int i, iter, ij, jj, j, jm1, jp1, k, l, nsing;
922 Real dxnorm, fp, gnorm, parc, parl, paru;
923 Real sum, temp;
924 static double zero = 0.0;
925 static double p1 = 0.1;
926 static double p001 = 0.001;
927
928
929 /*
930 * compute and store in x the gauss-newton direction. if the
931 * jacobian is rank-deficient, obtain a least squares solution.
932 */
933 nsing = n;
934 jj = 0;
935 for (j = 0; j < n; j++) {
936 wa1[j] = qtb[j];
937 if ((r[jj] == zero) && (nsing == n))
938 nsing = j;
939 if (nsing < n)
940 wa1[j] = zero;
941 jj += ldr + 1; /* [j+ldr*j] */
942 }
943if(nsing >= 1)
944 {
945 for( k=0; k<nsing; k++ )
946 {
947 j = nsing - k - 1;
948 wa1[j] = wa1[j]/r[j+ldr*j];
949 temp = wa1[j];
950 jm1 = j - 1;
951 if(jm1 >= 0)
952 {
953 ij = ldr * j;
954 for( i=0; i<=jm1; i++ )
955 {
956 wa1[i] -= r[ij]*temp;
957 ij += 1;
958 }
959 }
960 }
961 }
962
963for( j=0; j<n; j++ )
964 {
965 l = ipvt[j];
966 x[l] = wa1[j];
967 }
968/*
969* initialize the iteration counter.
970* evaluate the function at the origin, and test
971* for acceptance of the gauss-newton direction.
972*/
973iter = 0;
974for( j=0; j<n; j++ )
975 wa2[j] = diag[j]*x[j];
976dxnorm = enorm(n,wa2);
977fp = dxnorm - delta;
978if(fp <= p1*delta)
979 {
980 goto L220;
981 }
982/*
983* if the jacobian is not rank deficient, the newton
984* step provides a lower bound, parl, for the zero of
985* the function. otherwise set this bound to zero.
986*/
987parl = zero;
988if(nsing >= n)
989 {
990 for( j=0; j<n; j++ )
991 {
992 l = ipvt[j];
993 wa1[j] = diag[l]*(wa2[l]/dxnorm);
994 }
995 jj = 0;
996 for( j=0; j<n; j++ )
997 {
998 sum = zero;
999 jm1 = j - 1;
1000 if(jm1 >= 0)
1001 {
1002 ij = jj;
1003 for( i=0; i<=jm1; i++ )
1004 {
1005 sum += r[ij]*wa1[i];
1006 ij += 1;
1007 }
1008 }
1009 wa1[j] = (wa1[j] - sum)/r[j+ldr*j];
1010 jj += ldr; /* [i+ldr*j] */
1011 }
1012 temp = enorm(n,wa1);
1013 parl = ((fp/delta)/temp)/temp;
1014 }
1015/*
1016* calculate an upper bound, paru, for the zero of the function.
1017*/
1018jj = 0;
1019for( j=0; j<n; j++ )
1020 {
1021 sum = zero;
1022 ij = jj;
1023 for( i=0; i<=j; i++ )
1024 {
1025 sum += r[ij]*qtb[i];
1026 ij += 1;
1027 }
1028 l = ipvt[j];
1029 wa1[j] = sum/diag[l];
1030 jj += ldr; /* [i+ldr*j] */
1031 }
1032gnorm = enorm(n,wa1);
1033paru = gnorm/delta;
1034if(paru == zero)
1035 paru = DWARF/dmin1(delta,p1);
1036/*
1037* if the input par lies outside of the interval (parl,paru),
1038* set par to the closer endpoint.
1039*/
1040*par = dmax1( *par,parl);
1041*par = dmin1( *par,paru);
1042if( *par == zero)
1043 *par = gnorm/dxnorm;
1044/*
1045* beginning of an iteration.
1046*/
1047L150:
1048iter += 1;
1049/*
1050* evaluate the function at the current value of par.
1051*/
1052if( *par == zero)
1053 *par = dmax1(DWARF,p001*paru);
1054temp = std::sqrt( *par );
1055for( j=0; j<n; j++ )
1056 wa1[j] = temp*diag[j];
1057qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2);
1058for( j=0; j<n; j++ )
1059 wa2[j] = diag[j]*x[j];
1060dxnorm = enorm(n,wa2);
1061temp = fp;
1062fp = dxnorm - delta;
1063/*
1064* if the function is small enough, accept the current value
1065* of par. also test for the exceptional cases where parl
1066* is zero or the number of iterations has reached 10.
1067*/
1068if( (std::fabs(fp) <= p1*delta)
1069 || ((parl == zero) && (fp <= temp) && (temp < zero))
1070 || (iter == 10) )
1071 goto L220;
1072/*
1073* compute the newton correction.
1074*/
1075for( j=0; j<n; j++ )
1076 {
1077 l = ipvt[j];
1078 wa1[j] = diag[l]*(wa2[l]/dxnorm);
1079 }
1080jj = 0;
1081for( j=0; j<n; j++ )
1082 {
1083 wa1[j] = wa1[j]/sdiag[j];
1084 temp = wa1[j];
1085 jp1 = j + 1;
1086 if(jp1 < n)
1087 {
1088 ij = jp1 + jj;
1089 for( i=jp1; i<n; i++ )
1090 {
1091 wa1[i] -= r[ij]*temp;
1092 ij += 1; /* [i+ldr*j] */
1093 }
1094 }
1095 jj += ldr; /* ldr*j */
1096 }
1097temp = enorm(n,wa1);
1098parc = ((fp/delta)/temp)/temp;
1099/*
1100* depending on the sign of the function, update parl or paru.
1101*/
1102if(fp > zero)
1103 parl = dmax1(parl, *par);
1104if(fp < zero)
1105 paru = dmin1(paru, *par);
1106/*
1107* compute an improved estimate for par.
1108*/
1109*par = dmax1(parl, *par + parc);
1110/*
1111* end of an iteration.
1112*/
1113goto L150;
1114
1115L220:
1116/*
1117* termination.
1118*/
1119if(iter == 0)
1120 *par = zero;
1121/*
1122* last card of subroutine lmpar.
1123*/
1124}
1125
1126/*********************** lmdif.c ****************************/
1127
1128
1129
1130
1131void lmdif(int m,int n,Real* x,Real* fvec,Real ftol,
1132 Real xtol,Real gtol,int maxfev,Real epsfcn,
1133 Real* diag, int mode, Real factor,
1134 int nprint, int* info,int* nfev,Real* fjac,
1135 int ldfjac,int* ipvt,Real* qtf,
1136 Real* wa1,Real* wa2,Real* wa3,Real* wa4,
1139{
1140/*
1141* **********
1142*
1143* subroutine lmdif
1144*
1145* the purpose of lmdif is to minimize the sum of the squares of
1146* m nonlinear functions in n variables by a modification of
1147* the levenberg-marquardt algorithm. the user must provide a
1148* subroutine which calculates the functions. the jacobian is
1149* then calculated by a forward-difference approximation.
1150*
1151* the subroutine statement is
1152*
1153* subroutine lmdif(fcn,m,n,x,fvec,ftol,xtol,gtol,maxfev,epsfcn,
1154* diag,mode,factor,nprint,info,nfev,fjac,
1155* ldfjac,ipvt,qtf,wa1,wa2,wa3,wa4)
1156*
1157* where
1158*
1159* fcn is the name of the user-supplied subroutine which
1160* calculates the functions. fcn must be declared
1161* in an external statement in the user calling
1162* program, and should be written as follows.
1163*
1164* subroutine fcn(m,n,x,fvec,iflag)
1165* integer m,n,iflag
1166* double precision x(n),fvec(m)
1167* ----------
1168* calculate the functions at x and
1169* return this vector in fvec.
1170* ----------
1171* return
1172* end
1173*
1174* the value of iflag should not be changed by fcn unless
1175* the user wants to terminate execution of lmdif.
1176* in this case set iflag to a negative integer.
1177*
1178* m is a positive integer input variable set to the number
1179* of functions.
1180*
1181* n is a positive integer input variable set to the number
1182* of variables. n must not exceed m.
1183*
1184* x is an array of length n. on input x must contain
1185* an initial estimate of the solution vector. on output x
1186* contains the final estimate of the solution vector.
1187*
1188* fvec is an output array of length m which contains
1189* the functions evaluated at the output x.
1190*
1191* ftol is a nonnegative input variable. termination
1192* occurs when both the actual and predicted relative
1193* reductions in the sum of squares are at most ftol.
1194* therefore, ftol measures the relative error desired
1195* in the sum of squares.
1196*
1197* xtol is a nonnegative input variable. termination
1198* occurs when the relative error between two consecutive
1199* iterates is at most xtol. therefore, xtol measures the
1200* relative error desired in the approximate solution.
1201*
1202* gtol is a nonnegative input variable. termination
1203* occurs when the cosine of the angle between fvec and
1204* any column of the jacobian is at most gtol in absolute
1205* value. therefore, gtol measures the orthogonality
1206* desired between the function vector and the columns
1207* of the jacobian.
1208*
1209* maxfev is a positive integer input variable. termination
1210* occurs when the number of calls to fcn is at least
1211* maxfev by the end of an iteration.
1212*
1213* epsfcn is an input variable used in determining a suitable
1214* step length for the forward-difference approximation. this
1215* approximation assumes that the relative errors in the
1216* functions are of the order of epsfcn. if epsfcn is less
1217* than the machine precision, it is assumed that the relative
1218* errors in the functions are of the order of the machine
1219* precision.
1220*
1221* diag is an array of length n. if mode = 1 (see
1222* below), diag is internally set. if mode = 2, diag
1223* must contain positive entries that serve as
1224* multiplicative scale factors for the variables.
1225*
1226* mode is an integer input variable. if mode = 1, the
1227* variables will be scaled internally. if mode = 2,
1228* the scaling is specified by the input diag. other
1229* values of mode are equivalent to mode = 1.
1230*
1231* factor is a positive input variable used in determining the
1232* initial step bound. this bound is set to the product of
1233* factor and the euclidean norm of diag*x if nonzero, or else
1234* to factor itself. in most cases factor should lie in the
1235* interval (.1,100.). 100. is a generally recommended value.
1236*
1237* nprint is an integer input variable that enables controlled
1238* printing of iterates if it is positive. in this case,
1239* fcn is called with iflag = 0 at the beginning of the first
1240* iteration and every nprint iterations thereafter and
1241* immediately prior to return, with x and fvec available
1242* for printing. if nprint is not positive, no special calls
1243* of fcn with iflag = 0 are made.
1244*
1245* info is an integer output variable. if the user has
1246* terminated execution, info is set to the (negative)
1247* value of iflag. see description of fcn. otherwise,
1248* info is set as follows.
1249*
1250* info = 0 improper input parameters.
1251*
1252* info = 1 both actual and predicted relative reductions
1253* in the sum of squares are at most ftol.
1254*
1255* info = 2 relative error between two consecutive iterates
1256* is at most xtol.
1257*
1258* info = 3 conditions for info = 1 and info = 2 both hold.
1259*
1260* info = 4 the cosine of the angle between fvec and any
1261* column of the jacobian is at most gtol in
1262* absolute value.
1263*
1264* info = 5 number of calls to fcn has reached or
1265* exceeded maxfev.
1266*
1267* info = 6 ftol is too small. no further reduction in
1268* the sum of squares is possible.
1269*
1270* info = 7 xtol is too small. no further improvement in
1271* the approximate solution x is possible.
1272*
1273* info = 8 gtol is too small. fvec is orthogonal to the
1274* columns of the jacobian to machine precision.
1275*
1276* nfev is an integer output variable set to the number of
1277* calls to fcn.
1278*
1279* fjac is an output m by n array. the upper n by n submatrix
1280* of fjac contains an upper triangular matrix r with
1281* diagonal elements of nonincreasing magnitude such that
1282*
1283* t t t
1284* p *(jac *jac)*p = r *r,
1285*
1286* where p is a permutation matrix and jac is the final
1287* calculated jacobian. column j of p is column ipvt(j)
1288* (see below) of the identity matrix. the lower trapezoidal
1289* part of fjac contains information generated during
1290* the computation of r.
1291*
1292* ldfjac is a positive integer input variable not less than m
1293* which specifies the leading dimension of the array fjac.
1294*
1295* ipvt is an integer output array of length n. ipvt
1296* defines a permutation matrix p such that jac*p = q*r,
1297* where jac is the final calculated jacobian, q is
1298* orthogonal (not stored), and r is upper triangular
1299* with diagonal elements of nonincreasing magnitude.
1300* column j of p is column ipvt(j) of the identity matrix.
1301*
1302* qtf is an output array of length n which contains
1303* the first n elements of the vector (q transpose)*fvec.
1304*
1305* wa1, wa2, and wa3 are work arrays of length n.
1306*
1307* wa4 is a work array of length m.
1308*
1309* subprograms called
1310*
1311* user-supplied ...... fcn, jacFcn
1312*
1313* minpack-supplied ... dpmpar,enorm,fdjac2,lmpar,qrfac
1314*
1315* fortran-supplied ... dabs,dmax1,dmin1,dsqrt,mod
1316*
1317* argonne national laboratory. minpack project. march 1980.
1318* burton s. garbow, kenneth e. hillstrom, jorge j. more
1319*
1320* **********
1321*/
1322int i,iflag,ij,jj,iter,j,l;
1323Real actred,delta=0,dirder,fnorm,fnorm1,gnorm;
1324Real par,pnorm,prered,ratio;
1325Real sum,temp,temp1,temp2,temp3,xnorm=0;
1326static double one = 1.0;
1327static double p1 = 0.1;
1328static double p5 = 0.5;
1329static double p25 = 0.25;
1330static double p75 = 0.75;
1331static double p0001 = 1.0e-4;
1332static double zero = 0.0;
1333
1334*info = 0;
1335iflag = 0;
1336*nfev = 0;
1337/*
1338* check the input parameters for errors.
1339*/
1340if( (n <= 0) || (m < n) || (ldfjac < m) || (ftol < zero)
1341 || (xtol < zero) || (gtol < zero) || (maxfev <= 0)
1342 || (factor <= zero) )
1343 goto L300;
1344
1345if( mode == 2 )
1346 { /* scaling by diag[] */
1347 for( j=0; j<n; j++ )
1348 {
1349 if( diag[j] <= 0.0 )
1350 goto L300;
1351 }
1352 }
1353/*
1354* evaluate the function at the starting point
1355* and calculate its norm.
1356*/
1357iflag = 1;
1358fcn(m,n,x,fvec,&iflag);
1359*nfev = 1;
1360if(iflag < 0)
1361 goto L300;
1362fnorm = enorm(m,fvec);
1363/*
1364* initialize levenberg-marquardt parameter and iteration counter.
1365*/
1366par = zero;
1367iter = 1;
1368/*
1369* beginning of the outer loop.
1370*/
1371
1372L30:
1373
1374/*
1375* calculate the jacobian matrix.
1376*/
1377iflag = 2;
1378if (!jacFcn)
1379 fdjac2(m,n,x,fvec,fjac,ldfjac,&iflag,epsfcn,wa4, fcn);
1380else // use user supplied jacobian calculation
1381 jacFcn(m,n,x,fjac,&iflag);
1382*nfev += n;
1383if(iflag < 0)
1384 goto L300;
1385/*
1386* if requested, call fcn to enable printing of iterates.
1387*/
1388if( nprint > 0 )
1389 {
1390 iflag = 0;
1391 if(mod(iter-1,nprint) == 0)
1392 {
1393 fcn(m,n,x,fvec,&iflag);
1394 if(iflag < 0)
1395 goto L300;
1396 }
1397 }
1398/*
1399* compute the qr factorization of the jacobian.
1400*/
1401qrfac(m,n,fjac,ldfjac,1,ipvt,n,wa1,wa2,wa3);
1402/*
1403* on the first iteration and if mode is 1, scale according
1404* to the norms of the columns of the initial jacobian.
1405*/
1406if(iter == 1)
1407 {
1408 if(mode != 2)
1409 {
1410 for( j=0; j<n; j++ )
1411 {
1412 diag[j] = wa2[j];
1413 if( wa2[j] == zero )
1414 diag[j] = one;
1415 }
1416 }
1417
1418/*
1419* on the first iteration, calculate the norm of the scaled x
1420* and initialize the step bound delta.
1421*/
1422 for( j=0; j<n; j++ )
1423 wa3[j] = diag[j] * x[j];
1424
1425 xnorm = enorm(n,wa3);
1426 delta = factor*xnorm;
1427 if(delta == zero)
1428 delta = factor;
1429 }
1430
1431/*
1432* form (q transpose)*fvec and store the first n components in
1433* qtf.
1434*/
1435for( i=0; i<m; i++ )
1436 wa4[i] = fvec[i];
1437jj = 0;
1438for( j=0; j<n; j++ )
1439 {
1440 temp3 = fjac[jj];
1441 if(temp3 != zero)
1442 {
1443 sum = zero;
1444 ij = jj;
1445 for( i=j; i<m; i++ )
1446 {
1447 sum += fjac[ij] * wa4[i];
1448 ij += 1; /* fjac[i+m*j] */
1449 }
1450 temp = -sum / temp3;
1451 ij = jj;
1452 for( i=j; i<m; i++ )
1453 {
1454 wa4[i] += fjac[ij] * temp;
1455 ij += 1; /* fjac[i+m*j] */
1456 }
1457 }
1458 fjac[jj] = wa1[j];
1459 jj += m+1; /* fjac[j+m*j] */
1460 qtf[j] = wa4[j];
1461 }
1462
1463/*
1464* compute the norm of the scaled gradient.
1465*/
1466 gnorm = zero;
1467 if(fnorm != zero)
1468 {
1469 jj = 0;
1470 for( j=0; j<n; j++ )
1471 {
1472 l = ipvt[j];
1473 if(wa2[l] != zero)
1474 {
1475 sum = zero;
1476 ij = jj;
1477 for( i=0; i<=j; i++ )
1478 {
1479 sum += fjac[ij]*(qtf[i]/fnorm);
1480 ij += 1; /* fjac[i+m*j] */
1481 }
1482 gnorm = dmax1(gnorm,std::fabs(sum/wa2[l]));
1483 }
1484 jj += m;
1485 }
1486 }
1487
1488/*
1489* test for convergence of the gradient norm.
1490*/
1491 if(gnorm <= gtol)
1492 *info = 4;
1493 if( *info != 0)
1494 goto L300;
1495/*
1496* rescale if necessary.
1497*/
1498 if(mode != 2)
1499 {
1500 for( j=0; j<n; j++ )
1501 diag[j] = dmax1(diag[j],wa2[j]);
1502 }
1503
1504/*
1505* beginning of the inner loop.
1506*/
1507L200:
1508/*
1509* determine the levenberg-marquardt parameter.
1510*/
1511lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,&par,wa1,wa2,wa3,wa4);
1512/*
1513* store the direction p and x + p. calculate the norm of p.
1514*/
1515for( j=0; j<n; j++ )
1516 {
1517 wa1[j] = -wa1[j];
1518 wa2[j] = x[j] + wa1[j];
1519 wa3[j] = diag[j]*wa1[j];
1520 }
1521pnorm = enorm(n,wa3);
1522/*
1523* on the first iteration, adjust the initial step bound.
1524*/
1525if(iter == 1)
1526 delta = dmin1(delta,pnorm);
1527/*
1528* evaluate the function at x + p and calculate its norm.
1529*/
1530iflag = 1;
1531fcn(m,n,wa2,wa4,&iflag);
1532*nfev += 1;
1533if(iflag < 0)
1534 goto L300;
1535fnorm1 = enorm(m,wa4);
1536/*
1537* compute the scaled actual reduction.
1538*/
1539actred = -one;
1540if( (p1*fnorm1) < fnorm)
1541 {
1542 temp = fnorm1/fnorm;
1543 actred = one - temp * temp;
1544 }
1545/*
1546* compute the scaled predicted reduction and
1547* the scaled directional derivative.
1548*/
1549jj = 0;
1550for( j=0; j<n; j++ )
1551 {
1552 wa3[j] = zero;
1553 l = ipvt[j];
1554 temp = wa1[l];
1555 ij = jj;
1556 for( i=0; i<=j; i++ )
1557 {
1558 wa3[i] += fjac[ij]*temp;
1559 ij += 1; /* fjac[i+m*j] */
1560 }
1561 jj += m;
1562 }
1563temp1 = enorm(n,wa3)/fnorm;
1564temp2 = (std::sqrt(par)*pnorm)/fnorm;
1565prered = temp1*temp1 + (temp2*temp2)/p5;
1566dirder = -(temp1*temp1 + temp2*temp2);
1567/*
1568* compute the ratio of the actual to the predicted
1569* reduction.
1570*/
1571ratio = zero;
1572if(prered != zero)
1573 ratio = actred/prered;
1574/*
1575* update the step bound.
1576*/
1577if(ratio <= p25)
1578 {
1579 if(actred >= zero)
1580 temp = p5;
1581 else
1582 temp = p5*dirder/(dirder + p5*actred);
1583 if( ((p1*fnorm1) >= fnorm)
1584 || (temp < p1) )
1585 temp = p1;
1586 delta = temp*dmin1(delta,pnorm/p1);
1587 par = par/temp;
1588 }
1589else
1590 {
1591 if( (par == zero) || (ratio >= p75) )
1592 {
1593 delta = pnorm/p5;
1594 par = p5*par;
1595 }
1596 }
1597/*
1598* test for successful iteration.
1599*/
1600if(ratio >= p0001)
1601 {
1602/*
1603* successful iteration. update x, fvec, and their norms.
1604*/
1605 for( j=0; j<n; j++ )
1606 {
1607 x[j] = wa2[j];
1608 wa2[j] = diag[j]*x[j];
1609 }
1610 for( i=0; i<m; i++ )
1611 fvec[i] = wa4[i];
1612 xnorm = enorm(n,wa2);
1613 fnorm = fnorm1;
1614 iter += 1;
1615 }
1616/*
1617* tests for convergence.
1618*/
1619if( (std::fabs(actred) <= ftol)
1620 && (prered <= ftol)
1621 && (p5*ratio <= one) )
1622 *info = 1;
1623if(delta <= xtol*xnorm)
1624 *info = 2;
1625if( (std::fabs(actred) <= ftol)
1626 && (prered <= ftol)
1627 && (p5*ratio <= one)
1628 && ( *info == 2) )
1629 *info = 3;
1630if( *info != 0)
1631 goto L300;
1632/*
1633* tests for termination and stringent tolerances.
1634*/
1635if( *nfev >= maxfev)
1636 *info = 5;
1637if( (std::fabs(actred) <= MACHEP)
1638 && (prered <= MACHEP)
1639 && (p5*ratio <= one) )
1640 *info = 6;
1641if(delta <= MACHEP*xnorm)
1642 *info = 7;
1643if(gnorm <= MACHEP)
1644 *info = 8;
1645if( *info != 0)
1646 goto L300;
1647/*
1648* end of the inner loop. repeat if iteration unsuccessful.
1649*/
1650if(ratio < p0001)
1651 goto L200;
1652/*
1653* end of the outer loop.
1654*/
1655goto L30;
1656
1657L300:
1658/*
1659* termination, either normal or user imposed.
1660*/
1661if(iflag < 0)
1662 *info = iflag;
1663iflag = 0;
1664if(nprint > 0)
1665 fcn(m,n,x,fvec,&iflag);
1666/*
1667 last card of subroutine lmdif.
1668*/
1669}
1670}
1671}
1672
1673/************************fdjac2.c*************************/
1674
ext::function< Real(Real)> b
QL_REAL Real
real number
Definition: types.hpp:50
wrapper for MINPACK minimization routine
double MACHEP
Definition: lmdif.cpp:72
void lmdif(int m, int n, Real *x, Real *fvec, Real ftol, Real xtol, Real gtol, int maxfev, Real epsfcn, Real *diag, int mode, Real factor, int nprint, int *info, int *nfev, Real *fjac, int ldfjac, int *ipvt, Real *qtf, Real *wa1, Real *wa2, Real *wa3, Real *wa4, const QuantLib::MINPACK::LmdifCostFunction &fcn, const QuantLib::MINPACK::LmdifCostFunction &jacFcn)
Definition: lmdif.cpp:1131
int min0(int a, int b)
Definition: lmdif.cpp:229
ext::function< void(int, int, Real *, Real *, int *)> LmdifCostFunction
Definition: lmdif.hpp:38
void qrfac(int m, int n, Real *a, int, int pivot, int *ipvt, int, Real *rdiag, Real *acnorm, Real *wa)
Definition: lmdif.cpp:374
void fdjac2(int m, int n, Real *x, const Real *fvec, Real *fjac, int, int *iflag, Real epsfcn, Real *wa, const QuantLib::MINPACK::LmdifCostFunction &fcn)
Definition: lmdif.cpp:256
Real dmin1(Real a, Real b)
Definition: lmdif.cpp:221
int mod(int k, int m)
Definition: lmdif.cpp:238
Real enorm(int n, Real *x)
Definition: lmdif.cpp:80
double DWARF
Definition: lmdif.cpp:74
void qrsolv(int n, Real *r, int ldr, const int *ipvt, const Real *diag, const Real *qtb, Real *x, Real *sdiag, Real *wa)
Definition: lmdif.cpp:580
void lmpar(int n, Real *r, int ldr, int *ipvt, const Real *diag, Real *qtb, Real delta, Real *par, Real *x, Real *sdiag, Real *wa1, Real *wa2)
Definition: lmdif.cpp:814
Real dmax1(Real a, Real b)
Definition: lmdif.cpp:213
Definition: any.hpp:35
ext::shared_ptr< YieldTermStructure > r