misc cleaning

This commit is contained in:
Thomas Capricelli 2010-01-26 13:20:24 +01:00
parent 71f513e250
commit 55f328b264
6 changed files with 48 additions and 80 deletions

View File

@ -17,14 +17,13 @@ void ei_chkder(
const Scalar epsf = chkder_factor * epsilon<Scalar>(); const Scalar epsf = chkder_factor * epsilon<Scalar>();
const Scalar epslog = chkder_log10e * ei_log(eps); const Scalar epslog = chkder_log10e * ei_log(eps);
Scalar temp; Scalar temp;
int i,j;
const int m = fvec.size(), n = x.size(); const int m = fvec.size(), n = x.size();
if (mode != 2) { if (mode != 2) {
/* mode = 1. */
xp.resize(n); xp.resize(n);
/* mode = 1. */ for (int j = 0; j < n; ++j) {
for (j = 0; j < n; ++j) {
temp = eps * ei_abs(x[j]); temp = eps * ei_abs(x[j]);
if (temp == 0.) if (temp == 0.)
temp = eps; temp = eps;
@ -32,15 +31,15 @@ void ei_chkder(
} }
} }
else { else {
/* mode = 2. */ /* mode = 2. */
err.setZero(m); err.setZero(m);
for (j = 0; j < n; ++j) { for (int j = 0; j < n; ++j) {
temp = ei_abs(x[j]); temp = ei_abs(x[j]);
if (temp == 0.) if (temp == 0.)
temp = 1.; temp = 1.;
err += temp * fjac.col(j); err += temp * fjac.col(j);
} }
for (i = 0; i < m; ++i) { for (int i = 0; i < m; ++i) {
temp = 1.; temp = 1.;
if (fvec[i] != 0. && fvecp[i] != 0. && ei_abs(fvecp[i] - fvec[i]) >= epsf * ei_abs(fvec[i])) if (fvec[i] != 0. && fvecp[i] != 0. && ei_abs(fvecp[i] - fvec[i]) >= epsf * ei_abs(fvec[i]))
temp = eps * ei_abs((fvecp[i] - fvec[i]) / eps - err[i]) / (ei_abs(fvec[i]) + ei_abs(fvecp[i])); temp = eps * ei_abs((fvecp[i] - fvec[i]) / eps - err[i]) / (ei_abs(fvec[i]) + ei_abs(fvecp[i]));
@ -51,5 +50,5 @@ void ei_chkder(
err[i] = 0.; err[i] = 0.;
} }
} }
} /* chkder_ */ }

View File

@ -1,5 +1,5 @@
template <typename Scalar> template <typename Scalar>
void ei_covar( void ei_covar(
Matrix< Scalar, Dynamic, Dynamic > &r, Matrix< Scalar, Dynamic, Dynamic > &r,
const VectorXi &ipvt, const VectorXi &ipvt,
@ -17,7 +17,6 @@ void ei_covar(
assert(ipvt.size()==n); assert(ipvt.size()==n);
/* form the inverse of r in the full upper triangle of r. */ /* form the inverse of r in the full upper triangle of r. */
l = -1; l = -1;
for (k = 0; k < n; ++k) for (k = 0; k < n; ++k)
if (ei_abs(r(k,k)) > tolr) { if (ei_abs(r(k,k)) > tolr) {
@ -33,7 +32,6 @@ void ei_covar(
/* form the full upper triangle of the inverse of (r transpose)*r */ /* form the full upper triangle of the inverse of (r transpose)*r */
/* in the full upper triangle of r. */ /* in the full upper triangle of r. */
for (k = 0; k <= l; ++k) { for (k = 0; k <= l; ++k) {
for (j = 0; j <= k-1; ++j) { for (j = 0; j <= k-1; ++j) {
temp = r(j,k); temp = r(j,k);
@ -47,7 +45,6 @@ void ei_covar(
/* form the full lower triangle of the covariance matrix */ /* form the full lower triangle of the covariance matrix */
/* in the strict lower triangle of r and in wa. */ /* in the strict lower triangle of r and in wa. */
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
jj = ipvt[j]; jj = ipvt[j];
sing = j > l; sing = j > l;
@ -64,7 +61,6 @@ void ei_covar(
} }
/* symmetrize the covariance matrix in r. */ /* symmetrize the covariance matrix in r. */
for (j = 0; j < n; ++j) { for (j = 0; j < n; ++j) {
for (i = 0; i <= j; ++i) for (i = 0; i <= j; ++i)
r(i,j) = r(j,i); r(i,j) = r(j,i);

View File

@ -13,7 +13,7 @@ int ei_fdjac1(
int i, j, k; int i, j, k;
Scalar eps, temp; Scalar eps, temp;
int msum; int msum;
int iflag = 0; int iflag;
/* Function Body */ /* Function Body */
const Scalar epsmch = epsilon<Scalar>(); const Scalar epsmch = epsilon<Scalar>();
@ -42,7 +42,7 @@ int ei_fdjac1(
}else { }else {
/* computation of banded approximate jacobian. */ /* computation of banded approximate jacobian. */
for (k = 0; k < msum; ++k) { for (k = 0; k < msum; ++k) {
for (j = k; msum< 0 ? j > n: j < n; j += msum) { for (j = k; (msum<0) ? (j>n): (j<n); j += msum) {
wa2[j] = x[j]; wa2[j] = x[j];
h = eps * ei_abs(wa2[j]); h = eps * ei_abs(wa2[j]);
if (h == 0.) h = eps; if (h == 0.) h = eps;
@ -51,7 +51,7 @@ int ei_fdjac1(
iflag = Functor(x, wa1); iflag = Functor(x, wa1);
if (iflag < 0) if (iflag < 0)
return iflag; return iflag;
for (j = k; msum< 0 ? j > n: j < n; j += msum) { for (j = k; (msum<0) ? (j>n): (j<n); j += msum) {
x[j] = wa2[j]; x[j] = wa2[j];
h = eps * ei_abs(wa2[j]); h = eps * ei_abs(wa2[j]);
if (h == 0.) h = eps; if (h == 0.) h = eps;
@ -63,6 +63,6 @@ int ei_fdjac1(
} }
} }
} }
return iflag; return 0;
} /* fdjac1_ */ }

View File

@ -22,7 +22,6 @@ void ei_qrsolv(
/* copy r and (q transpose)*b to preserve input and initialize s. */ /* copy r and (q transpose)*b to preserve input and initialize s. */
/* in particular, save the diagonal elements of r in x. */ /* in particular, save the diagonal elements of r in x. */
x = s.diagonal(); x = s.diagonal();
wa = qtb; wa = qtb;
@ -35,7 +34,6 @@ void ei_qrsolv(
/* prepare the row of d to be eliminated, locating the */ /* prepare the row of d to be eliminated, locating the */
/* diagonal element using p from the qr factorization. */ /* diagonal element using p from the qr factorization. */
l = ipvt[j]; l = ipvt[j];
if (diag[l] == 0.) if (diag[l] == 0.)
break; break;
@ -45,7 +43,6 @@ void ei_qrsolv(
/* the transformations to eliminate the row of d */ /* the transformations to eliminate the row of d */
/* modify only a single element of (q transpose)*b */ /* modify only a single element of (q transpose)*b */
/* beyond the first n, which is initially zero. */ /* beyond the first n, which is initially zero. */
Scalar qtbpj = 0.; Scalar qtbpj = 0.;
for (k = j; k < n; ++k) { for (k = j; k < n; ++k) {
/* determine a givens rotation which eliminates the */ /* determine a givens rotation which eliminates the */
@ -55,7 +52,6 @@ void ei_qrsolv(
/* compute the modified diagonal element of r and */ /* compute the modified diagonal element of r and */
/* the modified element of ((q transpose)*b,0). */ /* the modified element of ((q transpose)*b,0). */
s(k,k) = givens.c() * s(k,k) + givens.s() * sdiag[k]; s(k,k) = givens.c() * s(k,k) + givens.s() * sdiag[k];
temp = givens.c() * wa[k] + givens.s() * qtbpj; temp = givens.c() * wa[k] + givens.s() * qtbpj;
qtbpj = -givens.s() * wa[k] + givens.c() * qtbpj; qtbpj = -givens.s() * wa[k] + givens.c() * qtbpj;
@ -72,7 +68,6 @@ void ei_qrsolv(
/* solve the triangular system for z. if the system is */ /* solve the triangular system for z. if the system is */
/* singular, then obtain a least squares solution. */ /* singular, then obtain a least squares solution. */
int nsing; int nsing;
for (nsing=0; nsing<n && sdiag[nsing]!=0; nsing++); for (nsing=0; nsing<n && sdiag[nsing]!=0; nsing++);
wa.segment(nsing,n-nsing).setZero(); wa.segment(nsing,n-nsing).setZero();

View File

@ -7,7 +7,7 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
int a_dim1, a_offset; int a_dim1, a_offset;
/* Local variables */ /* Local variables */
int i, j, nm1, nmj; int i, j, nmj;
Scalar cos__=0., sin__=0., temp; Scalar cos__=0., sin__=0., temp;
/* Parameter adjustments */ /* Parameter adjustments */
@ -18,12 +18,11 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
a -= a_offset; a -= a_offset;
/* Function Body */ /* Function Body */
nm1 = n - 1; if (n<=1)
if (nm1 < 1)
return; return;
/* apply the first set of givens rotations to a. */ /* apply the first set of givens rotations to a. */
for (nmj = 1; nmj <= nm1; ++nmj) { for (nmj = 1; nmj <= n-1; ++nmj) {
j = n - nmj; j = n - nmj;
if (ei_abs(v[j]) > 1.) { if (ei_abs(v[j]) > 1.) {
cos__ = 1. / v[j]; cos__ = 1. / v[j];
@ -40,7 +39,7 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
} }
} }
/* apply the second set of givens rotations to a. */ /* apply the second set of givens rotations to a. */
for (j = 1; j <= nm1; ++j) { for (j = 1; j <= n-1; ++j) {
if (ei_abs(w[j]) > 1.) { if (ei_abs(w[j]) > 1.) {
cos__ = 1. / w[j]; cos__ = 1. / w[j];
sin__ = ei_sqrt(1. - ei_abs2(cos__)); sin__ = ei_sqrt(1. - ei_abs2(cos__));
@ -56,5 +55,5 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
} }
} }
return; return;
} /* r1mpyq_ */ }

View File

@ -8,7 +8,6 @@ void ei_rwupdt(int n, Scalar *r__, int ldr,
int r_dim1, r_offset; int r_dim1, r_offset;
/* Local variables */ /* Local variables */
int i, j, jm1;
Scalar tan__, temp, rowj, cotan; Scalar tan__, temp, rowj, cotan;
/* Parameter adjustments */ /* Parameter adjustments */
@ -21,60 +20,40 @@ void ei_rwupdt(int n, Scalar *r__, int ldr,
r__ -= r_offset; r__ -= r_offset;
/* Function Body */ /* Function Body */
for (int j = 1; j <= n; ++j) {
rowj = w[j];
for (j = 1; j <= n; ++j) { /* apply the previous transformations to */
rowj = w[j]; /* r(i,j), i=1,2,...,j-1, and to w(j). */
jm1 = j - 1; if (j-1>=1)
for (int i = 1; i <= j-1; ++i) {
temp = cos__[i] * r__[i + j * r_dim1] + sin__[i] * rowj;
rowj = -sin__[i] * r__[i + j * r_dim1] + cos__[i] * rowj;
r__[i + j * r_dim1] = temp;
}
/* apply the previous transformations to */ /* determine a givens rotation which eliminates w(j). */
/* r(i,j), i=1,2,...,j-1, and to w(j). */ cos__[j] = 1.;
sin__[j] = 0.;
if (rowj != 0.) {
if (ei_abs(r__[j + j * r_dim1]) < ei_abs(rowj)) {
cotan = r__[j + j * r_dim1] / rowj;
sin__[j] = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
cos__[j] = sin__[j] * cotan;
}
else {
tan__ = rowj / r__[j + j * r_dim1];
cos__[j] = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
sin__[j] = cos__[j] * tan__;
}
if (jm1 < 1) { /* apply the current transformation to r(j,j), b(j), and alpha. */
goto L20; r__[j + j * r_dim1] = cos__[j] * r__[j + j * r_dim1] + sin__[j] * rowj;
} temp = cos__[j] * b[j] + sin__[j] * *alpha;
for (i = 1; i <= jm1; ++i) { *alpha = -sin__[j] * b[j] + cos__[j] * *alpha;
temp = cos__[i] * r__[i + j * r_dim1] + sin__[i] * rowj; b[j] = temp;
rowj = -sin__[i] * r__[i + j * r_dim1] + cos__[i] * rowj; }
r__[i + j * r_dim1] = temp;
/* L10: */
}
L20:
/* determine a givens rotation which eliminates w(j). */
cos__[j] = 1.;
sin__[j] = 0.;
if (rowj == 0.) {
goto L50;
}
if (ei_abs(r__[j + j * r_dim1]) >= ei_abs(rowj))
goto L30;
cotan = r__[j + j * r_dim1] / rowj;
/* Computing 2nd power */
sin__[j] = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
cos__[j] = sin__[j] * cotan;
goto L40;
L30:
tan__ = rowj / r__[j + j * r_dim1];
/* Computing 2nd power */
cos__[j] = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
sin__[j] = cos__[j] * tan__;
L40:
/* apply the current transformation to r(j,j), b(j), and alpha. */
r__[j + j * r_dim1] = cos__[j] * r__[j + j * r_dim1] + sin__[j] *
rowj;
temp = cos__[j] * b[j] + sin__[j] * *alpha;
*alpha = -sin__[j] * b[j] + cos__[j] * *alpha;
b[j] = temp;
L50:
/* L60: */
;
} }
return; return;
}
/* last card of subroutine rwupdt. */
} /* rwupdt_ */