covar : cleaning, removing goto's

This commit is contained in:
Thomas Capricelli 2009-08-24 16:49:38 +02:00
parent 312ab1abb3
commit f69869c42a

View File

@ -1,7 +1,7 @@
template <typename Scalar> template <typename Scalar>
void ei_covar(int n, Scalar *r__, int ldr, void ei_covar(int n, Scalar *r__, int ldr,
const int *ipvt, Scalar tol, Scalar *wa) const int *ipvt, Scalar tol, Scalar *wa)
{ {
/* System generated locals */ /* System generated locals */
int r_dim1, r_offset; int r_dim1, r_offset;
@ -21,98 +21,66 @@ void ei_covar(int n, Scalar *r__, int ldr,
/* Function Body */ /* Function Body */
/* 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 = 0; l = 0;
for (k = 1; k <= n; ++k) { for (k = 1; k <= n; ++k) {
if (ei_abs(r__[k + k * r_dim1]) <= tolr) { if (ei_abs(r__[k + k * r_dim1]) > tolr) {
goto L50; r__[k + k * r_dim1] = 1. / r__[k + k * r_dim1];
} km1 = k - 1;
r__[k + k * r_dim1] = 1. / r__[k + k * r_dim1]; if (km1 >= 1)
km1 = k - 1; for (j = 1; j <= km1; ++j) {
if (km1 < 1) { temp = r__[k + k * r_dim1] * r__[j + k * r_dim1];
goto L30; r__[j + k * r_dim1] = 0.;
} for (i = 1; i <= j; ++i) {
for (j = 1; j <= km1; ++j) { r__[i + k * r_dim1] -= temp * r__[i + j * r_dim1];
temp = r__[k + k * r_dim1] * r__[j + k * r_dim1]; }
r__[j + k * r_dim1] = 0.; }
for (i = 1; i <= j; ++i) { l = k;
r__[i + k * r_dim1] -= temp * r__[i + j * r_dim1]; }
/* L10: */
}
/* L20: */
}
L30:
l = k;
/* L40: */
} }
L50:
/* 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. */
if (l < 1) { if (l >= 1)
goto L110; for (k = 1; k <= l; ++k) {
} km1 = k - 1;
for (k = 1; k <= l; ++k) { if (km1 >= 1)
km1 = k - 1; for (j = 1; j <= km1; ++j) {
if (km1 < 1) { temp = r__[j + k * r_dim1];
goto L80; for (i = 1; i <= j; ++i)
} r__[i + j * r_dim1] += temp * r__[i + k * r_dim1];
for (j = 1; j <= km1; ++j) { }
temp = r__[j + k * r_dim1]; temp = r__[k + k * r_dim1];
for (i = 1; i <= j; ++i) { for (i = 1; i <= k; ++i)
r__[i + j * r_dim1] += temp * r__[i + k * r_dim1]; r__[i + k * r_dim1] = temp * r__[i + k * r_dim1];
/* L60: */ }
}
/* L70: */
}
L80:
temp = r__[k + k * r_dim1];
for (i = 1; i <= k; ++i) {
r__[i + k * r_dim1] = temp * r__[i + k * r_dim1];
/* L90: */
}
/* L100: */
}
L110:
/* 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 = 1; j <= n; ++j) { for (j = 1; j <= n; ++j) {
jj = ipvt[j]; jj = ipvt[j];
sing = j > l; sing = j > l;
for (i = 1; i <= j; ++i) { for (i = 1; i <= j; ++i) {
if (sing) { if (sing)
r__[i + j * r_dim1] = 0.; r__[i + j * r_dim1] = 0.;
} ii = ipvt[i];
ii = ipvt[i]; if (ii > jj)
if (ii > jj) { r__[ii + jj * r_dim1] = r__[i + j * r_dim1];
r__[ii + jj * r_dim1] = r__[i + j * r_dim1]; if (ii < jj)
} r__[jj + ii * r_dim1] = r__[i + j * r_dim1];
if (ii < jj) { }
r__[jj + ii * r_dim1] = r__[i + j * r_dim1]; wa[jj] = r__[j + j * r_dim1];
}
/* L120: */
}
wa[jj] = r__[j + j * r_dim1];
/* L130: */
} }
/* symmetrize the covariance matrix in r. */ /* symmetrize the covariance matrix in r. */
for (j = 1; j <= n; ++j) { for (j = 1; j <= n; ++j) {
for (i = 1; i <= j; ++i) { for (i = 1; i <= j; ++i)
r__[i + j * r_dim1] = r__[j + i * r_dim1]; r__[i + j * r_dim1] = r__[j + i * r_dim1];
/* L140: */ r__[j + j * r_dim1] = wa[j];
}
r__[j + j * r_dim1] = wa[j];
/* L150: */
} }
/*return 0;*/ }
/* last card of subroutine covar. */
} /* covar_ */