122Real agiant,floatn,s1,s2,s3,xabs,x1max,x3max;
124static double rdwarf = 3.834e-20;
125static double rgiant = 1.304e19;
126static double zero = 0.0;
127static double one = 1.0;
135agiant = rgiant/floatn;
139xabs = std::fabs(x[i]);
140if( (xabs > rdwarf) && (xabs < agiant) )
157 s1 = one + s1*temp*temp;
173 s3 = one + s3*temp*temp;
190 temp = s1 + (s2/x1max)/x1max;
191 ans = x1max*std::sqrt(temp);
197 temp = s2*(one+(x3max/s2)*(x3max*s3));
199 temp = x3max*((s2/x3max)+(x3max*s3));
200 ans = std::sqrt(temp);
204 ans = x3max*std::sqrt(s3);
345 static double zero = 0.0;
349 eps = std::sqrt(temp);
351 for (j = 0; j <
n; j++) {
353 h = eps * std::fabs(temp);
357 fcn(m,
n, x, wa, iflag);
361 for (i = 0; i < m; i++) {
362 fjac[ij] = (wa[i] - fvec[i]) / h;
454int i,ij,jj,j,jp1,k,kmax,minmn;
456static double zero = 0.0;
457static double one = 1.0;
458static double p05 = 0.05;
466 acnorm[j] =
enorm(m,&a[ij]);
467 rdiag[j] = acnorm[j];
477for( j=0; j<minmn; j++ )
487 if(rdiag[k] > rdiag[kmax])
503rdiag[kmax] = rdiag[j];
515ajnorm =
enorm(m-j,&a[jj]);
534for( k=jp1; k<
n; k++ )
554 if( (pivot != 0) && (rdiag[k] != zero) )
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)
562 rdiag[k] =
enorm(m-j-1,&a[jp1+m*k]);
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;
679 for (j = 0; j <
n; j++) {
682 for (i = j; i <
n; i++) {
721 if(std::fabs(
r[kk]) < std::fabs(sdiag[k]))
723 cotan =
r[kk]/sdiag[k];
724 sin = p5/std::sqrt(p25+p25*cotan*cotan);
729 tan = sdiag[k]/
r[kk];
730 cos = p5/std::sqrt(p25+p25*tan*tan);
737 r[kk] = cos*
r[kk] + sin*sdiag[k];
738 temp = cos*wa[k] + sin*qtbpj;
739 qtbpj = -sin*wa[k] + cos*qtbpj;
748 for( i=kp1; i<
n; i++ )
750 temp = cos*
r[ik] + sin*sdiag[i];
751 sdiag[i] = -sin*
r[ik] + cos*sdiag[i];
773 if( (sdiag[j] == zero) && (nsing ==
n) )
781for( k=0; k<nsing; k++ )
789 for( i=jp1; i<nsing; i++ )
795 wa[j] = (wa[j] - sum)/sdiag[j];
921 int i, iter, ij, jj, j, jm1, jp1, k, l, nsing;
922 Real dxnorm, fp, gnorm, parc, parl, paru;
924 static double zero = 0.0;
925 static double p1 = 0.1;
926 static double p001 = 0.001;
935 for (j = 0; j <
n; j++) {
937 if ((
r[jj] == zero) && (nsing ==
n))
945 for( k=0; k<nsing; k++ )
948 wa1[j] = wa1[j]/
r[j+ldr*j];
954 for( i=0; i<=jm1; i++ )
956 wa1[i] -=
r[ij]*temp;
975 wa2[j] = diag[j]*x[j];
993 wa1[j] = diag[l]*(wa2[l]/dxnorm);
1003 for( i=0; i<=jm1; i++ )
1005 sum +=
r[ij]*wa1[i];
1009 wa1[j] = (wa1[j] - sum)/
r[j+ldr*j];
1013 parl = ((fp/delta)/temp)/temp;
1023 for( i=0; i<=j; i++ )
1025 sum +=
r[ij]*qtb[i];
1029 wa1[j] = sum/diag[l];
1040*par =
dmax1( *par,parl);
1041*par =
dmin1( *par,paru);
1043 *par = gnorm/dxnorm;
1054temp = std::sqrt( *par );
1056 wa1[j] = temp*diag[j];
1057qrsolv(
n,
r,ldr,ipvt,wa1,qtb,x,sdiag,wa2);
1059 wa2[j] = diag[j]*x[j];
1068if( (std::fabs(fp) <= p1*delta)
1069 || ((parl == zero) && (fp <= temp) && (temp < zero))
1078 wa1[j] = diag[l]*(wa2[l]/dxnorm);
1083 wa1[j] = wa1[j]/sdiag[j];
1089 for( i=jp1; i<
n; i++ )
1091 wa1[i] -=
r[ij]*temp;
1098parc = ((fp/delta)/temp)/temp;
1103 parl =
dmax1(parl, *par);
1105 paru =
dmin1(paru, *par);
1109*par =
dmax1(parl, *par + parc);
1134 int nprint,
int* info,
int* nfev,
Real* fjac,
1135 int ldfjac,
int* ipvt,
Real* qtf,
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;
1340if( (
n <= 0) || (m <
n) || (ldfjac < m) || (ftol < zero)
1341 || (xtol < zero) || (gtol < zero) || (maxfev <= 0)
1342 || (factor <= zero) )
1347 for( j=0; j<
n; j++ )
1349 if( diag[j] <= 0.0 )
1358fcn(m,
n,x,fvec,&iflag);
1362fnorm =
enorm(m,fvec);
1379 fdjac2(m,
n,x,fvec,fjac,ldfjac,&iflag,epsfcn,wa4, fcn);
1381 jacFcn(m,
n,x,fjac,&iflag);
1391 if(
mod(iter-1,nprint) == 0)
1393 fcn(m,
n,x,fvec,&iflag);
1401qrfac(m,
n,fjac,ldfjac,1,ipvt,
n,wa1,wa2,wa3);
1410 for( j=0; j<
n; j++ )
1413 if( wa2[j] == zero )
1422 for( j=0; j<
n; j++ )
1423 wa3[j] = diag[j] * x[j];
1426 delta = factor*xnorm;
1445 for( i=j; i<m; i++ )
1447 sum += fjac[ij] * wa4[i];
1450 temp = -sum / temp3;
1452 for( i=j; i<m; i++ )
1454 wa4[i] += fjac[ij] * temp;
1470 for( j=0; j<
n; j++ )
1477 for( i=0; i<=j; i++ )
1479 sum += fjac[ij]*(qtf[i]/fnorm);
1482 gnorm =
dmax1(gnorm,std::fabs(sum/wa2[l]));
1500 for( j=0; j<
n; j++ )
1501 diag[j] =
dmax1(diag[j],wa2[j]);
1511lmpar(
n,fjac,ldfjac,ipvt,diag,qtf,delta,&par,wa1,wa2,wa3,wa4);
1518 wa2[j] = x[j] + wa1[j];
1519 wa3[j] = diag[j]*wa1[j];
1526 delta =
dmin1(delta,pnorm);
1531fcn(m,
n,wa2,wa4,&iflag);
1535fnorm1 =
enorm(m,wa4);
1540if( (p1*fnorm1) < fnorm)
1542 temp = fnorm1/fnorm;
1543 actred = one - temp * temp;
1556 for( i=0; i<=j; i++ )
1558 wa3[i] += fjac[ij]*temp;
1563temp1 =
enorm(
n,wa3)/fnorm;
1564temp2 = (std::sqrt(par)*pnorm)/fnorm;
1565prered = temp1*temp1 + (temp2*temp2)/p5;
1566dirder = -(temp1*temp1 + temp2*temp2);
1573 ratio = actred/prered;
1582 temp = p5*dirder/(dirder + p5*actred);
1583 if( ((p1*fnorm1) >= fnorm)
1586 delta = temp*
dmin1(delta,pnorm/p1);
1591 if( (par == zero) || (ratio >= p75) )
1605 for( j=0; j<
n; j++ )
1608 wa2[j] = diag[j]*x[j];
1610 for( i=0; i<m; i++ )
1619if( (std::fabs(actred) <= ftol)
1621 && (p5*ratio <= one) )
1623if(delta <= xtol*xnorm)
1625if( (std::fabs(actred) <= ftol)
1627 && (p5*ratio <= one)
1637if( (std::fabs(actred) <=
MACHEP)
1639 && (p5*ratio <= one) )
1665 fcn(m,
n,x,fvec,&iflag);
ext::function< Real(Real)> b
wrapper for MINPACK minimization routine
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)
ext::function< void(int, int, Real *, Real *, int *)> LmdifCostFunction
void qrfac(int m, int n, Real *a, int, int pivot, int *ipvt, int, Real *rdiag, Real *acnorm, Real *wa)
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)
Real dmin1(Real a, Real b)
Real enorm(int n, Real *x)
void qrsolv(int n, Real *r, int ldr, const int *ipvt, const Real *diag, const Real *qtb, Real *x, Real *sdiag, Real *wa)
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)
Real dmax1(Real a, Real b)
ext::shared_ptr< YieldTermStructure > r