37 fprintf(stderr,
"Numerical Recipes run-time error...\n");
38 fprintf(stderr,
"%s\n", error_text);
39 fprintf(stderr,
"...now exiting to system...\n");
42 #define NRSIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a)) 60 static long ix1, ix2, ix3;
66 if (*idum < 0 || iff == 0)
69 ix1 = (
IC1 - (*idum)) %
M1;
74 for (j = 1;j <= 97;j++)
85 j = 1 + ((97 * ix3) /
M3);
87 nrerror(
"RAN1: This cannot happen.");
110 double fac, r,
v1, v2;
116 v1 = 2.0 *
ran1(idum) - 1.0;
117 v2 = 2.0 *
ran1(idum) - 1.0;
118 r = v1 * v1 + v2 * v2;
140 double fac, r,
v1, v2;
146 v1 = 2.0 *
ran1(idum) - 1.0;
147 v2 = 2.0 *
ran1(idum) - 1.0;
148 r = v1 * v1 + v2 * v2;
151 fac =
sqrt(nu*(pow(r,-2.0/nu) -1.0)/r);
165 void ksone(
double data[],
int n,
double(*func)(
double),
double *
d,
double * prob)
168 double fn, ff, en, dt, fo=0.;
171 for (
int j=1;
j<=
n;
j++)
174 ff = (*func)(data[
j]);
175 dt =
XMIPP_MAX(fabs(fo - ff), fabs(fn - ff));
187 double a2, fac=2.0, sum=0.0, term, termbf=0.0;
188 double EPS1=0.001, EPS2=1.0e-8;
190 a2 = -2.0 * alam * alam;
191 for (j = 1; j<= 100; j++)
193 term = fac * exp(a2*j*j);
195 if (fabs(term) <= EPS1*termbf || fabs(term) <= EPS2*sum)
206 void indexx(
int n,
double arrin[],
int indx[])
208 int l,
j,
ir, indxt,
i;
211 for (j = 1;j <=
n;j++)
218 q = arrin[(indxt=indx[--l])];
221 q = arrin[(indxt=indx[
ir])];
233 if (j < ir && arrin[indx[j]] < arrin[indx[j+1]])
235 if (q < arrin[indx[j]])
253 double xx,
y, ans, ans1, ans2;
255 if ((ax = fabs(x)) < 8.0)
258 ans1 = 57568490574.0 + y * (-13362590354.0 +
260 + y * (-11214424.18 +
262 y * (-184.9052456)))));
263 ans2 = 57568490411.0 + y * (1029532985.0 +
274 xx = ax - 0.785398164;
275 ans1 = 1.0 + y * (-0.1098628627e-2 + y * (0.2734510407e-4
276 + y * (-0.2073370639e-5 + y * 0.2093887211e-6)));
277 ans2 = -0.1562499995e-1 + y * (0.1430488765e-3
278 + y * (-0.6911147651e-5 + y * (0.7621095161e-6
279 - y * 0.934935152e-7)));
280 ans =
sqrt(0.636619772 / ax) * (cos(xx) * ans1 - z * sin(xx) * ans2);
289 if ((ax = fabs(x)) < 3.75)
293 ans = 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
294 + y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
299 ans = (exp(ax) /
sqrt(ax)) * (0.39894228 + y * (0.1328592e-1
300 + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
301 + y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1
302 + y * 0.392377e-2))))))));
312 if ((ax = fabs(x)) < 3.75)
316 ans = ax * (0.5 + y * (0.87890594 + y * (0.51498869 + y * (0.15084934
317 + y * (0.2658733e-1 + y * (0.301532e-2 + y * 0.32411e-3))))));
322 ans = 0.2282967e-1 + y * (-0.2895312e-1 + y * (0.1787654e-1
324 ans = 0.39894228 + y * (-0.3988024e-1 + y * (-0.362018e-2
325 + y * (0.163801e-2 + y * (-0.1031555e-1 + y * ans))));
326 ans *= (exp(ax) /
sqrt(ax));
328 return x < 0.0 ? -ans : ans;
334 double d = 0.0, dd = 0.0, sv,
y, y2;
337 if ((x - a)*(x - b) > 0.0)
338 nrerror(
"x not in range in routine chebev");
339 y2 = 2.0 * (y = (2.0 * x - a -
b) / (b - a));
340 for (j = m - 1;j >= 1;j--)
343 d = y2 * d - dd + c[
j];
346 return y*d - dd + 0.5*c[0];
351 void beschb(
double x,
double *gam1,
double *gam2,
double *gampl,
double *gammi)
356 -1.142022680371172e0, 6.516511267076e-3,
357 3.08709017308e-4, -3.470626964e-6, 6.943764e-9,
358 3.6780e-11, -1.36e-13
362 1.843740587300906e0, -0.076852840844786e0,
363 1.271927136655e-3, -4.971736704e-6, -3.3126120e-8,
364 2.42310e-10, -1.70e-13, -1.0e-15
367 xx = 8.0 * x * x - 1.0;
370 *gampl = *gam2 - x * (*gam1);
371 *gammi = *gam2 + x * (*gam1);
378 #define FPMIN 1.0e-30 381 void bessjy(
double x,
double xnu,
double *rj,
double *ry,
double *rjp,
double *ryp)
384 double a,
b, br, bi,
c, cr, ci,
d,
del, del1, den,
di, dlr, dli, dr, e,
f,
fact, fact2,
385 fact3, ff, gam, gam1, gam2, gammi, gampl, h, p, pimu, pimu2, q, r, rjl,
386 rjl1, rjmu, rjp1, rjpl, rjtemp, ry1, rymu, rymup, rytemp, sum, sum1,
387 temp,
w, x2,
xi, xi2, xmu, xmu2;
389 if (x <= 0.0 || xnu < 0.0)
390 nrerror(
"bad arguments in bessjy");
391 nl = (x <
XMIN ? (int)(xnu + 0.5) :
XMIPP_MAX(0, (
int)(xnu - x + 1.5)));
404 for (i = 1;i <=
MAXIT;i++)
418 if (fabs(del - 1.0) <
EPS)
422 nrerror(
"x too large in bessjy; try asymptotic expansion");
428 for (l = nl;l >= 1;l--)
430 rjtemp = fact * rjl + rjpl;
432 rjpl = fact * rjtemp - rjl;
442 fact = (fabs(pimu) <
EPS ? 1.0 : pimu / sin(pimu));
445 fact2 = (fabs(e) <
EPS ? 1.0 : sinh(e) / e);
446 beschb(xmu, &gam1, &gam2, &gampl, &gammi);
447 ff = 2.0 /
PI * fact * (gam1 * cosh(e) + gam2 * fact2 *
d);
449 p = e / (gampl *
PI);
450 q = 1.0 / (e *
PI * gammi);
452 fact3 = (fabs(pimu2) <
EPS ? 1.0 : sin(pimu2) / pimu2);
453 r =
PI * pimu2 * fact3 * fact3;
458 for (i = 1;i <=
MAXIT;i++)
460 ff = (i * ff + p + q) / (i * i - xmu2);
464 del = c * (ff + r * q);
466 del1 = c * p - i *
del;
468 if (fabs(del) < (1.0 + fabs(sum))*
EPS)
472 nrerror(
"bessy series failed to converge");
475 rymup = xmu * xi * rymu - ry1;
476 rjmu = w / (rymup - f * rymu);
485 fact = a * xi / (p * p + q * q);
488 den = br * br + bi * bi;
491 dlr = cr * dr - ci *
di;
492 dli = cr * di + ci * dr;
493 temp = p * dlr - q * dli;
494 q = p * dli + q * dlr;
496 for (i = 2;i <=
MAXIT;i++)
502 if (fabs(dr) + fabs(di) <
FPMIN)
504 fact = a / (cr * cr + ci * ci);
507 if (fabs(cr) + fabs(ci) <
FPMIN)
509 den = dr * dr + di *
di;
512 dlr = cr * dr - ci *
di;
513 dli = cr * di + ci * dr;
514 temp = p * dlr - q * dli;
515 q = p * dli + q * dlr;
517 if (fabs(dlr - 1.0) + fabs(dli) <
EPS)
521 nrerror(
"cf2 failed in bessjy");
523 rjmu =
sqrt(w / ((p - f) * gam + q));
526 rymup = rymu * (p + q / gam);
527 ry1 = xmu * xi * rymu - rymup;
532 for (i = 1;i <= nl;i++)
534 rytemp = (xmu +
i) * xi2 * ry1 - rymu;
539 *ryp = xnu * xi * rymu - ry1;
549 return (x == 0) ? 0 :
sqrt(2 / (
PI*x))*sinh(x);
553 return (x == 0) ? 0 :
sqrt(2 / (
PI*x))*(cosh(x) - sinh(x) /
x);
557 return (x == 0) ? 0 :
bessi0(x) - ((2*1) / x) *
bessi1(x);
565 return (x == 0) ? 0 :
bessi1(x) - ((2*2) / x) *
bessi2(x);
573 return (x == 0) ? 0 :
bessi2(x) - ((2*3) / x) *
bessi3(x);
577 double rj, ry, rjp, ryp;
578 bessjy(x, 1.5, &rj, &ry, &rjp, &ryp);
583 double rj, ry, rjp, ryp;
584 bessjy(x, 3.5, &rj, &ry, &rjp, &ryp);
592 static double cof[6] =
594 76.18009173, -86.50532033, 24.01409822,
595 -1.231739516, 0.120858003e-2, -0.536382e-5
601 tmp -= (x + 0.5) *
log(tmp);
603 for (j = 0;j <= 5;j++)
608 return -tmp +
log(2.50662827465*ser);
615 if (x < 0.0 || x > 1.0)
616 nrerror(
"Bad x in routine BETAI");
617 if (x == 0.0 || x == 1.0)
621 if (x < (a + 1.0) / (a + b + 2.0))
624 return 1.0 -bt*
betacf(b, a, 1.0 - x) /
b;
632 double qap, qam, qab, em, tem,
d;
633 double bz, bm = 1.0, bp, bpp;
634 double az = 1.0, am = 1.0, ap, app, aold;
640 bz = 1.0 - qab * x / qap;
641 for (m = 1;m <=
ITMAX;m++)
645 d = em * (b - em) * x / ((qam + tem) * (a + tem));
648 d = -(a + em) * (qab + em) * x / ((qap + tem) * (a + tem));
656 if (fabs(az - aold) < (
EPS*fabs(az)))
659 nrerror(
"a or b too big, or ITMAX too small in BETACF");
668 #define GOLD 1.618034 671 #define MAX(a,b) ((a) > (b) ? (a) : (b)) 672 #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) 673 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); 674 #define F1DIM(x,f) {\ 675 for (int j = 1; j<=ncom; j++) \ 676 xt[j] = pcom[j] + x * xicom[j]; \ 677 f = (*func)(xt,prm);} 679 void mnbrak(
double *ax,
double *bx,
double *cx,
680 double *fa,
double *fb,
double *fc,
double(*func)(
double *,
void*),
681 void *
prm,
int ncom,
double *pcom,
double *xicom)
683 double ulim,
u, r, q, fu, dum;
684 std::vector<double>
buffer(ncom);
685 auto *xt= buffer.data()-1;
691 SHFT(dum, *ax, *bx, dum)
692 SHFT(dum, *fb, *fa, dum)
694 *cx = (*bx) +
GOLD * (*bx - *ax);
698 r = (*bx - *ax) * (*fb - *fc);
699 q = (*bx - *cx) * (*fb - *fa);
700 u = (*bx) - ((*bx - *cx) * q - (*bx - *ax) * r) /
702 ulim = (*bx) +
GLIMIT * (*cx - *bx);
703 if ((*bx - u)*(u - *cx) > 0.0)
720 u = (*cx) +
GOLD * (*cx - *bx);
723 else if ((*cx - u)*(u - ulim) > 0.0)
728 SHFT(*bx, *cx, u, *cx +
GOLD*(*cx - *bx))
731 SHFT(*fb, *fc, fu, aux)
734 else if ((u - ulim)*(ulim - *cx) >= 0.0)
741 u = (*cx) +
GOLD * (*cx - *bx);
744 SHFT(*ax, *bx, *cx, u)
745 SHFT(*fa, *fb, *fc, fu)
755 #define CGOLD 0.3819660 757 double brent(
double ax,
double bx,
double cx,
double(*func)(
double *,
void*),
758 void *
prm,
double tol,
double *xmin,
759 int ncom,
double *pcom,
double *xicom)
762 double a,
b,
d, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2,
u, v,
w,
x, xm;
764 std::vector<double>
buffer(ncom);
765 auto *xt= buffer.data()-1;
767 a = (ax < cx ? ax : cx);
768 b = (ax > cx ? ax : cx);
772 for (iter = 1;iter <=
ITMAX;iter++)
775 tol2 = 2.0 * (tol1 = tol * fabs(x) +
ZEPS);
776 if (fabs(x - xm) <= (tol2 - 0.5*(b -
a)))
783 r = (x -
w) * (fx - fv);
784 q = (x - v) * (fx - fw);
785 p = (x - v) * q - (x - w) * r;
792 if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a - x) || p >= q*(b - x))
793 d =
CGOLD * (e = (x >= xm ? a -
x : b -
x));
798 if (u - a < tol2 || b - u < tol2)
799 d =
SIGN(tol1, xm - x);
804 d =
CGOLD * (e = (x >= xm ? a -
x : b -
x));
806 u = (fabs(d) >= tol1 ? x +
d : x +
SIGN(tol1, d));
823 if (fu <= fw || w == x)
830 else if (fu <= fv || v == x || v == w)
837 nrerror(
"Too many iterations in brent");
848 void linmin(
double *p,
double *
xi,
int n,
double &fret,
849 double(*func)(
double *,
void*),
void *
prm)
852 double xx, xmin, fx, fb, fa, bx, ax;
855 std::vector<double>
buffer(2*n);
856 auto *pcom= buffer.data()-1;
857 auto *xicom= pcom +
n;
858 for (j = 1;j <=
n;j++)
866 mnbrak(&ax, &xx, &bx, &fa, &fx, &fb, func, prm, ncom, pcom, xicom);
867 fret =
brent(ax, xx, bx, func, prm,
TOL, &xmin, ncom, pcom, xicom);
868 for (j = 1;j <=
n;j++)
877 void powell(
double *p,
double *xi,
int n,
double ftol,
int &
iter,
878 double &fret,
double(*func)(
double *,
void *),
void *
prm,
882 double t, fptt, fp,
del;
883 std::vector<double>
buffer(3*n);
884 auto *pt= buffer.data()-1;
887 bool different_from_0;
889 fret = (*func)(p,
prm);
890 for (j = 1;j <=
n;j++)
893 for (iter = 1;;(
iter)++)
898 std::cout << iter <<
" (" << p[1];
899 for (
int co = 2; co <=
n; co++)
900 std::cout <<
"," << p[co];
901 std::cout <<
")--->" << fret << std::endl;
908 for (i = 1;i <=
n;i++)
910 different_from_0 =
false;
911 for (j = 1;j <=
n;j++)
915 different_from_0 =
true;
917 if (different_from_0)
920 linmin(p, xit, n, fret, func, prm);
921 if (fabs(fptt - fret) > del)
923 del = fabs(fptt - fret);
933 for (
int co = 2; co <=
n; co++)
940 std::cout <<
")--->" << fret << std::endl;
945 if (2.0*fabs(fp - fret) <= ftol*(fabs(fp) + fabs(fret)) || n==1)
return;
947 nrerror(
"Too many iterations in routine POWELL");
948 for (j = 1;j <=
n;j++)
950 ptt[
j] = 2.0 * p[
j] - pt[
j];
951 xit[
j] = p[
j] - pt[
j];
954 fptt = (*func)(ptt,
prm);
957 #define SQR(a) ((a)*(a)) 958 t = 2.0 * (fp - 2.0 * fret + fptt) *
SQR(fp - fret - del) - del *
SQR(fp - fptt);
961 linmin(p, xit, n, fret, func, prm);
962 for (j = 1;j <=
n;j++)
963 xi[j*n+ibig] = xit[j];
1011 void _nnls_g1(
double a,
double b,
double *cterm,
double *sterm,
double *sig)
1015 if (fabs(a) > fabs(b))
1019 yr =
sqrt(d1 * d1 + 1.);
1021 *cterm = (a >= 0.0 ? fabs(d1) : -fabs(d1));
1022 *sterm = (*cterm) * xr;
1023 *sig = fabs(a) * yr;
1029 yr =
sqrt(d1 * d1 + 1.);
1031 *sterm = (b >= 0.0 ? fabs(d1) : -fabs(d1));
1032 *cterm = (*sterm) * xr;
1033 *sig = fabs(b) * yr;
1059 double *
u,
int u_dim1,
double *up,
1076 double d1, d2,
b, clinv,
cl, sm;
1077 int incr,
k,
j, i2, i3, i4;
1080 if (mode != 1 && mode != 2)
1082 if (m < 1 || u == NULL || u_dim1 < 1 || cm == NULL)
1084 if (lpivot < 0 || lpivot >= l1 || l1 >= m)
1087 cl = (d1 = u[lpivot*u_dim1], fabs(d1));
1095 for (j = l1; j <
m; j++)
1097 d2 = (d1 = u[j*u_dim1], fabs(d1));
1105 d1 = u[lpivot*u_dim1] * clinv;
1107 for (j = l1; j <
m; j++)
1109 d1 = u[j*u_dim1] * clinv;
1113 if (u[lpivot*u_dim1] > 0.)
1115 *up = u[lpivot*u_dim1] -
cl;
1116 u[lpivot*u_dim1] =
cl;
1120 b = (*up) * u[lpivot*u_dim1];
1125 i2 = 1 - icv + ice * lpivot;
1126 incr = ice * (l1 - lpivot);
1127 for (j = 0; j < ncv; j++)
1132 sm = cm[i2-1] * (*up);
1133 for (k = l1; k <
m; k++)
1135 sm += cm[i3-1] * u[k*u_dim1];
1141 cm[i2-1] += sm * (*up);
1142 for (k = l1; k <
m; k++)
1144 cm[i4-1] += sm * u[k*u_dim1];
1166 double *
a,
int m,
int n,
1184 int pfeas, ret = 0, iz, jz, iz1, iz2, npp1, *
index;
1185 double d1, d2, sm, up, ss, *
w, *zz;
1186 int iter,
k,
j = 0, l, itmax, izmax = 0, nsetp, ii, jj = 0, ip;
1187 double temp, wmax, t, alpha, asave,
dummy, unorm, ztest, cc;
1191 if (m <= 0 || n <= 0 || a == NULL || b == NULL || x == NULL)
1197 w = (
double*)calloc(n,
sizeof(
double));
1201 zz = (
double*)calloc(m,
sizeof(
double));
1205 index = (
int*)calloc(n,
sizeof(
int));
1206 if (w == NULL || zz == NULL || index == NULL)
1210 for (k = 0; k <
n; k++)
1224 while (iz1 <= iz2 && nsetp < m)
1227 for (iz = iz1; iz <= iz2; iz++)
1231 for (l = npp1; l <
m; l++)
1232 sm += a[j*m+l] * b[l];
1239 for (wmax = 0., iz = iz1; iz <= iz2; iz++)
1259 asave = a[j*m+npp1];
1260 _nnls_h12(1, npp1, npp1 + 1, m, &a[j*m+0], 1, &up, &dummy, 1, 1, 0);
1263 for (l = 0; l < nsetp; l++)
1268 unorm =
sqrt(unorm);
1269 d2 = unorm + (d1 = a[j*m+npp1], fabs(d1)) * 0.01;
1270 if ((d2 - unorm) > 0.)
1274 for (l = 0; l <
m; l++)
1276 _nnls_h12(2, npp1, npp1 + 1, m, &a[j*m+0], 1, &up, zz, 1, 1, 1);
1277 ztest = zz[npp1] / a[j*m+npp1];
1285 a[j*m+npp1] = asave;
1294 for (l = 0; l <
m; ++l)
1296 index[iz] = index[iz1];
1302 for (jz = iz1; jz <= iz2; jz++)
1305 _nnls_h12(2, nsetp - 1, npp1, m, &a[j*m+0], 1, &up,
1306 &a[jj*m+0], 1, m, 1);
1309 for (l = npp1; l <
m; l++)
1313 for (l = 0; l < nsetp; l++)
1315 ip = nsetp - (l + 1);
1317 for (ii = 0; ii <= ip; ii++)
1318 zz[ii] -= a[jj*m+ii] * zz[ip+1];
1320 zz[ip] /= a[jj*m+ip];
1324 while (++iter < itmax)
1327 for (alpha = 2.0, ip = 0; ip < nsetp; ip++)
1332 t = -x[l] / (zz[ip] - x[l]);
1346 for (ip = 0; ip < nsetp; ip++)
1349 x[l] += alpha * (zz[ip] - x[l]);
1359 if (jj != (nsetp - 1))
1362 for (j = jj + 1; j < nsetp; j++)
1366 _nnls_g1(a[ii*m+j-1], a[ii*m+j], &cc, &ss, &a[ii*m+j-1]);
1367 for (a[ii*m+j] = 0., l = 0; l <
n; l++)
1372 a[l*m+j-1] = cc * temp + ss * a[l*m+
j];
1373 a[l*m+
j] = -ss * temp + cc * a[l*m+
j];
1377 b[j-1] = cc * temp + ss * b[
j];
1378 b[
j] = -ss * temp + cc * b[
j];
1390 for (jj = 0; jj < nsetp; jj++)
1403 for (k = 0; k <
m; k++)
1405 for (l = 0; l < nsetp; l++)
1407 ip = nsetp - (l + 1);
1409 for (ii = 0; ii <= ip; ii++)
1410 zz[ii] -= a[jj*m+ii] * zz[ip+1];
1412 zz[ip] /= a[jj*m+ip];
1420 for (ip = 0; ip < nsetp; ip++)
1429 for (k = npp1; k <
m; k++)
1430 sm += (b[k] * b[k]);
1432 for (j = 0; j <
n; j++)
1455 int nnlsWght(
int N,
int M,
double *A,
double *
b,
double *weight)
1461 if (N < 1 || M < 1 || A == NULL || b == NULL || weight == NULL)
1465 w = (
double*)malloc(M *
sizeof(
double));
1470 for (m = 0; m < M; m++)
1472 if (weight[m] <= 1.0e-20)
1475 w[
m] =
sqrt(weight[m]);
1479 for (m = 0; m < M; m++)
1481 for (n = 0; n < N; n++)
1501 return(absa *
sqrt(1.0 + absb * absb / (absa * absa)));
1503 return((absb == 0.0) ? (0.0) : (absb *
sqrt(1.0 + absa * absa / (absb * absb))));
1506 #define SVDMAXITER 1000 1507 void svdcmp(
double *U,
int Lines,
int Columns,
double *W,
double *V)
1510 double c,
f,
g, h, s;
1512 long i, its,
j, jj,
k, l = 0L, nm = 0L;
1516 std::vector<double>
buffer(Columns*Columns);
1517 auto *rv1= buffer.data();
1518 g = Scale = Norm = 0.0;
1519 for (i = 0L; (i < Columns); i++)
1523 g = s = Scale = 0.0;
1526 for (k = i; (k < Lines); k++)
1528 Scale += fabs(U[k * Columns + i]);
1532 for (k = i; (k < Lines); k++)
1534 U[k * Columns +
i] /= Scale;
1535 s += U[k * Columns +
i] * U[k * Columns +
i];
1537 f = U[i * Columns +
i];
1538 g = (0.0 <=
f) ? (-
sqrt(s)) : (
sqrt(s));
1540 U[i * Columns +
i] = f -
g;
1541 for (j = l; (j < Columns); j++)
1543 for (s = 0.0, k = i; (k < Lines); k++)
1545 s += U[k * Columns +
i] * U[k * Columns +
j];
1548 for (k = i; (k < Lines); k++)
1550 U[k * Columns +
j] += f * U[k * Columns +
i];
1553 for (k = i; (k < Lines); k++)
1555 U[k * Columns +
i] *= Scale;
1560 g = s = Scale = 0.0;
1561 if ((i < Lines) && (i != (Columns - 1L)))
1563 for (k = l; (k < Columns); k++)
1565 Scale += fabs(U[i * Columns + k]);
1569 for (k = l; (k < Columns); k++)
1571 U[i * Columns +
k] /= Scale;
1572 s += U[i * Columns +
k] * U[i * Columns +
k];
1574 f = U[i * Columns + l];
1575 g = (0.0 <=
f) ? (-
sqrt(s)) : (
sqrt(s));
1577 U[i * Columns + l] = f -
g;
1578 for (k = l; (k < Columns); k++)
1580 rv1[
k] = U[i * Columns +
k] / h;
1582 for (j = l; (j < Lines); j++)
1584 for (s = 0.0, k = l; (k < Columns); k++)
1586 s += U[j * Columns +
k] * U[i * Columns +
k];
1588 for (k = l; (k < Columns); k++)
1590 U[j * Columns +
k] += s * rv1[
k];
1593 for (k = l; (k < Columns); k++)
1595 U[i * Columns +
k] *= Scale;
1599 Norm = ((fabs(W[i]) + fabs(rv1[i])) < Norm) ? (Norm) : (fabs(W[i]) + fabs(rv1[i]));
1601 for (i = Columns - 1L; (0L <=
i); i--)
1603 if (i < (Columns - 1L))
1607 for (j = l; (j < Columns); j++)
1609 V[j * Columns +
i] = U[i * Columns +
j] / (U[i * Columns + l] *
g);
1611 for (j = l; (j < Columns); j++)
1613 for (s = 0.0, k = l; (k < Columns); k++)
1615 s += U[i * Columns +
k] * V[k * Columns +
j];
1617 for (k = l; (k < Columns); k++)
1621 V[k * Columns +
j] += s * V[k * Columns +
i];
1626 for (j = l; (j < Columns); j++)
1628 V[i * Columns +
j] = V[j * Columns +
i] = 0.0;
1631 V[i * Columns +
i] = 1.0;
1635 for (i = (Lines < Columns) ? (Lines - 1L) : (Columns - 1L); (0L <=
i); i--)
1639 for (j = l; (j < Columns); j++)
1641 U[i * Columns +
j] = 0.0;
1646 for (j = l; (j < Columns); j++)
1648 for (s = 0.0, k = l; (k < Lines); k++)
1650 s += U[k * Columns +
i] * U[k * Columns +
j];
1652 f = s * g / U[i * Columns +
i];
1653 for (k = i; (k < Lines); k++)
1657 U[k * Columns +
j] += f * U[k * Columns +
i];
1661 for (j = i; (j < Lines); j++)
1663 U[j * Columns +
i] *=
g;
1668 for (j = i; (j < Lines); j++)
1670 U[j * Columns +
i] = 0.0;
1673 U[i * Columns +
i] += 1.0;
1675 for (k = Columns - 1L; (0L <=
k); k--)
1677 for (its = 1L; (its <= MaxIterations); its++)
1680 for (l = k; (0L <= l); l--)
1683 if ((fabs(rv1[l]) + Norm) == Norm)
1688 if ((fabs(W[nm]) + Norm) == Norm)
1697 for (i = l; (i <=
k); i++)
1701 if ((fabs(f) + Norm) == Norm)
1711 for (j = 0L; (j < Lines); j++)
1713 y = U[j * Columns + nm];
1714 z = U[j * Columns +
i];
1715 U[j * Columns + nm] = y * c + z * s;
1716 U[j * Columns +
i] = z * c - y * s;
1726 for (j = 0L; (j < Columns); j++)
1728 V[j * Columns +
k] = -V[j * Columns +
k];
1733 if (its == MaxIterations)
return;
1739 f = ((y -
z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);
1741 f = ((x -
z) * (x + z) + h * ((y / (f + ((0.0 <=
f) ? (fabs(g))
1742 : (-fabs(g))))) - h)) /
x;
1744 for (j = l; (j <= nm); j++)
1759 for (jj = 0L; (jj < Columns); jj++)
1761 x = V[jj * Columns +
j];
1762 z = V[jj * Columns +
i];
1763 V[jj * Columns +
j] = x * c + z * s;
1764 V[jj * Columns +
i] = z * c - x * s;
1776 for (jj = 0L; (jj < Lines); jj++)
1778 y = U[jj * Columns +
j];
1779 z = U[jj * Columns +
i];
1780 U[jj * Columns +
j] = y * c + z * s;
1781 U[jj * Columns +
i] = z * c - y * s;
1791 void svbksb(
double *
u,
double *
w,
double *v,
int m,
int n,
double *
b,
double *
x)
1796 std::vector<double>
buffer(n);
1797 auto *tmp= buffer.data()-1;
1798 for (j = 1;j <=
n;j++)
1803 for (i = 1;i <=
m;i++)
1804 s += u[i*n+j] * b[i];
1809 for (j = 1;j <=
n;j++)
1812 for (jj = 1;jj <=
n;jj++)
1813 s += v[j*n+jj] * tmp[jj];
1973 #define Extern extern 2108 #define abs(x) ((x) >= 0 ? (x) : -(x)) 2109 #define dabs(x) (doublereal)abs(x) 2110 #define min(a,b) ((a) <= (b) ? (a) : (b)) 2111 #define max(a,b) ((a) >= (b) ? (a) : (b)) 2112 #define dmin(a,b) (doublereal)min(a,b) 2113 #define dmax(a,b) (doublereal)max(a,b) 2117 #define F2C_proc_par_types 1 2119 typedef int (*U_fp)(...);
2129 typedef int(*
S_fp)(...);
2131 typedef int (*U_fp)();
2152 #define Skip_f2c_Undefs 2154 #ifndef Skip_f2c_Undefs 2186 #define cmache_1 cmache_ 2221 double *
c,
double *
d,
double *
a,
double *
b,
double *
xl,
2227 int ql0001_(m, me, mmax, n, nmax, mnn, c, d, a, b, xl, xu, x,
2228 u, iout, ifail, iprint, war, lwar, iwar, liwar, eps1)
2265 a_offset = a_dim1 + 1;
2269 c_offset = c_dim1 + 1;
2279 if (fabs(c[*nmax + *nmax * c_dim1]) == 0.
e0)
2297 maxit = (*m + *
n) * 40;
2300 inw2 = inw1 + *
mmax;
2317 lw = *nmax * 3 * *nmax / 2 + *nmax * 10 + *
m;
2318 if (inw2 + lw > *lwar)
2326 if (*mnn < *m + *n + *n)
2334 ql0002_(n, m, me, mmax, &mn, mnn, nmax, &lql, &a[a_offset], &war[inw1], &
2335 d[1], &c[c_offset], &xl[1], &xu[1], &x[1], &nact, &iwar[1], &
2336 maxit, &qpeps, &info, &diag, &war[inw2], &lw);
2357 for (j = 1; j <=
i__1; ++
j)
2380 *ifail = -info + 10;
2469 int ql0002_(n, m, meq, mmax, mn, mnn, nmax, lql, a, b, grad,
2470 g, xl, xu, x, nact, iact, maxit, vsmall, info, diag, w, lw)
2490 static doublereal onha, xmag, suma, sumb, sumc, temp, step,
zero;
2495 static integer iflag, jflag, kflag, lflag;
2497 static integer ifinc, kfinc, jfinc, mflag, nflag;
2507 static integer ii, il, kk, jl,
ir, nm, is, iu, iw, ju, ix, iz, nu, iy;
2561 g_offset = g_dim1 + 1;
2566 a_offset = a_dim1 + 1;
2571 iwr = iwz + *nmax * *
nmax;
2572 iww = iwr + *nmax * (*nmax + 3) / 2;
2597 kfinc =
max(10, *n);
2612 for (i = 1; i <=
i__2; ++
i)
2616 d__1 = a[k + i *
a_dim1];
2641 sum = one /
sqrt(sum);
2649 for (k = 1; k <=
i__1; ++
k)
2664 for (i = 1; i <=
i__1; ++
i)
2667 w[id] = g[i + i *
g_dim1];
2669 d__1 = *diag, d__2 = *vsmall - w[id];
2670 *diag =
max(d__1, d__2);
2677 for (j = ii; j <=
i__2; ++
j)
2680 d__1 = w[id], d__2 = g[j + j *
g_dim1];
2681 ga = -
min(d__1, d__2);
2682 gb = (d__1 = w[id] - g[j + j *
g_dim1],
abs(d__1)) + (d__2 = g[i
2683 + j * g_dim1],
abs(d__2));
2687 d__1 = g[i + j *
g_dim1];
2688 ga += d__1 * d__1 / gb;
2691 *diag =
max(*diag, ga);
2701 *diag = diagr * *diag;
2703 for (i = 1; i <=
i__1; ++
i)
2707 g[i + i *
g_dim1] = *diag + w[id];
2716 for (j = 1; j <=
i__1; ++
j)
2721 for (i = 1; i <=
i__2; ++
i)
2723 temp = g[i + j *
g_dim1];
2729 for (k = irb; k <= i__3; ++
k)
2733 temp -= w[
k] * w[ira];
2740 w[
ir] = temp / w[ira];
2763 for (i = k; i <=
i__1; ++
i)
2765 sum -= w[ira] * w[
i];
2774 sumx += d__1 *
d__1;
2779 *diag = *diag + *vsmall - temp / sumx;
2788 for (i = 1; i <=
i__1; ++
i)
2791 for (j = 1; j <=
i__2; ++
j)
2804 for (i = 1; i <=
i__2; ++
i)
2812 for (j = 2; j <=
i__1; ++
j)
2819 ir = iwr + (i + i *
i) / 2;
2820 w[iz] = one / w[
ir];
2827 for (j = i; j <=
i__1; ++
j)
2833 for (k = iza; i__4 < 0 ? k >= i__3 : k <= i__3; k += i__4)
2835 sum += w[
k] * w[
ir];
2841 w[iz] = -sum / w[
ir];
2863 for (i = 1; i <=
i__2; ++
i)
2912 for (i = 1; i <=
i__2; ++
i)
2923 for (j = i; j <=
i__1; ++
j)
2926 w[id] += g[i + j *
g_dim1] * x[
j];
2929 for (j = 1; j <=
i__1; ++
j)
2933 w[iw] += g[j + i *
g_dim1] * w[id];
2938 for (j = 1; j <=
i__1; ++
j)
2941 w[iw] += g[i + j *
g_dim1] * x[
j];
2951 for (k = 1; k <=
i__2; ++
k)
2961 for (i = 1; i <=
i__1; ++
i)
2964 w[iw] -= w[
k] * a[kk + i *
a_dim1];
2966 w[is] -= x[
i] * a[kk + i *
a_dim1];
2977 w[is] = xl[k1] - x[k1];
2983 w[is] = -xu[k1] + x[k1];
2996 for (i = il; i <=
i__2; ++
i)
3005 for (j = il; j <=
i__1; ++
j)
3009 sum += w[
ir] * w[
j];
3014 w[
i] = (w[
i] - sum) / w[ir];
3021 for (i = 1; i <=
i__2; ++
i)
3026 for (j = il; j <=
i__1; ++
j)
3028 sum += w[
j] * w[iz];
3040 for (j = i; j <=
i__1; ++
j)
3043 w[id] += g[i + j *
g_dim1] * sum;
3047 for (j = 1; j <=
i__1; ++
j)
3051 w[iw] += g[j + i *
g_dim1] * w[id];
3056 for (j = 1; j <=
i__1; ++
j)
3060 w[iw] += sum * g[i + j *
g_dim1];
3081 il = iws + *nact + 1;
3082 iza = iwz + *nact * *
n;
3084 for (i = 1; i <=
i__2; ++
i)
3089 for (j = il; j <=
i__1; ++
j)
3091 sum += w[iz] * w[
j];
3111 for (k = 1; k <=
i__2; ++
k)
3137 if (w[kdrop] >= zero)
3141 if (iact[kdrop] <= *meq)
3165 for (k = 1; k <=
i__2; ++
k)
3174 for (i = 1; i <=
i__1; ++
i)
3177 sum += x[
i] * a[k + i *
a_dim1];
3179 sumx = -sum * w[ia];
3188 temp = (d__1 = b[
k],
abs(d__1));
3190 for (i = 1; i <=
i__1; ++
i)
3193 temp += (d__1 = x[
i] * a[k + i *
a_dim1],
abs(d__1));
3195 tempa = temp +
abs(sum);
3200 temp += onha *
abs(sum);
3213 for (k = 1; k <=
i__2; ++
k)
3257 if (cvmax <= *vsmall)
3277 for (i = 1; i <=
i__2; ++
i)
3279 sum = two * grad[
i];
3288 for (j = i; j <=
i__1; ++
j)
3292 w[id] += g[i + j *
g_dim1] * (w[ix] + x[
j]);
3295 for (j = 1; j <=
i__1; ++
j)
3298 temp = g[j + i *
g_dim1] * w[id];
3306 for (j = 1; j <=
i__1; ++
j)
3309 temp = g[i + j *
g_dim1] * (w[ix] + x[
j]);
3316 fdiff += sum * (x[
i] - w[ix]);
3318 fdiffa += sumx * (d__1 = x[
i] - w[ix],
abs(d__1));
3321 sum = fdiffa + fdiff;
3326 temp = fdiffa + onha * fdiff;
3335 for (i = 1; i <=
i__2; ++
i)
3348 if (iterc <= *maxit)
3355 iws = iwr + (*nact + *nact * *nact) / 2;
3361 for (i = 1; i <=
i__2; ++
i)
3365 w[iw] = a[knext + i *
a_dim1];
3370 for (i = 1; i <=
i__2; ++
i)
3385 for (i = 1; i <=
i__2; ++
i)
3399 for (i = 1; i <=
i__2; ++
i)
3435 iz = iwz + *nact * *
n;
3437 for (i = 1; i <=
i__2; ++
i)
3441 suma += w[iw] * w[iz];
3442 sumb += (d__1 = w[iw] * w[iz],
abs(d__1));
3446 sumc += d__1 *
d__1;
3448 temp = sumb +
abs(suma) * .1;
3449 tempa = sumb +
abs(suma) * .2;
3470 temp = sumc +
abs(suma) * .1;
3471 tempa = sumc +
abs(suma) * .2;
3501 for (i = 1; i <=
i__2; ++
i)
3503 suma = a[knext + i *
a_dim1];
3510 for (k = 1; k <=
i__1; ++
k)
3526 temp = -w[iww + kk];
3531 temp = w[iw] * a[kk + i *
a_dim1];
3538 if (suma <= *vsmall)
3542 temp = sumb +
abs(suma) * .1;
3543 tempa = sumb +
abs(suma) * .2;
3565 for (i = 1; i <=
i__2; ++
i)
3584 for (k = 1; k <=
i__1; ++
k)
3600 temp = -w[iww + kk];
3605 temp = w[iw] * a[kk + i *
a_dim1];
3612 temp = sumb +
abs(suma) * .1;
3613 tempa = sumb +
abs(suma) * .2;
3648 for (k = 1; k <=
i__2; ++
k)
3651 w[
k] -= parinc * w[iw];
3655 d__1 =
zero, d__2 = w[
k];
3656 w[
k] =
max(d__1, d__2);
3674 iws = iws - *nact - 1;
3677 for (i = 1; i <=
i__2; ++
i)
3694 parinc = step / sumy;
3710 temp = one - ratio / parinc;
3719 step = ratio * sumy;
3727 iwy = iwz + *nact * *
n;
3729 for (i = 1; i <=
i__2; ++
i)
3733 x[
i] += step * w[iy];
3746 iact[*nact] = knext;
3760 if (sum < xmagr * xmag)
3793 for (i = 1; i <=
i__2; ++
i)
3797 g[i + i *
g_dim1] = w[id];
3814 ir = iwr + (*nact + *nact * *nact) / 2;
3826 for (j = i; j <=
i__2; ++
j)
3829 sum += w[ira] * w[iw];
3839 w[iw] = (w[is] - sum) / w[ir];
3862 for (k = 1; k <=
i__2; ++
k)
3864 if (iact[k] <= *meq)
3869 if (res * w[iw] >= zero)
3873 temp = w[
k] / w[iw];
3878 if (
abs(temp) >=
abs(ratio))
3905 ia = iwa + iact[kdrop];
3906 if (iact[kdrop] > *mn)
3919 iz = iwz + kdrop * *
n;
3920 ir = iwr + (kdrop + kdrop * kdrop) / 2;
3923 ir = ir + kdrop + 1;
3925 d__3 = (d__1 = w[ir - 1],
abs(d__1)), d__4 = (d__2 = w[ir],
abs(d__2));
3926 temp =
max(d__3, d__4);
3928 d__1 = w[ir - 1] / temp;
3930 d__2 = w[
ir] / temp;
3931 sum = temp *
sqrt(d__1 * d__1 + d__2 * d__2);
3932 ga = w[ir - 1] / sum;
3938 for (i = 1; i <=
i__2; ++
i)
3954 for (i = kdrop; i <=
i__2; ++
i)
3956 temp = ga * w[ira] + gb * w[ira + 1];
3957 w[ira + 1] = ga * w[ira + 1] - gb * w[ira];
3966 for (i = 1; i <=
i__2; ++
i)
3970 temp = ga * w[
j] + gb * w[iz];
3971 w[iz] = ga * w[iz] - gb * w[
j];
3978 iact[kdrop - 1] = iact[kdrop];
3979 w[kdrop - 1] = w[kdrop];
4018 d__3 = (d__1 = w[is - 1],
abs(d__1)), d__4 = (d__2 = w[is],
abs(d__2));
4019 temp =
max(d__3, d__4);
4021 d__1 = w[is - 1] / temp;
4023 d__2 = w[is] / temp;
4024 sum = temp *
sqrt(d__1 * d__1 + d__2 * d__2);
4025 ga = w[is - 1] / sum;
4029 for (i = 1; i <=
i__2; ++
i)
4032 temp = ga * w[iz] + gb * w[
k];
4033 w[
k] = ga * w[
k] - gb * w[iz];
4058 for (i = 1; i <=
i__2; ++
i)
4060 sum += (d__1 = x[
i],
abs(d__1)) * vfact * ((d__2 = grad[i],
abs(d__2))
4061 + (d__3 = g[i + i *
g_dim1] * x[
i],
abs(d__3)));
4077 xmag =
max(xmag, sum);
4097 for (i = 1; i <=
i__2; ++
i)
4103 for (j = jl; j <=
i__1; ++
j)
4107 w[is] += w[iz] * w[
j];
4121 comments from the converter:
4222 static int *make_iv(
int);
4223 static double *make_dv(
int);
4224 static double **make_dm(
int,
int);
4226 static void free_dv(
double *);
4227 static void free_dm(
double **,
int);
4228 static double *convert(
double **,
int,
int);
4230 static int *make_iv();
4231 static double *make_dv();
4232 static double **make_dm();
4236 static double *convert();
4245 ql0001_(
int *,
int *,
int *,
int *,
int *,
int *,
double *,
double *,
4246 double *,
double *,
double *,
double *,
double *,
double *,
4247 int *,
int *,
int *,
double *,
int *,
int *,
int *,
double *);
4248 static void diagnl(
int,
double,
double **);
4249 static void error(
const char string[],
int *);
4251 estlam(
int,
int,
int *,
double,
double **,
double *,
double *,
double *,
4252 struct _constraint *,
double *,
double *,
double *,
double *);
4253 static double *colvec(
double **,
int,
int);
4254 static double scaprd(
int,
double *,
double *);
4255 static double smallNumber();
4256 static int fuscmp(
double,
double);
4257 static int indexs(
int,
int);
4258 static void matrcp(
int,
double **,
int,
double **);
4259 static void matrvc(
int,
int,
double **,
double *,
double *);
4260 static void nullvc(
int,
double *);
4262 resign(
int,
int,
double *,
double *,
double *,
struct _constraint *,
4263 double *,
int,
int);
4264 static void sbout1(FILE *,
int,
const char *,
double,
double *,
int,
int);
4265 static void sbout2(FILE *,
int,
int,
const char *,
const char *,
double *);
4266 static void shift(
int,
int,
int *);
4268 slope(
int,
int,
int,
int,
int,
struct _objective *,
double *,
double *,
4269 double *,
double,
double,
int,
double *,
int);
4270 static int element(
int *,
int,
int);
4274 static void error();
4275 static void estlam();
4276 static double *colvec();
4277 static double scaprd();
4278 static double smallNumber();
4279 static int fuscmp();
4280 static int indexs();
4284 static void resign();
4286 static void sbout2();
4287 static void shift();
4288 static double slope();
4289 static int element();
4297 void grobfd(
int,
int,
double *,
double *,
void(*)(
int,
int,
4298 double *,
double *,
void *),
void *);
4299 void grcnfd(
int,
int,
double *,
double *,
void(*)(
int,
int,
4300 double *,
double *,
void *),
void *);
4312 cfsqp1(
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int *,
int,
4313 int,
int,
int,
double,
double,
int *,
int *,
struct _parameter *,
4315 void(*)(
int,
int,
double *,
double *,
void *),
4316 void(*)(
int,
int,
double *,
double *,
void *),
4317 void(*)(
int,
int,
double *,
double *,
4318 void(*)(
int,
int,
double *,
double *,
void *),
void *),
4319 void(*)(
int,
int,
double *,
double *,
4320 void(*)(
int,
int,
double *,
double *,
void *),
void *));
4322 check(
int,
int,
int,
int *,
int,
int,
int,
int,
int,
int,
int,
int *,
double,
4325 initpt(
int,
int,
int,
int,
int,
int,
int,
struct _parameter *,
4326 struct _constraint *,
void(*)(
int,
int,
double *,
double *,
void *),
4327 void(*)(
int,
int,
double *,
double *,
4328 void(*)(
int,
int,
double *,
double *,
void *),
void *));
4330 dir(
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
double *,
4331 double,
double,
double *,
double *,
double,
double *,
double *,
int *,
4332 int *,
int *,
int *,
int *,
int *,
struct _parameter *,
double *,
4334 double *,
double *,
double *,
double *,
double *,
double **,
double *,
4335 double *,
double *,
double *,
double **,
double **,
double *,
4336 double *,
struct _violation *,
void(*)(
int,
int,
double *,
double *,
4337 void *),
void(*)(
int,
int,
double *,
double *,
void *));
4339 step1(
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int *,
int *,
int *,
4340 int *,
int *,
int *,
int *,
int *,
int,
double,
struct _objective *,
4341 double *,
double *,
double *,
double *,
double *,
double *,
double *,
4342 double *,
double *,
double *,
double *,
struct _constraint *,
4344 void(*)(
int,
int,
double *,
double *,
void *),
4345 void(*)(
int,
int,
double *,
double *,
void *),
void *);
4347 hessian(
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int *,
int,
4349 double,
double *,
double *,
double *,
double *,
double *,
4350 struct _constraint *,
double *,
int *,
int *,
double *,
4351 double *,
double *,
double **,
double *,
double,
int *,
4352 double *,
double *,
void(*)(
int,
int,
double *,
double *,
void *),
4353 void(*)(
int,
int,
double *,
double *,
void *),
4354 void(*)(
int,
int,
double *,
double *,
4355 void(*)(
int,
int,
double *,
double *,
void *),
void *),
4356 void(*)(
int,
int,
double *,
double *,
4357 void(*)(
int,
int,
double *,
double *,
void *),
void *),
4358 double **,
double *,
double *,
struct _violation *);
4360 out(
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int *,
double *,
4362 double,
double,
double,
double,
int);
4364 update_omega(
int,
int,
int,
int *,
int,
int,
int,
int,
double,
double,
4366 struct _violation *,
void(*)(
int,
int,
double *,
double *,
4367 void *),
void(*)(
int,
int,
double *,
double *,
void *),
4368 void(*)(
int,
int,
double *,
double *,
4369 void(*)(
int,
int,
double *,
double *,
void *),
void *),
4370 void(*)(
int,
int,
double *,
double *,
4371 void(*)(
int,
int,
double *,
double *,
void *),
void *),
4375 static void check();
4376 static void initpt();
4378 static void step1();
4379 static void hessian();
4381 static void update_omega();
4398 double *x,
double *
f,
double *g,
double *
lambda,
4399 void(*
obj)(
int,
int,
double *,
double *,
void *),
4400 void(*
constr)(
int,
int,
double *,
double *,
void *),
4401 void(*
gradob)(
int,
int,
double *,
double *,
4402 void(*)(
int,
int,
double *,
double *,
void *),
void *),
4403 void(*
gradcn)(
int,
int,
double *,
double *,
4404 void(*)(
int,
int,
double *,
double *,
void *),
void *),
4408 cfsqp(nparam, nf, nfsr, nineqn, nineq, neqn, neq, ncsrl, ncsrn, mesh_pts,
4409 mode, iprint, miter, inform, bigbnd, eps, epseqn, udelta, bl, bu, x,
4411 int nparam, nf, nfsr, neqn, nineqn, nineq, neq, ncsrl, ncsrn, mode,
4412 iprint, miter, *mesh_pts, *inform;
4413 double bigbnd, eps, epseqn, udelta;
4414 double *bl, *bu, *x, *f, *g, *lambda;
4630 int i,
ipp,
j,
ncnstr,
nclin,
nctotl,
nob,
nobL,
modem=0,
nn,
4643 mesh_pts = mesh_pts - 1;
4650 for (i = 1; i <=
nfsip1; i++)
4651 nfsr = nfsr + mesh_pts[i];
4653 nineqn = nineqn -
ncsrn;
4654 nineq = nineq - ncsrl -
ncsrn;
4660 for (i = 1; i <=
ncsipn1; i++)
4661 ncsrn = ncsrn + mesh_pts[nfsip1+i];
4663 for (i = 1; i <=
ncsipl1; i++)
4664 ncsrl = ncsrl + mesh_pts[nfsip1+ncsipn1+i];
4665 nineqn = nineqn +
ncsrn;
4666 nineq = nineq + ncsrn +
ncsrl;
4668 cs = (
struct _constraint *)calloc(nineq + neq + 1,
4670 for (i = 1; i <= nineq +
neq; i++)
4672 cs[
i].
grad = make_dv(nparam);
4677 _param.
x = make_dv(nparam + 1);
4687 for (i = 1; i <=
nparam; i++)
4690 param->
bl[
i] = bl[
i];
4691 param->
bu[
i] = bu[
i];
4697 lambda = lambda - 1;
4708 signeq = make_dv(neqn);
4722 "\n\n CFSQP Version 2.5b (Released June 1997) \n");
4724 " Copyright (c) 1993 --- 1997 \n");
4726 " C.T. Lawrence, J.L. Zhou \n");
4728 " and A.L. Tits \n");
4730 " All Rights Reserved \n\n");
4735 check(nparam, nf, nfsr, &Linfty, nineq, nineqn, neq, neqn,
4736 ncsrl, ncsrn, mode, &modem, eps, bgbnd, param);
4747 nppram = nparam + 1;
4752 for (i = 1; i <=
nparam; i++)
4755 if (param->
bl[i] <= xi && param->
bu[i] >= xi)
4766 for (i = 1; i <=
nclin; i++)
4771 constr(nparam, j, (param->
x) + 1, &gi, param->
cd);
4777 constr(nparam, j + neqn, (param->
x) + 1, &gi, param->
cd);
4792 " The given initial point is infeasible for inequality\n");
4794 " constraints and linear equality constraints:\n");
4799 nctotl = nparam +
nclin;
4800 lenw = 2 * nparam * nparam + 10 * nparam + 2 * nctotl + 1;
4801 leniw =
DMAX1(2 * nparam + 2 * nctotl + 3, 2 * nclin + 2 * nparam + 6);
4806 nrowa =
DMAX1(nclin, 1);
4807 iw = make_iv(leniw);
4809 initpt(nparam, nineqn, neq, neqn, nclin, nctotl, nrowa, param,
4810 &cs[nineqn], constr, gradcn);
4819 indxob = make_iv(
DMAX1(nineq + neq, nf));
4820 indxcn = make_iv(nineq + neq);
4826 constr(nparam, i, (param->
x) + 1, &(cs[i].
val), param->
cd);
4827 if (cs[i].
val > 0.
e0)
4836 for (i = 1; i <=
nineqn; i++)
4838 ob[
i].
grad = make_dv(nparam);
4841 for (i = 1; i <=
nineqn; i++)
4848 for (i = 1; i <= nineq -
nineqn; i++)
4849 indxcn[i] = nineqn + i;
4850 for (i = 1; i <= neq -
neqn; i++)
4851 indxcn[i+nineq-nineqn] = nineq + neqn + i;
4859 for (i = 1; i <=
nf; i++)
4861 ob[
i].
grad = make_dv(nparam);
4864 for (i = 1; i <=
nineqn; i++)
4868 for (i = 1; i <= neq -
neqn; i++)
4869 cs[i+nineq+neqn].
val = cs[i+nineq].
val;
4870 for (i = 1; i <=
neqn; i++)
4873 constr(nparam, j, (param->
x) + 1, &(cs[j].
val), param->
cd);
4874 indxcn[nineqn+
i] =
j;
4876 for (i = 1; i <= nineq -
nineqn; i++)
4877 indxcn[i+nn] = nineqn + i;
4878 for (i = 1; i <= neq -
neqn; i++)
4879 indxcn[i+nineq+neqn] = nineq + neqn + i;
4886 "The given initial point is feasible for inequality\n");
4888 " constraints and linear equality constraints:\n");
4900 "To generate a feasible point for nonlinear inequality\n");
4902 "constraints and linear equality constraints, ");
4913 if (feasb && !feasbl)
4916 "Starting from the generated point feasible for\n");
4918 "inequality constraints and linear equality constraints:\n");
4920 dummy, param->
x, 2, 1);
4926 "Starting from the generated point feasible for\n");
4928 "inequality constraints and linear equality constraints:\n");
4930 dummy, param->
x, 2, 1);
4936 if (ipp > 0 && !feasb && !prnt)
4939 " The given initial point is infeasible for inequality\n");
4941 " constraints and linear equality constraints:\n");
4961 if (nob != 0 || neqn != 0)
4964 "current feasible iterate with no objective specified\n");
4966 for (i = 1; i <= nineq +
neq; i++)
4968 dealloc(nineq, neq, signeq, indxcn, indxob, cs, param);
4978 nctotl = nppram + ncnstr +
DMAX1(nobL, 1);
4979 leniw = 2 * (ncnstr +
DMAX1(nobL, 1)) + 2 * nppram + 6;
4980 lenw = 2 * nppram * nppram + 10 * nppram + 6 * (ncnstr +
DMAX1(nobL, 1) + 1);
4982 if (modem == 1 && nn == 0)
4985 param->
x[nparam+1] =
gmax;
4988 for (i = 1; i <=
neqn; i++)
4990 if (cs[i+nineq].
val > 0.
e0)
5001 mesh_pts1 = &mesh_pts[
nfsr];
5016 iw = make_iv(leniw);
5018 cfsqp1(miter, nparam, nob, nobL, nfsip1, nineqn, neq, neqn, ncsipl1, ncsipn1,
5019 mesh_pts1, ncnstr, nctotl, nrowa, feasb, epskt, epseqn, indxob,
5020 indxcn, param, cs, ob, signeq, obj, constr, gradob, gradcn);
5026 for (i = 1; i <=
nob; i++)
5029 for (i = 1; i <=
nineqn; i++)
5038 for (i = 1; i <=
nparam; i++)
5040 for (i = 1; i <= nineq +
neq; i++)
5043 dealloc(nineq, neq, signeq, indxcn, indxob, cs, param);
5044 for (i = 1; i <=
nf; i++)
5054 "Error: No feasible point is found for nonlinear inequality\n");
5056 "constraints and linear equality constraints\n");
5058 dealloc(nineq, neq, signeq, indxcn, indxob, cs, param);
5059 for (i = 1; i <=
nineqn; i++)
5066 for (i = 1; i <=
nparam; i++)
5069 lambda[
i] = param->
mult[
i];
5071 for (i = 1; i <= nineq +
neq; i++)
5076 for (i = 1; i <=
nf; i++)
5079 lambda[i+nparam+nineq+
neq] = ob[
i].
mult;
5084 lambda[1+nparam+nineq+
neq] = 1.e0;
5086 dealloc(nineq, neq, signeq, indxcn, indxob, cs, param);
5096 dealloc(
int nineq,
int neq,
double *signeq,
int *indxob,
5100 dealloc(nineq, neq, signeq, indxob, indxcn, cs, param)
5117 for (i = 1; i <= nineq +
neq; i++)
5128 dealloc1(
int,
int,
double **,
double **,
double **,
double *,
double *,
5129 double *,
double *,
double *,
double *,
double *,
double *,
5130 double *,
double *,
double *,
double *,
int *,
int *,
int *);
5132 static void dealloc1();
5137 cfsqp1(
int miter,
int nparam,
int nob,
int nobL,
int nfsip,
int nineqn,
5138 int neq,
int neqn,
int ncsipl,
int ncsipn,
int *mesh_pts,
int ncnstr,
5139 int nctotl,
int nrowa,
int feasb,
double epskt,
double epseqn,
5140 int *indxob,
int *indxcn,
struct _parameter *param,
5142 double *signeq,
void(*obj)(
int,
int,
double *,
double *,
void *),
5143 void(*constr)(
int,
int,
double *,
double *,
void *),
5144 void(*gradob)(
int,
int,
double *,
double *,
5145 void(*)(
int,
int,
double *,
double *,
void *),
void *),
5146 void(*gradcn)(
int,
int,
double *,
double *,
5147 void(*)(
int,
int,
double *,
double *,
void *),
void *))
5150 cfsqp1(miter, nparam, nob, nobL, nfsip, nineqn, neq, neqn, ncsipl, ncsipn,
5151 mesh_pts, ncnstr, nctotl, nrowa, feasb, epskt, epseqn, indxob,
5152 indxcn, param, cs, ob, signeq, obj, constr, gradob, gradcn)
5153 int miter, nparam, nob, nobL, nfsip, nineqn, neq, neqn, ncnstr,
5154 nctotl, nrowa, feasb, ncsipl, ncsipn, *mesh_pts;
5155 int *indxob, *indxcn;
5156 double epskt, epseqn;
5161 void(* obj)(), (* constr)(), (* gradob)(), (* gradcn)();
5177 hess = make_dm(nparam, nparam);
5178 hess1 = make_dm(nparam + 1, nparam + 1);
5179 a = make_dm(nrowa, nparam + 2);
5180 di = make_dv(nparam + 1);
5181 d = make_dv(nparam + 1);
5182 gm = make_dv(4 * neqn);
5183 grdpsf = make_dv(nparam);
5184 penp = make_dv(neqn);
5185 bl = make_dv(nctotl);
5186 bu = make_dv(nctotl);
5187 clamda = make_dv(nctotl + nparam + 1);
5188 cvec = make_dv(nparam + 1);
5189 psmu = make_dv(neqn);
5191 backup = make_dv(nob + ncnstr);
5192 iact = make_iv(nob + nineqn + neqn);
5194 istore = make_iv(nineqn + nob + neqn);
5237 for (i = 1; i <=
ncnst1; i++)
5241 if (feasb && i > nineqn && i <= nn)
5246 if (feasb && i <= nn)
5248 iact[
i] = indxcn[
i];
5253 gradcn(nparam, indxcn[i], (param->
x) + 1, (cs[indxcn[i]].
grad) + 1,
5258 if (feasb && neqn != 0)
5259 resign(nparam, neqn, &psf, grdpsf, penp, cs, signeq, 12, 12);
5261 for (i = 1; i <=
nob; i++)
5268 gradcn(nparam, indxob[i], (param->
x) + 1, (ob[i].
grad) + 1, constr,
5275 obj(nparam, i, (param->
x) + 1, &(ob[i].
val), param->
cd);
5278 gradob(nparam, i, (param->
x) + 1, (ob[i].
grad) + 1, obj, param->
cd);
5281 fmax =
DMAX1(fmax, -ob[i].val);
5283 fmax =
DMAX1(fmax, ob[i].val);
5285 if (feasb && nob == 0)
5293 for (i = 1; i <=
nob; i++)
5300 sbout2(
glob_prnt.
io, nparam, i,
"gradf(j,",
")", tempv);
5306 dummy, tempv, 2, 2);
5311 sbout2(
glob_prnt.
io, nparam, indxob[i],
"gradg(j,",
")", tempv);
5315 for (i = 1; i <=
ncnst1; i++)
5317 tempv = cs[indxcn[
i]].
grad;
5318 sbout2(
glob_prnt.
io, nparam, indxcn[i],
"gradg(j,",
")", tempv);
5328 for (i = 1; i <=
nparam; i++)
5330 tempv = colvec(hess, i, nparam);
5331 sbout2(
glob_prnt.
io, nparam, i,
"hess (j,",
")", tempv);
5343 out(miter, nparam, nob, nobL, nfsip, nineqn, nn, nineqn, ncnst1,
5344 ncsipl, ncsipn, mesh_pts, param->
x, cs, ob, fM, fmax, steps,
5345 sktnom, d0nm, feasb);
5350 dealloc1(nparam, nrowa, a, hess, hess1, di, d, gm,
5351 grdpsf, penp, bl, bu, clamda, cvec, psmu, span, backup,
5352 iact, iskip, istore);
5355 for (i = 1; i <=
ncnst1; i++)
5356 cs[i].val = backup[i];
5357 for (i = 1; i <=
nob; i++)
5358 ob[i].val = backup[i+ncnst1];
5359 for (i = 1; i <=
neqn; i++)
5361 dealloc1(nparam, nrowa, a, hess, hess1, di, d, gm,
5362 grdpsf, penp, bl, bu, clamda, cvec, psmu, span, backup,
5363 iact, iskip, istore);
5369 if ((ncsipl + ncsipn) != 0 || nfsip)
5370 update_omega(nparam, ncsipl, ncsipn, mesh_pts, nineqn, nob, nobL,
5371 nfsip, steps, fmax, cs, ob, param->
x, viol,
5372 constr, obj, gradob, gradcn, param->
cd, feasb);
5374 dir(nparam, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn,
5375 ncnst1, feasb, &steps, epskt, epseqn, &sktnom, &scvneq, Ck, &d0nm,
5376 &grdftd, indxob, indxcn, iact, &iskp, iskip, istore, param, di, d,
5377 cs, ob, &fM, &fMp, &fmax, &psf, grdpsf, penp, a, bl, bu, clamda, cvec,
5378 hess, hess1, backup, signeq, viol, obj, constr);
5385 step1(nparam, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn,
5386 ncnst1, &ncg, &ncf, indxob, indxcn, iact, &iskp, iskip, istore,
5387 feasb, grdftd, ob, &fM, &fMp, &fmax, &psf, penp, &steps, &scvneq,
5388 bu, param->
x, di, d, cs, backup, signeq, viol, obj, constr,
5394 hessian(nparam, nob, nfsip, nobL, nineqn, neq, neqn, nn, ncsipn, ncnst1,
5395 nfs, &nstart, feasb, bu, param, ob, fmax, &fM, &fMp, &psf, grdpsf,
5396 penp, cs, gm, indxob, indxcn, bl, clamda, di, hess, d, steps, &nrst,
5397 signeq, span, obj, constr, gradob, gradcn, hess1, cvec, psmu, viol);
5414 dealloc1(
int nparam,
int nrowa,
double **a,
double **hess,
double **hess1,
5415 double *di,
double *d,
double *gm,
double *grdpsf,
double *penp,
5416 double *bl,
double *bu,
double *clamda,
double *cvec,
double *psmu,
5417 double *span,
double *backup,
int *iact,
int *iskip,
int *istore)
5420 dealloc1(nparam, nrowa, a, hess, hess1, di, d, gm, grdpsf, penp, bl, bu, clamda,
5421 cvec, psmu, span, backup, iact, iskip, istore)
5424 double *
di, *
d, *
gm, *
grdpsf, *
penp, *
bl, *
bu, *
clamda, *
cvec, *
psmu, *
span,
5455 check(
int nparam,
int nf,
int nfsip,
int *Linfty,
int nineq,
5456 int nnl,
int neq,
int neqn,
int ncsipl,
int ncsipn,
int mode,
5457 int *modem,
double eps,
double bigbnd,
struct _parameter *param)
5460 check(nparam, nf, nfsip, Linfty, nineq, nnl, neq, neqn, ncsipl, ncsipn,
5461 mode, modem, eps, bigbnd, param)
5462 int nparam, nf, nfsip, nineq, nnl, neq, neqn, ncsipl, ncsipn, mode, *modem,
5472 error(
"nparam should be positive! ",
5475 error(
"nf should not be negative! ",
5478 error(
"nineq should not be negative! ",
5480 if (nineq >= 0 && nnl >= 0 && nineq < nnl)
5481 error(
"nineq should be no smaller then nnl! ",
5484 error(
"neqn should not be negative! ",
5487 error(
"neq should not be smaller then neqn ",
5490 error(
"nf should not be smaller then nfsip ",
5492 if (nineq < ncsipn + ncsipl)
5493 error(
"ncsrl+ncsrn should not be larger then nineq",
5495 if (nparam <= neq - neqn)
5496 error(
"Must have nparam > number of linear equalities",
5499 error(
"iprint mod 10 should be 0,1,2 or 3! ",
5503 error(
"eps should be bigger than epsmac! ",
5506 "epsmac = %22.14e which is machine dependent.\n",
5509 if (!(mode == 100 || mode == 101 || mode == 110 || mode == 111
5510 || mode == 200 || mode == 201 || mode == 210 || mode == 211))
5511 error(
"mode is not properly specified! ",
5516 "Error: Input parameters are not consistent.\n");
5519 for (i = 1; i <=
nparam; i++)
5526 "lower bounds should be smaller than upper bounds\n");
5531 if (bli < (-bigbnd))
5566 initpt(
int nparam,
int nnl,
int neq,
int neqn,
int nclin,
int nctotl,
5568 void(*constr)(
int,
int,
double *,
double *,
void *),
5569 void(*gradcn)(
int,
int,
double *,
double *,
5570 void(*)(
int,
int,
double *,
double *,
void *),
void *))
5573 initpt(nparam, nnl, neq, neqn, nclin, nctotl, nrowa, param, cs,
5578 void(* constr)(), (* gradcn)();
5586 hess = make_dm(nparam, nparam);
5587 a = make_dm(nrowa, nparam);
5588 x = make_dv(nparam);
5589 bl = make_dv(nctotl);
5590 bu = make_dv(nctotl);
5591 cvec = make_dv(nparam);
5592 clamda = make_dv(nctotl + nparam + 1);
5593 bj = make_dv(nclin);
5596 for (i = 1; i <=
nclin; i++)
5601 gradcn(nparam, j, (param->
x) + 1, cs[i].
grad + 1, constr, param->
cd);
5603 gradcn(nparam, j + neqn, (param->
x) + 1, cs[i].
grad + 1, constr,
5606 for (i = 1; i <=
nparam; i++)
5613 for (i = nclin; i >= 1; i--)
5614 bj[nclin-i+1] = -cs[i].val;
5615 for (i = nclin; i >= 1; i--)
5616 for (j = 1; j <=
nparam; j++)
5617 a[nclin-i+1][j] = -cs[i].grad[j];
5623 mnn = nrowa + 2 *
nparam;
5626 htemp = convert(hess, nparam, nparam);
5627 atemp = convert(a, nrowa, nparam);
5629 ql0001_(&nclin, &temp1, &nrowa, &nparam, &nparam, &mnn, (htemp + 1),
5630 (cvec + 1), (atemp + 1), (bj + 1), (bl + 1), (bu + 1), (x + 1), (clamda + 1),
5631 &iout, &infoql, &zero, (w + 1), &lenw, (iw + 1), &leniw,
5638 for (i = 1; i <=
nparam; i++)
5639 param->
x[i] = param->
x[i] + x[i];
5645 constr(nparam, j, (param->
x) + 1,
5646 &(cs[i].
val), param->
cd);
5648 constr(nparam, j + neqn, (param->
x) + 1, &(cs[i].
val), param->
cd);
5655 "\n Error: No feasible point is found for the");
5676 update_omega(
int nparam,
int ncsipl,
int ncsipn,
int *mesh_pts,
5677 int nineqn,
int nob,
int nobL,
int nfsip,
double steps,
5680 void(*constr)(
int,
int,
double *,
double *,
void *),
5681 void(*obj)(
int,
int,
double *,
double *,
void *),
5682 void(*gradob)(
int,
int,
double *,
double *,
5683 void(*)(
int,
int,
double *,
double *,
void *),
void *),
5684 void(*gradcn)(
int,
int,
double *,
double *,
5685 void(*)(
int,
int,
double *,
double *,
void *),
void *),
5686 void *cd,
int feasb)
5689 update_omega(nparam, ncsipl, ncsipn, mesh_pts, nineqn, nob, nobL, nfsip,
5690 steps, fmax, cs, ob, x, viol, constr, obj, gradob, gradcn, cd, feasb)
5714 for (i = 1; i <=
ncsipl; i++)
5715 cs[nineq-ncsipl+i].act_sip =
FALSE;
5716 for (i = 1; i <=
ncsipn; i++)
5717 cs[nineqn-ncsipn+i].act_sip =
FALSE;
5719 for (i = nob - nfsip + 1; i <=
nob; i++)
5720 ob[i].act_sip =
FALSE;
5728 offset = nineqn -
ncsipn;
5736 if (cs[offset].val >= -epsilon &&
5737 cs[offset].val >= cs[offset+1].val)
5744 gradcn(nparam, offset, x + 1, cs[offset].grad + 1, constr,
5752 if (cs[offset].val >= -epsilon &&
5753 cs[offset].val > cs[offset-1].val)
5760 gradcn(nparam, offset, x + 1, cs[offset].grad + 1, constr,
5768 if (cs[offset].val >= -epsilon && cs[offset].val >
5769 cs[offset-1].val && cs[offset].val >=
5777 gradcn(nparam, offset, x + 1, cs[offset].grad + 1, constr,
5790 gradcn(nparam, offset, x + 1, cs[offset].grad + 1, constr, cd);
5794 if (cs[offset].
mult > 0.
e0)
5800 if (cs[offset].d1bind)
5807 gradcn(nparam, offset, x + 1, cs[offset].grad + 1, constr, cd);
5824 for (j = 1; j <= mesh_pts[
index]; j++)
5829 if (cs[offset].val >= -epsilon &&
5830 cs[offset].val >= cs[offset+1].val)
5838 if (j == mesh_pts[index])
5840 if (cs[offset].val >= -epsilon &&
5841 cs[offset].val > cs[offset-1].val)
5850 if (cs[offset].val >= -epsilon && cs[offset].val >
5851 cs[offset-1].val && cs[offset].val >=
5875 offset = nineqn -
ncsipn;
5879 g_max = cs[i_max].
val;
5880 if (!cs[i_max].act_sip)
5888 if (cs[offset].val > g_max)
5891 g_max = cs[i_max].
val;
5894 if (!cs[i_max].act_sip)
5899 if (!cs[offset].act_sip)
5910 g_max = cs[i_max].
val;
5911 if (!cs[i_max].act_sip)
5920 for (j = 2;j <= mesh_pts[
index]; j++)
5923 if (cs[offset].val > g_max)
5926 g_max = cs[i_max].
val;
5929 if (!cs[i_max].act_sip)
5934 if (!cs[offset].act_sip)
5956 for (i = 1; i <=
ncsipl; i++)
5957 cs[nineq-ncsipl+i].d1bind =
FALSE;
5958 for (i = 1; i <=
ncsipn; i++)
5959 cs[nineqn-ncsipn+i].d1bind =
FALSE;
5967 offset = nob -
nfsip;
5972 for (i = 1; i <=
index; i++)
5974 for (j = 1; j <= mesh_pts[
i]; j++)
5979 fnow = fabs(ob[offset].val);
5981 fabs(ob[offset].mult_L));
5985 fnow = ob[offset].
val;
5986 fmult = ob[offset].
mult;
5991 fnext = fabs(ob[offset+1].val);
5993 fnext = ob[offset+1].
val;
5994 if ((fnow >= fmax - epsilon) && fnow >= fnext)
6002 gradob(nparam, offset, x + 1,
6003 ob[offset].grad + 1, obj, cd);
6005 gradcn(nparam, offset, x + 1, ob[offset].grad + 1,
6011 else if (j == mesh_pts[i])
6014 fprev = fabs(ob[offset-1].val);
6016 fprev = ob[offset-1].
val;
6017 if ((fnow >= fmax - epsilon) && fnow > fprev)
6025 gradob(nparam, offset, x + 1,
6026 ob[offset].grad + 1, obj, cd);
6028 gradcn(nparam, offset, x + 1, ob[offset].grad + 1,
6038 fprev = fabs(ob[offset-1].val);
6039 fnext = fabs(ob[offset+1].val);
6043 fprev = ob[offset-1].
val;
6044 fnext = ob[offset+1].
val;
6046 if ((fnow >= fmax - epsilon) && fnow > fprev &&
6055 gradob(nparam, offset, x + 1,
6056 ob[offset].grad + 1, obj, cd);
6058 gradcn(nparam, offset, x + 1, ob[offset].grad + 1,
6072 gradob(nparam, offset, x + 1,
6073 ob[offset].grad + 1, obj, cd);
6075 gradcn(nparam, offset, x + 1, ob[offset].grad + 1,
6080 if (fmult != 0.
e0 && !ob[offset].act_sip)
6092 offset = nob -
nfsip;
6097 for (i = 1; i <=
index; i++)
6101 g_max = ob[i_max].
val;
6103 g_max = fabs(ob[i_max].val);
6104 if (!ob[i_max].act_sip)
6109 for (j = 2;j <= mesh_pts[
i];j++)
6113 fnow = ob[offset].
val;
6115 fnow = fabs(ob[offset].val);
6122 if (!ob[i_max].act_sip)
6127 if (!ob[offset].act_sip)
6160 dqp(
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
6161 int,
int,
int *,
struct _parameter *,
double *,
int,
6163 double **,
double *,
double *,
double *,
double *,
6164 double **,
double **,
double *,
double,
int);
6166 di1(
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int,
int *,
6169 double *,
double *,
double *,
double **,
double *,
double);
6177 dir(
int nparam,
int nob,
int nobL,
int nfsip,
int nineqn,
int neq,
int neqn,
6178 int nn,
int ncsipl,
int ncsipn,
int ncnstr,
6179 int feasb,
double *steps,
double epskt,
double epseqn,
6180 double *sktnom,
double *scvneq,
double Ck,
double *d0nm,
6181 double *grdftd,
int *indxob,
int *indxcn,
int *iact,
int *iskp,
6182 int *iskip,
int *istore,
struct _parameter *param,
double *di,
6184 double *fM,
double *fMp,
double *fmax,
double *psf,
double *grdpsf,
6185 double *penp,
double **a,
double *bl,
double *bu,
double *clamda,
6186 double *cvec,
double **hess,
double **hess1,
6187 double *backup,
double *signeq,
struct _violation *viol,
6188 void(*obj)(
int,
int,
double *,
double *,
void *),
6189 void(*constr)(
int,
int,
double *,
double *,
void *))
6192 dir(nparam, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn, ncnstr,
6193 feasb, steps, epskt, epseqn, sktnom, scvneq, Ck, d0nm,
6194 grdftd, indxob, indxcn, iact, iskp, iskip, istore, param, di, d, cs, ob,
6195 fM, fMp, fmax, psf, grdpsf, penp, a, bl, bu, clamda, cvec, hess, hess1,
6196 backup, signeq, viol, obj, constr)
6197 int nparam, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn, ncnstr,
6199 int *indxob, *indxcn, *iact, *iskip, *istore;
6200 double *steps, epskt, epseqn, *sktnom, Ck, *d0nm, *grdftd, *fM, *fMp,
6201 *fmax, *psf, *scvneq;
6202 double *di, *d, *grdpsf, *penp, **a, *bl, *bu, *clamda, *cvec, **hess,
6203 **hess1, *backup, *signeq;
6208 void(* obj)(), (* constr)();
6211 int i,
j,
k, kk,
ncg,
ncf, nqprm0,
nclin0,
nctot0,
infoqp, nqprm1,
ncl,
6214 double fmxl,
vv,
dx,
dmx,
dnm1,
dnm,
v0,
v1,
vk=0.,
temp1,
temp2,
theta,
6218 ncg = ncf = *iskp = 0;
6223 adummy = make_dv(1);
6226 temp1 = temp2 = 0.e0;
6240 nqprm0 = nparam + 1;
6246 for (i = 1; i <=
ncnstr; i++)
6257 iskip[ncl+2-
i] = nineqn +
i;
6261 iw[
i] = nineqn + neqn +
i;
6264 for (i = 1; i <=
nob; i++)
6268 dqp(nparam, nqprm0, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn,
6269 ncnstr, nctot0, nrowa0, nineqn, &infoqp, param, di, feasb, ob,
6270 *fmax, grdpsf, cs, a, cvec, bl, bu, clamda, hess, hess1, di, vv, 0);
6287 for (i = nn; i >= 1; i--)
6289 if (fuscmp(cs[indxcn[i]].
mult, thrshd))
6291 iact[
j] = indxcn[
i];
6296 iact[
k] = indxcn[
i];
6305 for (i = nob; i >= 1; i--)
6308 ltem1 = fuscmp(ob[i].
mult, thrshd);
6309 ltem2 = (nobL !=
nob) && (fuscmp(ob[i].mult_L, thrshd));
6323 vv = ob[iact[nn+1]].
val;
6324 *d0nm =
sqrt(scaprd(nparam, di, di));
6327 dx =
sqrt(scaprd(nparam, param->
x, param->
x));
6331 for (i = 1; i <=
nparam; i++)
6332 di[i] = di[i] * dmx / (*d0nm);
6336 matrvc(nparam, nparam, hess, di, w);
6338 *grdftd = -scaprd(nparam, w, di);
6339 *sktnom =
sqrt(scaprd(nparam, w, w));
6340 if (((*d0nm <= epskt) || ((gLgeps > 0.
e0) && (*sktnom <= gLgeps))) &&
6341 (neqn == 0 || *scvneq <= epseqn))
6366 " for g ", cs[1].
mult);
6367 for (j = 2; j <=
ncnstr; j++)
6373 " for f ", ob[1].
mult);
6374 for (j = 2; j <=
nob; j++)
6386 && *scvneq > epseqn)
6398 *grdftd = slope(nob, nobL, neqn, nparam, feasb, ob, grdpsf,
6399 di, d, *fmax, dummy, 0, adummy, 0);
6401 if (nn == 0 && nobL <= 1)
6403 for (i = 1; i <=
nparam; i++)
6423 vk =
DMIN1(Ck * (*d0nm) * (*d0nm), *d0nm);
6425 for (i = 1; i <=
nn; i++)
6427 tempv = cs[indxcn[
i]].
grad;
6428 grdgd0 = scaprd(nparam, tempv, di);
6439 nqprm1 = nparam + 1;
6441 nclin1 = ncnstr +
DMAX1(nobL, 1);
6444 nrowa1 =
DMAX1(nclin1, 1);
6446 di1(nparam, nqprm1, nob, nobL, nfsip, nineqn, neq, neqn, ncnstr,
6448 param, di, ob, *fmax, grdpsf, cs, cvec, bl, bu, clamda,
6459 dnm1 =
sqrt(scaprd(nparam, d, d));
6469 for (i = 1; i <=
nparam; i++)
6474 v0 = pow(*d0nm, 2.1);
6475 v1 =
DMAX1(0.5
e0, pow(dnm1, 2.5));
6476 rho = v0 / (v0 +
v1);
6484 for (i = 1; i <=
nn; i++)
6486 tempv = cs[indxcn[
i]].
grad;
6487 grdgd0 = scaprd(nparam, tempv, di);
6488 grdgd1 = scaprd(nparam, tempv, d);
6499 rhol =
DMAX1(rhol, -temp1 / temp2);
6500 if (temp2 < 0.
e0 && rhol < 1.
e0)
6516 rhog = slope(nob, nobL, neqn, nparam, feasb, ob, grdpsf,
6518 rhog =
DMIN1(rhol, rhog);
6524 grdfd1 = scaprd(nparam, ob[1].grad, d);
6527 grdfd1 = grdfd1 - scaprd(nparam, grdpsf, d);
6529 temp2 = (theta - 1.e0) * grdfd0 / temp1;
6533 rhog =
DMIN1(rhol, temp2);
6536 if (*steps == 1.
e0 && rhol < 0.5
e0)
6539 for (i = 1; i <=
nparam; i++)
6543 di[
i] = (1.e0 -
rho) * di[i] + rho * d[i];
6545 dnm =
sqrt(scaprd(nparam, di, di));
6568 step1(nparam, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn,
6569 ncnstr, &ncg, &ncf, indxob, indxcn, iact, iskp, iskip, istore,
6570 feasb, *grdftd, ob, fM, fMp, fmax, psf, penp, steps, scvneq, bu,
6571 param->
x, di, d, cs, backup, signeq, viol, obj, constr, param->
cd);
6580 if (rho != rhog && nn != 0)
6581 for (i = 1; i <=
nparam; i++)
6582 di[i] = (1 - rhog) * cvec[
i] + rhog * d[
i];
6583 dnm =
sqrt(scaprd(nparam, di, di));
6594 *grdftd = slope(nob, nobL, neqn, nparam, feasb, ob,
6595 grdpsf, di, d, *fmax, theta, 0, bl, 1);
6597 for (i = 1; i <=
nparam; i++)
6598 bu[i] = param->
x[i] + di[i];
6620 iskip[ncl+1-
j] = kk;
6625 tempv = cs[kk].
grad;
6626 temp1 = dnm *
sqrt(scaprd(nparam, tempv, tempv));
6627 temp2 = cs[kk].
mult;
6629 if (temp2 != 0.
e0 || cs[kk].val >= (-0.2
e0*temp1) ||
6634 if (feasb && kk <= nineqn)
6636 constr(nparam, kk, bu + 1, &(cs[kk].val), param->
cd);
6641 fmxl =
DMAX1(fmxl, cs[kk].val);
6645 if (fabs(fmxl) > bgbnd)
6647 for (i = 1; i <=
nparam; i++)
6662 if ((neqn != 0) && (feasb))
6663 resign(nparam, neqn, psf, grdpsf, penp, cs, signeq, 10, 20);
6667 for (i = 1; i <=
ncg; i++)
6670 istore[iact[
i]] = 1;
6671 fmxl =
DMAX1(fmxl, cs[iact[i]].val);
6672 if (fabs(fmxl) > bgbnd)
6674 for (i = 1; i <=
nparam; i++)
6694 if (ob[iact[nn+1]].
mult < 0.
e0)
6696 for (i = nff; i <=
nob; i++)
6705 for (j = 1; j <=
nparam; j++)
6706 w[nparam+j] = sign * ob[iact[k]].grad[j] - ob[kk].grad[j];
6707 temp1 = fabs(ob[kk].val - sign * vv);
6708 temp2 = dnm *
sqrt(scaprd(nparam, &w[nparam], &w[nparam]));
6709 if (temp1 != 0.
e0 && temp2 != 0.
e0)
6711 temp1 = temp1 /
temp2;
6712 temp2 = ob[kk].
mult;
6713 if (temp2 == 0.
e0 && temp1 > 0.2
e0)
6717 iw[nobb+ninq+
neq] = kk;
6719 istore[nineqn+kk] = 1;
6724 constr(nparam, indxob[kk], bu + 1, &(ob[kk].val), param->
cd);
6729 obj(nparam, kk, bu + 1, &(ob[kk].val), param->
cd);
6732 fmxl =
DMAX1(fmxl, -ob[kk].val);
6734 fmxl =
DMAX1(fmxl, ob[kk].val);
6735 if (fabs(fmxl) > bgbnd)
6737 for (i = 1; i <=
nparam; i++)
6747 for (i = 1; i <=
ncf; i++)
6749 iw[ninq+neq+
i] = iact[i+
nn];
6750 istore[nineqn+iact[i+
nn]] = 1;
6751 fmxl =
DMAX1(fmxl, ob[iact[i+nn]].val);
6753 fmxl =
DMAX1(fmxl, -ob[iact[i+nn]].val);
6754 if (fabs(fmxl) > bgbnd)
6756 for (i = 1; i <=
nparam; i++)
6766 matrvc(nparam, nparam, hess, di, cvec);
6767 vv = -
DMIN1(0.01
e0 * dnm, pow(dnm, 2.5));
6778 nclin0 = ninq +
neq;
6782 nqprm0 = nparam + 1;
6783 nclin0 = ninq + neq + nobbL;
6785 nctot0 = nqprm0 +
nclin0;
6786 nrowa0 =
DMAX1(nclin0, 1);
6789 dqp(nparam, nqprm0, nobb, nobbL, nfsip, nncn, neq, neqn, nn, ncsipl, ncsipn, i,
6790 nctot0, nrowa0, nineqn, &infoqp, param, di, feasb, ob, fmxl,
6791 grdpsf, cs, a, cvec, bl, bu, clamda, hess, hess1, d, vv, 1);
6792 dnmtil =
sqrt(scaprd(nparam, d, d));
6793 if (infoqp != 0 || dnmtil > dnm)
6795 for (i = 1; i <=
nparam; i++)
6803 for (i = 1; i <= nineqn +
nob; i++)
6822 dqp(
int nparam,
int nqpram,
int nob,
int nobL,
int nfsip,
int nineqn,
6823 int neq,
int neqn,
int nn,
int ncsipl,
int ncsipn,
int ncnstr,
6824 int nctotl,
int nrowa,
int nineqn_tot,
int *infoqp,
6826 double fmax,
double *grdpsf,
struct _constraint *cs,
6827 double **a,
double *cvec,
double *bl,
double *bu,
double *clamda,
6828 double **hess,
double **hess1,
double *x,
6832 dqp(nparam, nqpram, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn,
6833 ncnstr, nctotl, nrowa, nineqn_tot, infoqp, param, di, feasb, ob,
6834 fmax, grdpsf, cs, a, cvec, bl, bu, clamda, hess, hess1, x, vv, job)
6835 int nparam,
nqpram,
nob,
nobL,
nfsip,
nineqn,
neq,
neqn,
nn,
ncsipl,
ncsipn,
6844 int i, ii,
j, jj, ij,
k,
iout,
mnn, nqnp,
zero,
temp1,
temp2,
ncnstr_used,
6850 bj = make_dv(nrowa);
6851 iw_hold = make_iv(nrowa);
6852 for (i = 1; i <=
nparam; i++)
6859 bl[
i] = param->
bl[
i] - x0i - xdi;
6860 bu[
i] = param->
bu[
i] - x0i - xdi;
6861 cvec[
i] = cvec[
i] - grdpsf[
i];
6873 for (i = 1; i <=
ncnstr; i++)
6875 jj = iw[ncnstr+1-
i];
6882 if (i <= (neq - neqn) || (i > neq && i <= (ncnstr - nineqn)))
6886 bj[
k] = x0i - cs[jj].
val;
6887 for (j = 1; j <=
nparam; j++)
6888 a[k][j] = -cs[jj].grad[j];
6900 for (i = 1; i <=
nparam; i++)
6901 cvec[i] = cvec[i] + ob[1].grad[i];
6909 for (i = 1; i <=
nob; i++)
6910 if (ob[iw[ncnstr+i]].act_sip)
6913 for (i = 1; i <=
nob; i++)
6916 if ((i <= nob - nfsip) || ob[iw[ij]].act_sip)
6919 iw_hold[
k] = iw[ij];
6920 bj[
k] = fmax - ob[iw[ij]].
val;
6922 bj[k+numf_used] = fmax + ob[iw[ij]].
val;
6923 for (j = 1; j <=
nparam; j++)
6925 a[
k][
j] = -ob[iw[ij]].
grad[
j];
6927 a[k+numf_used][
j] = ob[iw[ij]].
grad[
j];
6931 a[k+numf_used][
nqpram] = 1.e0;
6938 matrcp(nparam, hess, nparam + 1, hess1);
6946 htemp = convert(hess1, nparam + 1, nparam + 1);
6947 atemp = convert(a, nrowa, nqpram);
6949 ql0001_(&k, &temp1, &nrowa, &nqpram, &temp2, &mnn, (htemp + 1),
6950 (cvec + 1), (atemp + 1), (bj + 1), (bl + 1), (bu + 1), (x + 1),
6951 (clamda + 1), &iout, infoqp, &zero, (w + 1), &lenw, (iw + 1), &leniw,
6956 if (*infoqp != 0 || job == 1)
6967 if (ncsipl + ncsipn)
6968 for (i = 1; i <=
ncnstr; i++)
6971 for (i = 1; i <=
nob; i++)
6976 for (i = 1; i <=
nqpram; i++)
6979 if (clamda[ii] == 0.
e0 && clamda[ii+nqpram] == 0.
e0)
6981 else if (clamda[ii] != 0.
e0)
6982 clamda[ii] = -clamda[ii];
6984 clamda[ii] = clamda[ii+
nqpram];
6987 for (i = 1; i <=
nqpram; i++)
6991 for (i = 1; i <= numf_used; i++)
6993 ij = ncnstr_used +
i;
6996 ii = k - 2 * numf_used +
i;
6997 ob[iw_hold[ij]].
mult = clamda[ii] - clamda[ii+numf_used];
6998 ob[iw_hold[ij]].
mult_L = clamda[ii+numf_used];
7002 ii = k - numf_used +
i;
7003 ob[iw_hold[ij]].
mult = clamda[ii];
7008 cs[iw_hold[i]].
mult = clamda[i];
7019 di1(
int nparam,
int nqpram,
int nob,
int nobL,
int nfsip,
int nineqn,
7020 int neq,
int neqn,
int ncnstr,
int ncsipl,
int ncsipn,
7021 int nrowa,
int *infoqp,
int mode,
struct _parameter *param,
7023 *grdpsf,
struct _constraint *cs,
double *cvec,
double *bl,
double *bu,
7024 double *clamda,
double **hess1,
double *x,
double steps)
7027 di1(nparam, nqpram, nob, nobL, nfsip, nineqn, neq, neqn, ncnstr, ncsipl,
7028 ncsipn, nrowa, infoqp, mode, param, d0, ob, fmax, grdpsf, cs,
7029 cvec, bl, bu, clamda, hess1, x, steps)
7030 int nparam, nqpram, nob, nobL, nfsip, nineqn, neq, neqn, ncnstr,
7031 nrowa, *infoqp, mode, ncsipl, ncsipn;
7033 double *d0, *grdpsf, *cvec, *bl, *bu, *clamda, **hess1, *x;
7039 int i,
k, ii, jj,
iout,
j,
mnn,
zero,
temp1,
temp3,
ncnstr_used, numf_used=0;
7043 if ((ncsipl + ncsipn) != 0)
7052 nrowa =
DMAX1(nrowa, 1);
7053 a = make_dm(nrowa, nqpram);
7054 bj = make_dv(nrowa);
7055 iw_hold = make_iv(nrowa);
7061 for (i = 1; i <=
nparam; i++)
7067 cvec[
i] = -eta * d0[
i];
7076 for (i = 1; i <=
ncnstr; i++)
7078 jj = ncnstr + 1 -
i;
7084 bj[
k] = -cs[jj].
val;
7085 for (j = 1; j <=
nparam; j++)
7086 a[k][j] = -cs[jj].grad[j];
7088 if ((i > (neq - neqn) && i <= neq) || i > ii)
7098 for (i = 1; i <=
nob; i++)
7100 if ((i <= nob - nfsip) || ob[i].act_sip)
7103 bj[
k] = fmax - ob[
i].
val;
7104 for (j = 1; j <=
nparam; j++)
7108 a[k+numf_used][
j] = ob[
i].
grad[
j] + grdpsf[
j];
7112 a[k+numf_used][
nqpram] = 1.e0;
7119 for (j = 1; j <=
nparam; j++)
7120 a[k][j] = grdpsf[j];
7124 diagnl(nqpram, eta, hess1);
7132 temp3 = k + numf_used;
7135 mnn = temp3 + 2 *
nqpram;
7136 htemp = convert(hess1, nparam + 1, nparam + 1);
7137 atemp = convert(a, nrowa, nqpram);
7139 ql0001_(&temp3, &temp1, &nrowa, &nqpram, &nqpram, &mnn, (htemp + 1),
7140 (cvec + 1), (atemp + 1), (bj + 1), (bl + 1), (bu + 1), (x + 1),
7141 (clamda + 1), &iout, infoqp, &zero, (w + 1), &lenw, (iw + 1), &leniw,
7149 if (ncsipl + ncsipn)
7152 if (clamda[i] > 0.
e0)
7167 step1(
int nparam,
int nob,
int nobL,
int nfsip,
int nineqn,
int neq,
int neqn,
7168 int nn,
int ncsipl,
int ncsipn,
int ncnstr,
int *ncg,
int *ncf,
7169 int *indxob,
int *indxcn,
int *iact,
int *iskp,
int *iskip,
7170 int *istore,
int feasb,
double grdftd,
struct _objective *ob,
7171 double *fM,
double *fMp,
double *fmax,
double *psf,
double *penp,
7172 double *steps,
double *scvneq,
double *
xnew,
double *x,
double *di,
7173 double *d,
struct _constraint *cs,
double *backup,
double *signeq,
7175 void(*obj)(
int,
int,
double *,
double *,
void *),
7176 void(*constr)(
int,
int,
double *,
double *,
void *),
void *cd)
7179 step1(nparam, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn, ncnstr,
7180 ncg, ncf, indxob, indxcn, iact, iskp, iskip, istore, feasb, grdftd, ob,
7181 fM, fMp, fmax, psf, penp, steps, scvneq, xnew, x, di, d, cs, backup,
7182 signeq, sip_viol, obj, constr, cd)
7183 int nparam, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn, ncnstr,
7184 *ncg, *ncf, feasb, *iskp;
7185 int *indxob, *indxcn, *iact, *iskip, *istore;
7186 double grdftd, *fM, *fMp, *fmax, *steps, *scvneq, *psf;
7187 double *xnew, *x, *di, *d, *penp, *backup, *signeq;
7191 void(* obj)(), (* constr)();
7195 int i, ii, ij, jj,
itry,
ikeep,
j, job,
nlin, mnm, ltem1, ltem2, reform,
7201 ostep = *steps = 1.e0;
7204 sipldone = (ncsipl == 0);
7208 prod1 = (0.1e0) * grdftd;
7210 adummy = make_dm(1, 1);
7211 adummy[1][1] = 0.e0;
7226 prod = prod1 * (*steps);
7227 if (!feasb || (nobL > 1))
7228 prod = prod +
tolfe;
7229 for (i = 1; i <=
nparam; i++)
7232 xnew[
i] = x[
i] + (*steps) * di[
i];
7234 xnew[
i] = x[
i] + (*steps) * di[
i] + d[
i] * (*steps) * (*steps);
7240 *(adummy + 1), 1, 2);
7250 for (i = ii; i <= *
iskp; i++)
7253 constr(nparam, ij, xnew + 1, &(cs[ij].val), cd);
7258 "\t\t\t trial constraints %d \t %22.14e\n", ij,
7262 "\t\t\t\t\t %d \t %22.14e\n", ij, cs[ij].
val);
7264 if (cs[ij].val <= tolfe)
7270 sip_viol->
index = ij;
7275 sip_viol->
index = 0;
7286 for (i = jj; i <=
ncsipl; i++)
7289 if (cs[ij].act_sip || element(iskip, nlin - ikeep, ij))
7291 constr(nparam, ij, xnew + 1, &(cs[ij].val), cd);
7296 "\t\t\t trial constraints %d \t %22.14e\n", ij,
7300 "\t\t\t\t\t %d \t %22.14e\n", ij, cs[ij].
val);
7302 if (cs[ij].val <= tolfe)
7306 sip_viol->
index = ij;
7319 for (i = 1; i <=
nn; i++)
7330 constr(nparam, ii, xnew + 1, &(cs[ii].val), cd);
7336 if (i == 1 && ikeep == nlin)
7338 "\t\t\t trial constraints %d \t %22.14e\n", ii,
7340 if (i != 1 || ikeep != nlin)
7342 "\t\t\t\t\t %d \t %22.14e\n", ii, cs[ii].
val);
7346 shift(nn, ii, iact);
7347 if (ncsipn && ii > nineqn - ncsipn)
7350 sip_viol->
index = ii;
7355 sip_viol->
index = 0;
7361 if (ncsipn && ii > nineqn - ncsipn)
7364 sip_viol->
index = ii;
7369 sip_viol->
index = 0;
7375 cdone = eqdone =
TRUE;
7385 for (i = 0; i <=
nob; i++)
7387 if (nob != 0 && i == 0)
7393 if (!(eqdone || neqn == 0))
7395 for (j = 1; j <=
neqn; j++)
7406 resign(nparam, neqn, psf, *(adummy + 1), penp, cs, signeq,
7409 if (istore[nineqn+ii] != 1 && i != 0)
7411 obj(nparam, ii, xnew + 1, &(ob[ii].val), cd);
7420 "\t\t\t trial penalty term \t %22.14e\n", -*psf);
7423 "\t\t\t trial objectives %d \t %22.14e\n",
7427 "\t\t\t\t\t %d \t %22.14e\n", ii, fii - *psf);
7431 if (istore[ii] != 1)
7433 constr(nparam, indxob[ii], xnew + 1, &(ob[ii].val), cd);
7436 if (ob[ii].val > tolfe)
7440 "\t\t\t trial objectives %d \t %22.14e\n",
7441 indxob[ii], ob[ii].
val);
7444 "\t\t\t\t\t %d \t %22.14e\n", indxob[ii],
7448 fmax1 =
DMAX1(fmax1, fii);
7450 fmax1 =
DMAX1(fmax1, -fii);
7451 if (!feasb && reform)
7455 if ((fii - *psf) > (*fMp + prod))
7458 shift(nob, ii, &iact[nn]);
7459 if (nfsip && ii > nob - nfsip)
7462 sip_viol->
index = ii;
7467 sip_viol->
index = 0;
7471 if (nobL == nob || (-fii - *psf) <= (*fMp + prod))
7474 shift(nob, ii, &iact[nn]);
7475 if (nfsip && ii > nob - nfsip)
7478 sip_viol->
index = ii;
7483 sip_viol->
index = 0;
7487 ltem1 = (fii - *
psf) > (*fMp + prod);
7488 ltem2 = (nobL !=
nob) && ((-fii - *psf) > (*fMp +
prod));
7493 fdone = eqdone =
TRUE;
7496 if (ostep == *steps)
7497 mnm = ikeep + neq -
neqn;
7498 if (ostep != *steps)
7500 for (i = 1; i <= mnm; i++)
7503 if (ikeep != nlin && ostep == *steps)
7506 ii = iskip[nlin+2-
i];
7508 ii = indxcn[nn+i-ikeep+
nlin];
7510 constr(nparam, ii, xnew + 1, &(cs[ii].val), cd);
7513 for (i = 1; i <=
ncnstr; i++)
7516 *scvneq = *scvneq - cs[
i].
val;
7517 backup[
i] = cs[
i].
val;
7519 for (i = 1; i <=
nob; i++)
7520 backup[i+ncnstr] = ob[i].val;
7521 if (!feasb && reform)
7523 for (i = 1; i <=
nparam; i++)
7533 *fMp = fmax1 - *
psf;
7535 for (i = 1; i <=
nn; i++)
7536 iact[i] = indxcn[i];
7537 for (i = 1; i <=
nob; i++)
7541 cdone = fdone = eqdone = reform =
FALSE;
7547 for (i = 1; i <= nob +
nineqn; i++)
7549 *steps = *steps * 0.5e0;
7559 for (i = 1; i <= nob +
nineqn; i++)
7571 hessian(
int nparam,
int nob,
int nfsip,
int nobL,
int nineqn,
int neq,
7572 int neqn,
int nn,
int ncsipn,
int ncnstr,
int nfs,
int *nstart,
7573 int feasb,
double *xnew,
struct _parameter *param,
7574 struct _objective *ob,
double fmax,
double *fM,
double *fMp,
7575 double *psf,
double *grdpsf,
double *penp,
struct _constraint *cs,
7576 double *gm,
int *indxob,
int *indxcn,
double *
delta,
double *eta,
7577 double *
gamma,
double **hess,
double *
hd,
double steps,
int *nrst,
7578 double *signeq,
double *span,
7579 void(*obj)(
int,
int,
double *,
double *,
void *),
7580 void(*constr)(
int,
int,
double *,
double *,
void *),
7581 void(*gradob)(
int,
int,
double *,
double *,
7582 void(*)(
int,
int,
double *,
double *,
void *),
void *),
7583 void(*gradcn)(
int,
int,
double *,
double *,
7584 void(*)(
int,
int,
double *,
double *,
void *),
void *),
7585 double **
phess,
double *
psb,
double *psmu,
7589 hessian(nparam, nob, nfsip, nobL, nineqn, neq, neqn, nn, ncsipn, ncnstr,
7590 nfs, nstart, feasb, xnew, param, ob, fmax, fM, fMp, psf, grdpsf, penp,
7591 cs, gm, indxob, indxcn, delta, eta, gamma, hess, hd, steps, nrst, signeq,
7592 span, obj, constr, gradob, gradcn, phess, psb, psmu, sip_viol)
7593 int nparam,
nob,
nobL,
nineqn,
neq,
neqn,
nn,
nfsip,
ncsipn,
ncnstr,
7603 void(* obj)(), (* constr)(), (* gradob)(), (* gradcn)();
7616 if (feasb && nstop && !neqn)
7617 if ((fabs(w[1] - fmax) <=
objeps) ||
7618 (fabs(w[1] - fmax) <= objrep*fabs(w[1])))
7622 for (i = 1; i <=
nparam; i++)
7623 param->
x[i] = xnew[i];
7640 for (j = 1; j <= 2; j++)
7645 for (i = 1; i <=
nparam; i++)
7648 for (k = 1; k <=
nob; k++)
7649 hd[i] = hd[i] + ob[k].grad[i] * ob[k].
mult;
7656 for (i = 1; i <=
nparam; i++)
7659 for (k = 1; k <=
nineqn; k++)
7660 gamma[i] = gamma[i] + cs[k].grad[i] * cs[k].
mult;
7665 for (i = 1; i <=
nparam; i++)
7668 for (k = 1; k <=
neqn; k++)
7674 for (i = 1; i <=
nparam; i++)
7679 psb[
i] = hd[
i] + param->
mult[
i] + gamma[
i];
7680 gamma[
i] = gamma[
i] + hd[
i] - grdpsf[
i] + eta[
i];
7686 gamma[
i] = gamma[
i] + ob[1].
grad[
i] - grdpsf[
i] + eta[
i];
7691 psb[
i] = param->
mult[
i] + gamma[
i];
7692 gamma[
i] = gamma[
i] - grdpsf[
i] + eta[
i];
7695 delta[
i] = gamma[
i];
7701 for (i = 1; i <=
nn; i++)
7703 if ((feasb) && (i > nineqn))
7704 signgj = signeq[i-
nineqn];
7705 if ((!feasb) || (i <= nineqn))
7707 if ((feasb) && (ncsipn) && (i > nineqn - ncsipn) &&
7708 (cs[indxcn[i]].
mult == 0.
e0))
7711 gradcn(nparam, indxcn[i], xnew + 1, cs[indxcn[i]].grad + 1,
7714 resign(nparam, neqn, psf, grdpsf, penp, cs, signeq, 11, 11);
7716 for (i = 1; i <=
nob; i++)
7719 if ((i <= nob - nfsip) || (i > nob - nfsip &&
7720 ((ob[i].
mult != 0.
e0) || (ob[i].mult_L != 0.
e0))))
7723 gradob(nparam, i, xnew + 1, ob[i].grad + 1, obj, param->
cd);
7725 gradcn(nparam, indxob[i], xnew + 1, ob[i].grad + 1,
7736 if (!(feasb && steps < delta_s && ((sip_viol->
type ==
OBJECT &&
7740 if (*nrst < (5*nparam) || steps > 0.1
e0)
7743 for (i = 1; i <=
nparam; i++)
7745 gamma[
i] = gamma[
i] - delta[
i];
7746 delta[
i] = xnew[
i] - param->
x[
i];
7748 matrvc(nparam, nparam, hess, delta, hd);
7749 dhd = scaprd(nparam, delta, hd);
7757 gammd = scaprd(nparam, delta, gamma);
7758 if (gammd >= (0.2
e0*dhd))
7761 theta = 0.8e0 * dhd / (dhd -
gammd);
7762 for (i = 1; i <=
nparam; i++)
7763 eta[i] = hd[i] * (1.
e0 - theta) + theta * gamma[
i];
7764 etad = theta * gammd + (1.e0 -
theta) * dhd;
7765 for (i = 1; i <=
nparam; i++)
7767 for (j = i; j <=
nparam; j++)
7769 hess[
i][
j] = hess[
i][
j] - hd[
i] * hd[
j] / dhd +
7771 hess[
j][
i] = hess[
i][
j];
7781 for (i = 1; i <=
nparam; i++)
7782 param->
x[i] = xnew[i];
7785 if (neqn != 0 && (feasb))
7790 for (j = 1; j <=
nparam; j++)
7793 for (k = 1; k <=
i; k++)
7794 gamma[j] = gamma[j] + cs[k+nineqn].grad[j] *
7797 for (i = 1; i <=
nparam; i++)
7798 psb[i] = psb[i] + gamma[i];
7803 for (j = 1; j <=
nparam; j++)
7806 for (k = 1; k <=
i; k++)
7810 for (i = 1; i <=
nparam; i++)
7811 psb[i] = psb[i] + gamma[i];
7814 estlam(nparam, neqn, &ifail, bgbnd, phess, delta, eta,
7815 gamma, cs, psb, hd, xnew, psmu);
7818 for (i = 1; i <=
neqn; i++)
7821 penp[
i] = 2.e0 * penp[
i];
7824 etad = psmu[
i] + penp[
i];
7826 penp[
i] =
DMAX1((1.
e0 - psmu[i]), (2.
e0 * penp[i]));
7828 if (penp[i] > bgbnd)
7835 resign(nparam, neqn, psf, grdpsf, penp, cs, signeq, 20, 12);
7841 np = indexs(*nstart, nfs);
7843 for (i = 1; i <=
neqn; i++)
7848 for (i = 1; i <=
neqn; i++)
7853 mnm =
DMIN1(*nstart, nfs);
7854 for (i = 2; i <= mnm; i++)
7859 for (j = 1; j <=
neqn; j++)
7860 psfnew = psfnew + gm[(i-1)*neqn +
j]*penp[
j];
7862 *fM =
DMAX1(*fM, span[i]);
7863 *fMp =
DMAX1(*fMp, span[i] - psfnew);
7868 for (i = 1; i <=
nob; i++)
7872 sbout2(
glob_prnt.
io, nparam, indxob[i],
"gradg(j,",
7881 dummy, ob[1].
grad, 2, 2);
7885 for (i = 1; i <=
ncnstr; i++)
7887 tempv = cs[indxcn[
i]].
grad;
7888 sbout2(
glob_prnt.
io, nparam, indxcn[i],
"gradg(j,",
")", tempv);
7893 dummy, grdpsf, 2, 2);
7905 " for g ", cs[1].
mult);
7906 for (j = 2; j <=
ncnstr; j++)
7912 " for f ", ob[1].
mult);
7913 for (j = 2; j <=
nob; j++)
7916 for (i = 1; i <=
nparam; i++)
7918 tempv = colvec(hess, i, nparam);
7919 sbout2(
glob_prnt.
io, nparam, i,
"hess (j,",
")", tempv);
7931 out(
int miter,
int nparam,
int nob,
int nobL,
int nfsip,
int ncn,
7932 int nn,
int nineqn,
int ncnstr,
int ncsipl,
int ncsipn,
7934 struct _objective *ob,
double fM,
double fmax,
7935 double steps,
double sktnom,
double d0norm,
int feasb)
7938 out(miter, nparam, nob, nobL, nfsip, ncn, nn, nineqn, ncnstr, ncsipl, ncsipn,
7939 mesh_pts, x, cs, ob, fM, fmax, steps, sktnom, d0norm, feasb)
7940 int miter, nparam, nob, nobL, nfsip, ncn, nn, ncnstr, feasb,
7941 ncsipl, ncsipn, nineqn, *mesh_pts;
7942 double fM, fmax, steps, sktnom, d0norm;
7951 adummy = make_dv(1);
7973 if (feasb && nob > 0)
7976 for (i = 1; i <= nob -
nfsip; i++)
7985 offset = nob -
nfsip;
7989 gmax = ob[++offset].
val;
7991 gmax = fabs(ob[++offset].val);
7992 for (j = 2; j <= mesh_pts[
i]; j++)
7995 if (nob == nobL && ob[offset].val > gmax)
7996 gmax = ob[offset].
val;
7997 else if (nob != nobL && fabs(ob[offset].val) > gmax)
7998 gmax = fabs(ob[offset].val);
8013 for (i = 1; i <= nineqn -
ncsipn; i++)
8017 offset = nineqn -
ncsipn;
8020 gmax = cs[++offset].
val;
8024 if (cs[offset].val > gmax)
8025 gmax = cs[offset].
val;
8037 gmax = cs[++offset].
val;
8042 for (j = 2; j <= mesh_pts[
index]; j++)
8045 if (cs[offset].val > gmax)
8046 gmax = cs[offset].
val;
8074 display = (nstop == 0);
8085 for (i = 1; i <= nob -
nfsip; i++)
8096 offset = nob -
nfsip;
8101 for (i = 1; i <=
index; i++)
8104 gmax = ob[++offset].
val;
8106 gmax = fabs(ob[++offset].val);
8107 for (j = 2; j <= mesh_pts[
i]; j++)
8110 if (nob == nobL && ob[offset].val > gmax)
8111 gmax = ob[offset].
val;
8112 else if (nob != nobL && fabs(ob[offset].val) > gmax)
8113 gmax = fabs(ob[offset].val);
8121 if (nob > 1 && display)
8123 if (ncnstr != 0 && display)
8126 for (i = 1; i <= nineqn -
ncsipn; i++)
8130 offset = nineqn -
ncsipn;
8133 gmax = cs[++offset].
val;
8137 if (cs[offset].val > gmax)
8138 gmax = cs[offset].
val;
8150 gmax = cs[++offset].
val;
8155 for (j = 2; j <= mesh_pts[
index];
8159 if (cs[offset].val > gmax)
8160 gmax = cs[offset].
val;
8171 SNECV = SNECV + fabs(cs[i].val);
8174 SNECV, adummy, 1, 1);
8201 "\n The following was calculated during iteration %5d:\n",
8215 "\nNormal termination: You have obtained a solution !!\n");
8218 "Warning: Norm of Kuhn-Tucker vector is large !!\n");
8222 "\nWarning: Maximum iterations have been reached before\n");
8228 "\nError : Step size has been smaller than the computed\n");
8233 "\nError: Failure in constructing d0 !!\n\n");
8236 "\nError: Failure in constructing d1 !!\n\n");
8240 "\nError: The new iterate is numerically equivalent to the\n");
8242 "previous iterate, though the stopping criterion is not \n");
8248 "\nError: Could not satisfy nonlinear equality constraints -\n");
8264 void grobfd(
int nparam,
int j,
double *x,
double *
gradf,
8265 void(*obj)(
int,
int,
double *,
double *,
void *),
void *cd)
8267 void grobfd(nparam, j, x, gradf, obj, cd)
8277 for (i = 0; i <= nparam - 1; i++)
8308 void grcnfd(
int nparam,
int j,
double *x,
double *
gradg,
8309 void(*constr)(
int,
int,
double *,
double *,
void *),
void *cd)
8311 void grcnfd(nparam, j, x, gradg, constr, cd)
8321 for (i = 0; i <= nparam - 1; i++)
8359 static void fool(
double,
double,
double *);
8369 static void diagnl(
int nrowa,
double diag,
double **a)
8371 static void diagnl(nrowa, diag, a)
8378 for (i = 1; i <=
nrowa; i++)
8380 for (j = i; j <=
nrowa; j++)
8395 static void error(
const char string[],
int *inform)
8397 static void error(
string, inform)
8398 const char string[];
8403 fprintf(stderr,
"%s\n",
string);
8415 estlam(
int nparam,
int neqn,
int *ifail,
double bigbnd,
double **hess,
8416 double *cvec,
double *a,
double *b,
struct _constraint *cs,
8417 double *psb,
double *bl,
double *bu,
double *x)
8420 estlam(nparam, neqn, ifail, bigbnd, hess, cvec, a, b, cs, psb, bl, bu, x)
8429 for (i = 1; i <=
neqn; i++)
8435 for (j = i; j <=
neqn; j++)
8439 hess[
j][
i] = hess[
i][
j];
8446 ctemp = convert(hess, neqn, neqn);
8450 ql0001_(&zero, &zero, &one, &neqn, &neqn, &mnn, (ctemp + 1), (cvec + 1),
8451 (a + 1), (b + 1), (bl + 1), (bu + 1), (x + 1), (w + 1), &iout, ifail,
8463 static double *colvec(
double **a,
int col,
int nrows)
8465 static double *colvec(a, col, nrows)
8473 temp = make_dv(nrows);
8474 for (i = 1;i <=
nrows;i++)
8475 temp[i] = a[i][col];
8484 static double scaprd(
int n,
double *x,
double *
y)
8486 static double scaprd(n, x, y)
8495 for (i = 1;i <=
n;i++)
8496 z = z + x[i] * y[i];
8505 static void fool(
double x,
double y,
double *
z)
8507 static void fool(x, y, z)
8519 static double smallNumber()
8521 double one, two,
z, tsmall;
8528 tsmall = tsmall / two;
8529 fool(tsmall, one, &z);
8532 return tsmall*two*two;
8540 static int fuscmp(
double val,
double thrshd)
8542 static int fuscmp(val, thrshd)
8548 if (fabs(val) <= thrshd)
8560 static int indexs(
int i,
int nfs)
8562 static int indexs(i, nfs)
8580 static void matrcp(ndima, a, ndimb, b)
8587 for (i = 1; i <=
ndima; i++)
8588 for (j = 1; j <=
ndima; j++)
8592 for (i = 1; i <=
ndimb; i++)
8605 static void matrvc(
int la,
int na,
double **a,
double *x,
double *y)
8607 static void matrvc(la, na, a, x, y)
8615 for (i = 1; i <= la; i++)
8618 for (j = 1; j <= na; j++)
8619 yi = yi + a[i][j] * x[j];
8630 static void nullvc(
int nparam,
double *x)
8632 static void nullvc(nparam, x)
8639 for (i = 1; i <=
nparam; i++)
8655 resign(
int n,
int neqn,
double *psf,
double *grdpsf,
double *penp,
8656 struct _constraint *cs,
double *signeq,
int job1,
int job2)
8659 resign(n, neqn, psf, grdpsf, penp, cs, signeq, job1, job2)
8660 int job1, job2,
n,
neqn;
8668 if (job2 == 10 || job2 == 12)
8670 for (i = 1; i <=
neqn; i++)
8672 if (job1 == 10 || job1 == 12)
8675 if (job2 == 10 || job2 == 12)
8676 *psf = *psf + cs[i+
nineq].
val * penp[
i];
8677 if (job1 == 10 || job1 == 20)
8679 for (j = 1; j <=
n; j++)
8680 cs[i+nineq].grad[j] = cs[i+nineq].grad[j] * signeq[i];
8682 if (job2 == 10 || job2 == 20)
8685 for (i = 1; i <=
n; i++)
8686 for (j = 1; j <=
neqn; j++)
8687 grdpsf[i] = grdpsf[i] + cs[j+nineq].grad[i] * penp[j];
8697 sbout1(FILE *
io,
int n,
const char *s1,
double z,
double *z1,
int job,
int level)
8699 static void sbout1(io, n, s1, z, z1, job, level)
8711 fprintf(io,
" %s\t %22.14e\n", s1, z);
8713 fprintf(io,
"\t\t\t %s\t %22.14e\n", s1, z);
8719 fprintf(io,
" %s\t %22.14e\n", s1, z1[1]);
8721 fprintf(io,
"\t\t\t %s\t %22.14e\n", s1, z1[1]);
8722 for (j = 2; j <=
n; j++)
8725 fprintf(io,
" \t\t\t %22.14e\n", z1[j]);
8727 fprintf(io,
" \t\t\t\t\t\t %22.14e\n", z1[j]);
8738 sbout2(FILE *io,
int n,
int i,
const char *s1,
const char *s2,
double *z)
8740 static void sbout2(io, n, i, s1, s2, z)
8744 const char *s1, *s2;
8749 fprintf(io,
"\t\t\t %8s %5d %1s\t %22.14e\n", s1, i, s2, z[1]);
8750 for (j = 2; j <=
n; j++)
8751 fprintf(io,
"\t\t\t\t\t\t %22.14e\n", z[j]);
8760 static void shift(
int n,
int ii,
int *iact)
8762 static void shift(n, ii, iact)
8770 for (j = 1; j <=
n; j++)
8774 for (k = j; k >= 2; k--)
8775 iact[k] = iact[k-1];
8790 slope(
int nob,
int nobL,
int neqn,
int nparam,
int feasb,
8791 struct _objective *ob,
double *grdpsf,
double *x,
double *y,
8792 double fmax,
double theta,
int job,
double *prev,
int old)
8795 slope(nob, nobL, neqn, nparam, feasb, ob, grdpsf, x, y, fmax, theta, job,
8804 double slope1, rhs,
rhog, grdftx, grdfty, diff, grpstx, grpsty;
8808 if (feasb && nob == 0)
8810 if (neqn == 0 || !feasb)
8817 grpstx = scaprd(nparam, grdpsf, x);
8818 grpsty = scaprd(nparam, grdpsf, y);
8820 for (i = 1; i <=
nob; i++)
8823 slope1 = prev[
i] + scaprd(nparam, ob[i].grad, x);
8825 slope1 = ob[
i].
val + scaprd(nparam, ob[i].grad, x);
8826 tslope =
DMAX1(tslope, slope1);
8828 tslope =
DMAX1(tslope, -slope1);
8830 tslope = tslope - fmax - grpstx;
8833 rhs = theta * tslope +
fmax;
8835 for (i = 1; i <=
nob; i++)
8837 grdftx = scaprd(nparam, ob[i].grad, x) - grpstx;
8838 grdfty = scaprd(nparam, ob[i].grad, y) - grpsty;
8839 diff = grdfty - grdftx;
8842 rhog =
DMIN1(rhog, (rhs - ob[i].val - grdftx) / diff);
8844 rhog =
DMIN1(rhog, -(rhs + ob[i].val + grdftx) / diff);
8855 static int element(
int *
set,
int length,
int index)
8857 static int element(
set, length, index)
8865 for (i = 1; i <=
length; i++)
8869 if (
set[i] == index)
8903 v = (
double *)calloc(len,
sizeof(
double));
8906 fprintf(stderr,
"Run-time error in make_dv");
8929 v = (
int *)calloc(len,
sizeof(
int));
8932 fprintf(stderr,
"Run-time error in make_iv");
8944 make_dm(
int rows,
int cols)
8958 temp = (
double **)calloc(rows,
sizeof(
double *));
8961 fprintf(stderr,
"Run-time error in make_dm");
8965 for (i = 1; i <= rows; i++)
8967 temp[
i] = (
double *)calloc(cols,
sizeof(
double));
8970 fprintf(stderr,
"Run-time error in make_dm");
8991 free((
char *)(v + 1));
9007 free((
char *)(v + 1));
9028 for (i = 1; i <= rows; i++)
9029 free((
char *)(m[i] + 1));
9030 free((
char *)(m + 1));
9040 convert(
double **a,
int m,
int n)
9051 temp = make_dv(m * n);
9053 for (i = 1; i <=
n; i++)
9054 for (j = 1; j <=
m; j++)
9055 temp[(m*(i-1)+j)] = a[
j][
i];
9063 void wtn(
double a[],
unsigned long nn[],
int ndim,
int isign,
9064 void(*wtstep)(
double [],
unsigned long,
int))
9066 unsigned long i1, i2, i3,
k,
n, nnew, nprev = 1, nt,
ntot = 1;
9069 for (idim = 1;idim <= ndim;idim++)
9071 std::vector<double>
buffer(ntot);
9072 auto *wksp= buffer.data()-1;
9073 for (idim = 1;idim <= ndim;idim++)
9079 for (i2 = 0;i2 <
ntot;i2 += nnew)
9081 for (i1 = 1;i1 <= nprev;i1++)
9083 for (i3 = i1 + i2, k = 1;k <=
n;++
k, i3 += nprev)
9087 for (nt = n;nt >= 4;nt >>= 1)
9088 (*wtstep)(wksp, nt, isign);
9092 for (nt = 4;nt <=
n;nt <<= 1)
9093 (*wtstep)(wksp, nt, isign);
9096 for (i3 = i1 + i2, k = 1;k <=
n;++
k, i3 += nprev)
9107 unsigned int ncof, ioff, joff;
9118 static double c2[3] =
9120 0.0, 0.707106781186547, 0.707106781186547
9122 static double c4[5] =
9124 0.0, 0.4829629131445341, 0.8365163037378079,
9125 0.2241438680420134, -0.1294095225512604
9127 static double c12[13] =
9129 0.0, 0.111540743350, 0.494623890398, 0.751133908021,
9130 0.315250351709, -0.226264693965, -0.129766867567,
9131 0.097501605587, 0.027522865530, -0.031582039318,
9132 0.000553842201, 0.004777257511, -0.001077301085
9134 static double c20[21] =
9136 0.0, 0.026670057901, 0.188176800078, 0.527201188932,
9137 0.688459039454, 0.281172343661, -0.249846424327,
9138 -0.195946274377, 0.127369340336, 0.093057364604,
9139 -0.071394147166, -0.029457536822, 0.033212674059,
9140 0.003606553567, -0.010733175483, 0.001395351747,
9141 0.001992405295, -0.000685856695, -0.000116466855,
9142 0.000093588670, -0.000013264203
9144 static double c2r[2], c4r[5], c12r[13], c20r[21];
9168 nrerror(
"unimplemented value n in pwtset");
9169 for (k = 1;k <=
n;k++)
9171 wfilt.cr[wfilt.ncof+1-
k] = sig * wfilt.cc[
k];
9174 wfilt.ioff = wfilt.joff = -(n >> 1);
9177 void pwt(
double a[],
unsigned long n,
int isign)
9180 unsigned long i, ii, jf, jr,
k, n1,
ni, nj, nh, nmod;
9184 std::vector<double>
buffer(n, 0.0);
9185 auto *wksp= buffer.data()-1;
9186 nmod = wfilt.ncof *
n;
9191 for (ii = 1, i = 1;i <=
n;i += 2, ii++)
9193 ni = i + nmod + wfilt.ioff;
9194 nj = i + nmod + wfilt.joff;
9195 double &aux1=wksp[ii];
9196 double &aux2=wksp[ii+nh];
9197 unsigned long kmax=4*(wfilt.ncof/4);
9199 for (k = 1;k <= kmax;k+=4)
9201 unsigned long k_1=k+1;
9202 unsigned long k_2=k+2;
9203 unsigned long k_3=k+3;
9206 aux1 += wfilt.cc[
k] * a[jf+1];
9207 aux2 += wfilt.cr[
k] * a[jr+1];
9208 unsigned long jf_1 = n1 & (ni + k_1);
9209 unsigned long jr_1 = n1 & (nj + k_1);
9210 aux1 += wfilt.cc[k_1] * a[jf_1+1];
9211 aux2 += wfilt.cr[k_1] * a[jr_1+1];
9212 unsigned long jf_2 = n1 & (ni + k_2);
9213 unsigned long jr_2 = n1 & (nj + k_2);
9214 aux1 += wfilt.cc[k_2] * a[jf_2+1];
9215 aux2 += wfilt.cr[k_2] * a[jr_2+1];
9216 unsigned long jf_3 = n1 & (ni + k_3);
9217 unsigned long jr_3 = n1 & (nj + k_3);
9218 aux1 += wfilt.cc[k_3] * a[jf_3+1];
9219 aux2 += wfilt.cr[k_3] * a[jr_3+1];
9222 for (k = kmax+1;k <= wfilt.ncof;++
k)
9226 aux1 += wfilt.cc[
k] * a[jf+1];
9227 aux2 += wfilt.cr[
k] * a[jr+1];
9233 for (ii = 1, i = 1;i <=
n;i += 2, ii++)
9237 ni = i + nmod + wfilt.ioff;
9238 nj = i + nmod + wfilt.joff;
9239 for (k = 1;k <= wfilt.ncof;k++)
9241 jf = (n1 & (ni +
k)) + 1;
9242 jr = (n1 & (nj +
k)) + 1;
9243 wksp[jf] += wfilt.cc[
k] * ai;
9244 wksp[jr] += wfilt.cr[
k] * ai1;
9248 memcpy(&a[1],&wksp[1],n*
sizeof(
double));
9255 void gser(
double *gamser,
double a,
double x,
double *gln)
9258 double sum,
del, ap;
9264 nrerror(
"x less than 0 in routine gser");
9271 del = sum = 1.0 /
a;
9272 for (n = 1;n <=
ITMAX;n++)
9277 if (fabs(del) < fabs(sum)*
EPS)
9279 *gamser = sum * exp(-x + a *
log(x) - (*gln));
9283 nrerror(
"a too large, ITMAX too small in routine gser");
9292 #define FPMIN 1.0e-30 9294 void gcf(
double *gammcf,
double a,
double x,
double *gln)
9304 for (i = 1;i <=
ITMAX;i++)
9309 if (fabs(d) <
FPMIN)
9312 if (fabs(c) <
FPMIN)
9317 if (fabs(del - 1.0) <
EPS)
9321 nrerror(
"a too large, ITMAX too small in gcf");
9322 *gammcf = exp(-x + a *
log(x) - (*gln)) *
h;
9330 double gamser, gammcf, gln;
9332 if (x < 0.0 || a <= 0.0)
9333 nrerror(
"Invalid arguments in routine gammp");
9336 gser(&gamser, a, x, &gln);
9341 gcf(&gammcf, a, x, &gln);
9352 for (i = 1;i <=
n;i++)
9354 for (j = i;j <=
n;j++)
9356 for (sum = a[i*n+j], k = i - 1;k >= 1;k--)
9357 sum -= a[i*n+k] * a[j*n+k];
9365 a[j*n+
i] = sum / p[
i];
9370 void cholsl(
double *a,
int n,
double *p,
double *b,
double *x)
9375 for (i = 1;i <=
n;i++)
9377 for (sum = b[i], k = i - 1;k >= 1;k--)
9378 sum -= a[i*n+k] * x[k];
9381 for (i = n;i >= 1;i--)
9383 for (sum = x[i], k = i + 1;k <=
n;k++)
9384 sum -= a[k*n+i] * x[k];
9390 void polint(
double *xa,
double *ya,
int n,
double x,
double &y,
double &dy)
9393 double den, dif, dift, ho, hp,
w;
9394 dif = fabs(x - xa[1]);
9395 std::vector<double>
buffer(2*n);
9396 auto *c= buffer.data()-1;
9398 for (i = 1;i <=
n;i++)
9400 if ((dift = fabs(x - xa[i])) < dif)
9409 for (m = 1;m <
n;m++)
9411 for (i = 1;i <= n -
m;i++)
9416 if ((den = ho - hp) == 0.0)
9418 nrerror(
"error in routine polint\n");
9424 y += (dy = (2 * ns < (n -
m) ? c[ns+1] : d[ns--]));
void cholsl(double *a, int n, double *p, double *b, double *x)
int ql0001_(m, me, mmax, n, nmax, mnn, c, d, a, b, xl, xu, x, u, iout, ifail, iprint, war, lwar, iwar, liwar, eps1) integer *m
void svbksb(double *u, double *w, double *v, int m, int n, double *b, double *x)
struct Tglob_info glob_info
void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc, double(*func)(double *, void *), void *prm, int ncom, double *pcom, double *xicom)
double betai(double a, double b, double x)
void wtn(double a[], unsigned long nn[], int ndim, int isign, void(*wtstep)(double [], unsigned long, int))
struct _violation * sip_viol
void ksone(double data[], int n, double(*func)(double), double *d, double *prob)
double probks(double alam)
void svdcmp(double *U, int Lines, int Columns, double *W, double *V)
double brent(double ax, double bx, double cx, double(*func)(double *, void *), void *prm, double tol, double *xmin, int ncom, double *pcom, double *xicom)
void sqrt(Image< double > &op)
sbout1(glob_prnt.io, nparam, "multipliers for x ", dummy, param->mult, 2, 2)
diagnl(nqpram, eta, hess1)
void choldc(double *a, int n, double *p)
int nnlsWght(int N, int M, double *A, double *b, double *weight)
double betacf(double a, double b, double x)
void nrerror(const char error_text[])
void(*)(*)(*)(*) gradcn()
void gcf(double *gammcf, double a, double x, double *gln)
dealloc(nineq, neq, signeq, indxcn, indxob, cs, param)
void _nnls_g1(double a, double b, double *cterm, double *sterm, double *sig)
i<=ncnstr;i++) cs[i].mult=0.e0;if(nfsip) for(i=1;i<=nob;i++) { ob[i].mult=0.e0;ob[i].mult_L=0.e0;} for(i=1;i<=nqpram;i++) { ii=k+i;if(clamda[ii]==0.e0 &&clamda[ii+nqpram]==0.e0) continue;else if(clamda[ii] !=0.e0) clamda[ii]=-clamda[ii];else clamda[ii]=clamda[ii+nqpram];} nqnp=nqpram+ncnstr;for(i=1;i<=nqpram;i++) param-> mult[i]
matrvc(nparam, nparam, hess, di, w)
double chebev(double a, double b, double c[], int m, double x)
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
matrcp(nparam, hess, nparam+1, hess1)
cfsqp1(miter, nparam, nob, nobL, nfsip1, nineqn, neq, neqn, ncsipl1, ncsipn1, mesh_pts1, ncnstr, nctotl, nrowa, feasb, epskt, epseqn, indxob, indxcn, param, cs, ob, signeq, obj, constr, gradob, gradcn)
void linmin(double *p, double *xi, int n, double &fret, double(*func)(double *, void *), void *prm)
void log(Image< double > &op)
void cfsqp(nparam, nf, nfsr, nineqn, nineq, neqn, neq, ncsrl, ncsrn, mesh_pts, mode, iprint, miter, inform, bigbnd, eps, epseqn, udelta, bl, bu, x, f, g, lambda, obj, constr, gradob, gradcn, cd) int nparam
struct Tglob_prnt glob_prnt
if(fabs(c[*nmax+ *nmax *c_dim1])==0.e0)
double bessi1_5(double x)
void indexx(int n, double arrin[], int indx[])
void beschb(double x, double *gam1, double *gam2, double *gampl, double *gammi)
void powell(double *p, double *xi, int n, double ftol, int &iter, double &fret, double(*func)(double *, void *), void *prm, bool show)
void polint(double *xa, double *ya, int n, double x, double &y, double &dy)
__host__ __device__ float length(float2 v)
double gammp(double a, double x)
void sort(struct DCEL_T *dcel)
void bessjy(double x, double xnu, double *rj, double *ry, double *rjp, double *ryp)
double bessj1_5(double x)
double bessi2_5(double x)
struct _parameter * param
double tdev(double nu, int *idum)
double bessj3_5(double x)
struct Tglob_log glob_log
int nnls(double *a, int m, int n, double *b, double *x, double *rnorm, double *wp, double *zzp, int *indexp)
dqp(nparam, nqprm0, nob, nobL, nfsip, nineqn, neq, neqn, nn, ncsipl, ncsipn, ncnstr, nctot0, nrowa0, nineqn, &infoqp, param, di, feasb, ob, *fmax, grdpsf, cs, a, cvec, bl, bu, clamda, hess, hess1, di, vv, 0)
struct Tglob_grd glob_grd
void gser(double *gamser, double a, double x, double *gln)
fprintf(glob_prnt.io, "\)
double bessi3_5(double x)
double bessi0_5(double x)
int _nnls_h12(int mode, int lpivot, int l1, int m, double *u, int u_dim1, double *up, double *cm, int ice, int icv, int ncv)
void pwt(double a[], unsigned long n, int isign)
double Pythag(double a, double b)
ql0001_ & zero(ctemp+1),(cvec+1),(a+1),(b+1),(bl+1),(bu+1),(x+1),(w+1), &iout, ifail, &zero,(w+3), &lwar2,(iw+1), &leniw, &glob_grd.epsmac
check(nparam, nf, nfsr, &Linfty, nineq, nineqn, neq, neqn, ncsrl, ncsrn, mode, &modem, eps, bgbnd, param)