1	/************************lmdif*************************/
2
3	/*
4	The original Fortran version is Copyright (C) 1999 University of Chicago.
5	All rights reserved.
6
7	Redistribution and use in source and binary forms, with or
8	without modification, are permitted provided that the
9	following conditions are met:
10
11	1. Redistributions of source code must retain the above
12	copyright notice, this list of conditions and the following
13	disclaimer.
14
15	2. Redistributions in binary form must reproduce the above
16	copyright notice, this list of conditions and the following
17	disclaimer in the documentation and/or other materials
18	provided with the distribution.
19
20	3. The end-user documentation included with the
21	redistribution, if any, must include the following
22	acknowledgment:
23
24	"This product includes software developed by the
25	University of Chicago, as Operator of Argonne National
26	Laboratory.
27
28	Alternately, this acknowledgment may appear in the software
29	itself, if and wherever such third-party acknowledgments
30	normally appear.
31
32	4. WARRANTY DISCLAIMER. THE SOFTWARE IS SUPPLIED "AS IS"
33	WITHOUT WARRANTY OF ANY KIND. THE COPYRIGHT HOLDER, THE
34	UNITED STATES, THE UNITED STATES DEPARTMENT OF ENERGY, AND
35	THEIR EMPLOYEES: (1) DISCLAIM ANY WARRANTIES, EXPRESS OR
36	IMPLIED, INCLUDING BUT NOT LIMITED TO ANY IMPLIED WARRANTIES
37	OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE
38	OR NON-INFRINGEMENT, (2) DO NOT ASSUME ANY LEGAL LIABILITY
39	OR RESPONSIBILITY FOR THE ACCURACY, COMPLETENESS, OR
40	USEFULNESS OF THE SOFTWARE, (3) DO NOT REPRESENT THAT USE OF
41	THE SOFTWARE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS, (4)
42	DO NOT WARRANT THAT THE SOFTWARE WILL FUNCTION
43	UNINTERRUPTED, THAT IT IS ERROR-FREE OR THAT ANY ERRORS WILL
44	BE CORRECTED.
45
46	5. LIMITATION OF LIABILITY. IN NO EVENT WILL THE COPYRIGHT
47	HOLDER, THE UNITED STATES, THE UNITED STATES DEPARTMENT OF
48	ENERGY, OR THEIR EMPLOYEES: BE LIABLE FOR ANY INDIRECT,
49	INCIDENTAL, CONSEQUENTIAL, SPECIAL OR PUNITIVE DAMAGES OF
50	ANY KIND OR NATURE, INCLUDING BUT NOT LIMITED TO LOSS OF
51	PROFITS OR LOSS OF DATA, FOR ANY REASON WHATSOEVER, WHETHER
52	SUCH LIABILITY IS ASSERTED ON THE BASIS OF CONTRACT, TORT
53	(INCLUDING NEGLIGENCE OR STRICT LIABILITY), OR OTHERWISE,
54	EVEN IF ANY OF SAID PARTIES HAS BEEN WARNED OF THE
55	POSSIBILITY OF SUCH LOSS OR DAMAGES.
56
57
58	C translation Copyright (C) Steve Moshier
59
60	What you see here may be used freely but it comes with no support
61	or guarantee.
62	*/
63
64	#include <ql/math/optimization/lmdif.hpp>
65	#include <cmath>
66	#include <cstdio>
67
68	namespace QuantLib {
69	namespace MINPACK {
70	#define BUG 0
71	/* resolution of arithmetic */
72	double MACHEP = 1.2e-16;
73	/* smallest nonzero number */
74	double DWARF = 1.0e-38;
75
76
77
78
79
80	Real enorm(int n,Real* x)
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	*/
121	int i;
122	Real agiant,floatn,s1,s2,s3,xabs,x1max,x3max;
123	Real ans, temp;
124	static double rdwarf = 3.834e-20;
125	static double rgiant = 1.304e19;
126	static double zero = 0.0;
127	static double one = 1.0;
128
129	s1 = zero;
130	s2 = zero;
131	s3 = zero;
132	x1max = zero;
133	x3max = zero;
134	floatn = n;
135	agiant = rgiant/floatn;
136
137	for( i=0; i<n; i++ )
138	{
139	xabs = std::fabs(x: x[i]);
140	if( (xabs > rdwarf) && (xabs < agiant) )
141	{
142	/*
143	* sum for intermediate components.
144	*/
145	s2 += xabs*xabs;
146	continue;
147	}
148
149	if(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	*/
170	if(xabs > x3max)
171	{
172	temp = x3max/xabs;
173	s3 = one + s3*temp*temp;
174	x3max = xabs;
175	}
176	else
177	{
178	if(xabs != zero)
179	{
180	temp = xabs/x3max;
181	s3 += temp*temp;
182	}
183	}
184	}
185	/*
186	* calculation of norm.
187	*/
188	if(s1 != zero)
189	{
190	temp = s1 + (s2/x1max)/x1max;
191	ans = x1max*std::sqrt(x: temp);
192	return(ans);
193	}
194	if(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(x: temp);
201	}
202	else
203	{
204	ans = x3max*std::sqrt(x: s3);
205	}
206	return(ans);
207	/*
208	* last card of function enorm.
209	*/
210	}
211	/************************lmmisc.c*************************/
212
213	Real dmax1(Real a,Real b)
214	{
215	if( a >= b )
216	return(a);
217	else
218	return(b);
219	}
220
221	Real dmin1(Real a,Real b)
222	{
223	if( a <= b )
224	return(a);
225	else
226	return(b);
227	}
228
229	int min0(int a,int b)
230
231	{
232	if( a <= b )
233	return(a);
234	else
235	return(b);
236	}
237
238	int mod( int k, int m )
239	{
240	return( 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
256	void 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,
265	const QuantLib::MINPACK::LmdifCostFunction& fcn) {
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(a: epsfcn, b: MACHEP);
349	eps = std::sqrt(x: temp);
350	ij = 0;
351	for (j = 0; j < n; j++) {
352	temp = x[j];
353	h = eps * std::fabs(x: 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
373	void
374	qrfac(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	*/
454	int i,ij,jj,j,jp1,k,kmax,minmn;
455	Real ajnorm,sum,temp;
456	static double zero = 0.0;
457	static double one = 1.0;
458	static double p05 = 0.05;
459
460	/*
461	* compute the initial column norms and initialize several arrays.
462	*/
463	ij = 0;
464	for( j=0; j<n; j++ )
465	{
466	acnorm[j] = enorm(n: m,x: &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	*/
476	minmn = min0(a: m,b: n);
477	for( j=0; j<minmn; j++ )
478	{
479	if(pivot == 0)
480	goto L40;
481	/*
482	* bring the column of largest norm into the pivot position.
483	*/
484	kmax = j;
485	for( k=j; k<n; k++ )
486	{
487	if(rdiag[k] > rdiag[kmax])
488	kmax = k;
489	}
490	if(kmax == j)
491	goto L40;
492
493	ij = m * j;
494	jj = m * kmax;
495	for( 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	}
503	rdiag[kmax] = rdiag[j];
504	wa[kmax] = wa[j];
505	k = ipvt[j];
506	ipvt[j] = ipvt[kmax];
507	ipvt[kmax] = k;
508
509	L40:
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	*/
514	jj = j + m*j;
515	ajnorm = enorm(n: m-j,x: &a[jj]);
516	if(ajnorm == zero)
517	goto L100;
518	if(a[jj] < zero)
519	ajnorm = -ajnorm;
520	ij = jj;
521	for( i=j; i<m; i++ )
522	{
523	a[ij] /= ajnorm;
524	ij += 1; /* [i+m*j] */
525	}
526	a[jj] += one;
527	/*
528	* apply the transformation to the remaining columns
529	* and update the norms.
530	*/
531	jp1 = j + 1;
532	if(jp1 < n )
533	{
534	for( 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( a: zero, b: one-temp*temp );
558	rdiag[k] *= std::sqrt(x: temp);
559	temp = rdiag[k]/wa[k];
560	if( (p05*temp*temp) <= MACHEP)
561	{
562	rdiag[k] = enorm(n: m-j-1,x: &a[jp1+m*k]);
563	wa[k] = rdiag[k];
564	}
565	}
566	}
567	}
568
569	L100:
570	rdiag[j] = -ajnorm;
571	}
572	/*
573	* last card of subroutine qrfac.
574	*/
575	}
576
577	/************************qrsolv.c*************************/
578
579
580	void 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	*/
694	for( 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	*/
700	l = ipvt[j];
701	if(diag[l] == zero)
702	goto L90;
703	for( k=j; k<n; k++ )
704	sdiag[k] = zero;
705	sdiag[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	*/
711	qtbpj = zero;
712	for( 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(x: r[kk]) < std::fabs(x: sdiag[k]))
722	{
723	cotan = r[kk]/sdiag[k];
724	sin = p5/std::sqrt(x: p25+p25*cotan*cotan);
725	cos = sin*cotan;
726	}
727	else
728	{
729	tan = sdiag[k]/r[kk];
730	cos = p5/std::sqrt(x: 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	}
757	L90:
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	*/
770	nsing = n;
771	for( 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	}
778	if(nsing < 1)
779	goto L150;
780
781	for( 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	}
797	L150:
798	/*
799	* permute the components of z back to components of x.
800	*/
801	for( 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
814	void 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	}
943	if(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
963	for( 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	*/
973	iter = 0;
974	for( j=0; j<n; j++ )
975	wa2[j] = diag[j]*x[j];
976	dxnorm = enorm(n,x: wa2);
977	fp = dxnorm - delta;
978	if(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	*/
987	parl = zero;
988	if(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,x: wa1);
1013	parl = ((fp/delta)/temp)/temp;
1014	}
1015	/*
1016	* calculate an upper bound, paru, for the zero of the function.
1017	*/
1018	jj = 0;
1019	for( 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	}
1032	gnorm = enorm(n,x: wa1);
1033	paru = gnorm/delta;
1034	if(paru == zero)
1035	paru = DWARF/dmin1(a: delta,b: 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( a: *par,b: parl);
1041	*par = dmin1( a: *par,b: paru);
1042	if( *par == zero)
1043	*par = gnorm/dxnorm;
1044	/*
1045	* beginning of an iteration.
1046	*/
1047	L150:
1048	iter += 1;
1049	/*
1050	* evaluate the function at the current value of par.
1051	*/
1052	if( *par == zero)
1053	*par = dmax1(a: DWARF,b: p001*paru);
1054	temp = std::sqrt( x: *par );
1055	for( j=0; j<n; j++ )
1056	wa1[j] = temp*diag[j];
1057	qrsolv(n,r,ldr,ipvt,diag: wa1,qtb,x,sdiag,wa: wa2);
1058	for( j=0; j<n; j++ )
1059	wa2[j] = diag[j]*x[j];
1060	dxnorm = enorm(n,x: wa2);
1061	temp = fp;
1062	fp = 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	*/
1068	if( (std::fabs(x: fp) <= p1*delta)
1069	|| ((parl == zero) && (fp <= temp) && (temp < zero))
1070	|| (iter == 10) )
1071	goto L220;
1072	/*
1073	* compute the newton correction.
1074	*/
1075	for( j=0; j<n; j++ )
1076	{
1077	l = ipvt[j];
1078	wa1[j] = diag[l]*(wa2[l]/dxnorm);
1079	}
1080	jj = 0;
1081	for( 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	}
1097	temp = enorm(n,x: wa1);
1098	parc = ((fp/delta)/temp)/temp;
1099	/*
1100	* depending on the sign of the function, update parl or paru.
1101	*/
1102	if(fp > zero)
1103	parl = dmax1(a: parl, b: *par);
1104	if(fp < zero)
1105	paru = dmin1(a: paru, b: *par);
1106	/*
1107	* compute an improved estimate for par.
1108	*/
1109	*par = dmax1(a: parl, b: *par + parc);
1110	/*
1111	* end of an iteration.
1112	*/
1113	goto L150;
1114
1115	L220:
1116	/*
1117	* termination.
1118	*/
1119	if(iter == 0)
1120	*par = zero;
1121	/*
1122	* last card of subroutine lmpar.
1123	*/
1124	}
1125
1126	/*********************** lmdif.c ****************************/
1127
1128
1129
1130
1131	void 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,
1137	const QuantLib::MINPACK::LmdifCostFunction& fcn,
1138	const QuantLib::MINPACK::LmdifCostFunction& jacFcn)
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	*/
1322	int i,iflag,ij,jj,iter,j,l;
1323	Real actred,delta=0,dirder,fnorm,fnorm1,gnorm;
1324	Real par,pnorm,prered,ratio;
1325	Real sum,temp,temp1,temp2,temp3,xnorm=0;
1326	static double one = 1.0;
1327	static double p1 = 0.1;
1328	static double p5 = 0.5;
1329	static double p25 = 0.25;
1330	static double p75 = 0.75;
1331	static double p0001 = 1.0e-4;
1332	static double zero = 0.0;
1333
1334	*info = 0;
1335	iflag = 0;
1336	*nfev = 0;
1337	/*
1338	* check the input parameters for errors.
1339	*/
1340	if( (n <= 0) || (m < n) || (ldfjac < m) || (ftol < zero)
1341	|| (xtol < zero) || (gtol < zero) || (maxfev <= 0)
1342	|| (factor <= zero) )
1343	goto L300;
1344
1345	if( 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	*/
1357	iflag = 1;
1358	fcn(m,n,x,fvec,&iflag);
1359	*nfev = 1;
1360	if(iflag < 0)
1361	goto L300;
1362	fnorm = enorm(n: m,x: fvec);
1363	/*
1364	* initialize levenberg-marquardt parameter and iteration counter.
1365	*/
1366	par = zero;
1367	iter = 1;
1368	/*
1369	* beginning of the outer loop.
1370	*/
1371
1372	L30:
1373
1374	/*
1375	* calculate the jacobian matrix.
1376	*/
1377	iflag = 2;
1378	if (!jacFcn)
1379	fdjac2(m,n,x,fvec,fjac,ldfjac,iflag: &iflag,epsfcn,wa: wa4, fcn);
1380	else // use user supplied jacobian calculation
1381	jacFcn(m,n,x,fjac,&iflag);
1382	*nfev += n;
1383	if(iflag < 0)
1384	goto L300;
1385	/*
1386	* if requested, call fcn to enable printing of iterates.
1387	*/
1388	if( nprint > 0 )
1389	{
1390	iflag = 0;
1391	if(mod(k: iter-1,m: 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	*/
1401	qrfac(m,n,a: fjac,ldfjac,pivot: 1,ipvt,n,rdiag: wa1,acnorm: wa2,wa: 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	*/
1406	if(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,x: 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	*/
1435	for( i=0; i<m; i++ )
1436	wa4[i] = fvec[i];
1437	jj = 0;
1438	for( 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(a: gnorm,b: std::fabs(x: 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(a: diag[j],b: wa2[j]);
1502	}
1503
1504	/*
1505	* beginning of the inner loop.
1506	*/
1507	L200:
1508	/*
1509	* determine the levenberg-marquardt parameter.
1510	*/
1511	lmpar(n,r: fjac,ldr: ldfjac,ipvt,diag,qtb: qtf,delta,par: &par,x: wa1,sdiag: wa2,wa1: wa3,wa2: wa4);
1512	/*
1513	* store the direction p and x + p. calculate the norm of p.
1514	*/
1515	for( 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	}
1521	pnorm = enorm(n,x: wa3);
1522	/*
1523	* on the first iteration, adjust the initial step bound.
1524	*/
1525	if(iter == 1)
1526	delta = dmin1(a: delta,b: pnorm);
1527	/*
1528	* evaluate the function at x + p and calculate its norm.
1529	*/
1530	iflag = 1;
1531	fcn(m,n,wa2,wa4,&iflag);
1532	*nfev += 1;
1533	if(iflag < 0)
1534	goto L300;
1535	fnorm1 = enorm(n: m,x: wa4);
1536	/*
1537	* compute the scaled actual reduction.
1538	*/
1539	actred = -one;
1540	if( (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	*/
1549	jj = 0;
1550	for( 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	}
1563	temp1 = enorm(n,x: wa3)/fnorm;
1564	temp2 = (std::sqrt(x: par)*pnorm)/fnorm;
1565	prered = temp1*temp1 + (temp2*temp2)/p5;
1566	dirder = -(temp1*temp1 + temp2*temp2);
1567	/*
1568	* compute the ratio of the actual to the predicted
1569	* reduction.
1570	*/
1571	ratio = zero;
1572	if(prered != zero)
1573	ratio = actred/prered;
1574	/*
1575	* update the step bound.
1576	*/
1577	if(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(a: delta,b: pnorm/p1);
1587	par = par/temp;
1588	}
1589	else
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	*/
1600	if(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,x: wa2);
1613	fnorm = fnorm1;
1614	iter += 1;
1615	}
1616	/*
1617	* tests for convergence.
1618	*/
1619	if( (std::fabs(x: actred) <= ftol)
1620	&& (prered <= ftol)
1621	&& (p5*ratio <= one) )
1622	*info = 1;
1623	if(delta <= xtol*xnorm)
1624	*info = 2;
1625	if( (std::fabs(x: actred) <= ftol)
1626	&& (prered <= ftol)
1627	&& (p5*ratio <= one)
1628	&& ( *info == 2) )
1629	*info = 3;
1630	if( *info != 0)
1631	goto L300;
1632	/*
1633	* tests for termination and stringent tolerances.
1634	*/
1635	if( *nfev >= maxfev)
1636	*info = 5;
1637	if( (std::fabs(x: actred) <= MACHEP)
1638	&& (prered <= MACHEP)
1639	&& (p5*ratio <= one) )
1640	*info = 6;
1641	if(delta <= MACHEP*xnorm)
1642	*info = 7;
1643	if(gnorm <= MACHEP)
1644	*info = 8;
1645	if( *info != 0)
1646	goto L300;
1647	/*
1648	* end of the inner loop. repeat if iteration unsuccessful.
1649	*/
1650	if(ratio < p0001)
1651	goto L200;
1652	/*
1653	* end of the outer loop.
1654	*/
1655	goto L30;
1656
1657	L300:
1658	/*
1659	* termination, either normal or user imposed.
1660	*/
1661	if(iflag < 0)
1662	*info = iflag;
1663	iflag = 0;
1664	if(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
1675