From a3034ee0790fb65f69f024335cf35acbf6de5a44 Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Tue, 26 Jan 2010 06:05:01 +0100 Subject: cleaning --- .../Eigen/src/NonLinearOptimization/r1updt.h | 179 +++++++++------------ 1 file changed, 75 insertions(+), 104 deletions(-) (limited to 'unsupported/Eigen/src/NonLinearOptimization/r1updt.h') diff --git a/unsupported/Eigen/src/NonLinearOptimization/r1updt.h b/unsupported/Eigen/src/NonLinearOptimization/r1updt.h index e5e907582..71a704d24 100644 --- a/unsupported/Eigen/src/NonLinearOptimization/r1updt.h +++ b/unsupported/Eigen/src/NonLinearOptimization/r1updt.h @@ -25,119 +25,90 @@ void ei_r1updt(int m, int n, Matrix< Scalar, Dynamic, Dynamic > &s, const Scalar /* rotate the vector v into a multiple of the n-th unit vector */ /* in such a way that a spike is introduced into w. */ - nm1 = n - 1; - if (nm1 < 1) { - goto L70; - } - for (nmj = 1; nmj <= nm1; ++nmj) { - j = n - nmj; - w[j] = 0.; - if (v[j] == 0.) { - goto L50; - } - - /* determine a givens rotation which eliminates the */ - /* j-th element of v. */ - - if (ei_abs(v[n]) >= ei_abs(v[j])) - goto L20; - cotan = v[n] / v[j]; - /* Computing 2nd power */ - sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); - cos__ = sin__ * cotan; - tau = 1.; - if (ei_abs(cos__) * giant > 1.) { - tau = 1. / cos__; - } - goto L30; -L20: - tan__ = v[j] / v[n]; - /* Computing 2nd power */ - cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); - sin__ = cos__ * tan__; - tau = sin__; -L30: - - /* apply the transformation to v and store the information */ - /* necessary to recover the givens rotation. */ - - v[n] = sin__ * v[j] + cos__ * v[n]; - v[j] = tau; - - /* apply the transformation to s and extend the spike in w. */ - - for (i = j; i <= m; ++i) { - temp = cos__ * s(j-1,i-1) - sin__ * w[i]; - w[i] = sin__ * s(j-1,i-1) + cos__ * w[i]; - s(j-1,i-1) = temp; + if (nm1 >= 1) + for (nmj = 1; nmj <= nm1; ++nmj) { + j = n - nmj; + w[j] = 0.; + if (v[j] != 0.) { + /* determine a givens rotation which eliminates the */ + /* j-th element of v. */ + if (ei_abs(v[n]) < ei_abs(v[j])) { + cotan = v[n] / v[j]; + /* Computing 2nd power */ + sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); + cos__ = sin__ * cotan; + tau = 1.; + if (ei_abs(cos__) * giant > 1.) { + tau = 1. / cos__; + } + } else { + tan__ = v[j] / v[n]; + /* Computing 2nd power */ + cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); + sin__ = cos__ * tan__; + tau = sin__; + } + + /* apply the transformation to v and store the information */ + /* necessary to recover the givens rotation. */ + v[n] = sin__ * v[j] + cos__ * v[n]; + v[j] = tau; + + /* apply the transformation to s and extend the spike in w. */ + for (i = j; i <= m; ++i) { + temp = cos__ * s(j-1,i-1) - sin__ * w[i]; + w[i] = sin__ * s(j-1,i-1) + cos__ * w[i]; + s(j-1,i-1) = temp; + } + } } -L50: - ; - } -L70: /* add the spike from the rank 1 update to w. */ - - for (i = 1; i <= m; ++i) { + for (i = 1; i <= m; ++i) w[i] += v[n] * u[i]; - /* L80: */ - } /* eliminate the spike. */ - *sing = false; - if (nm1 < 1) { - goto L140; - } - for (j = 1; j <= nm1; ++j) { - if (w[j] == 0.) { - goto L120; - } - - /* determine a givens rotation which eliminates the */ - /* j-th element of the spike. */ - - if (ei_abs(s(j-1,j-1)) >= ei_abs(w[j])) - goto L90; - cotan = s(j-1,j-1) / w[j]; - /* Computing 2nd power */ - sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); - cos__ = sin__ * cotan; - tau = 1.; - if (ei_abs(cos__) * giant > 1.) { - tau = 1. / cos__; - } - goto L100; -L90: - tan__ = w[j] / s(j-1,j-1); - /* Computing 2nd power */ - cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); - sin__ = cos__ * tan__; - tau = sin__; -L100: - - /* apply the transformation to s and reduce the spike in w. */ - - for (i = j; i <= m; ++i) { - temp = cos__ * s(j-1,i-1) + sin__ * w[i]; - w[i] = -sin__ * s(j-1,i-1) + cos__ * w[i]; - s(j-1,i-1) = temp; + if (nm1 >= 1) + for (j = 1; j <= nm1; ++j) { + if (w[j] != 0.) { + /* determine a givens rotation which eliminates the */ + /* j-th element of the spike. */ + if (ei_abs(s(j-1,j-1)) < ei_abs(w[j])) { + cotan = s(j-1,j-1) / w[j]; + /* Computing 2nd power */ + sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); + cos__ = sin__ * cotan; + tau = 1.; + if (ei_abs(cos__) * giant > 1.) { + tau = 1. / cos__; + } + } else { + tan__ = w[j] / s(j-1,j-1); + /* Computing 2nd power */ + cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); + sin__ = cos__ * tan__; + tau = sin__; + } + + /* apply the transformation to s and reduce the spike in w. */ + for (i = j; i <= m; ++i) { + temp = cos__ * s(j-1,i-1) + sin__ * w[i]; + w[i] = -sin__ * s(j-1,i-1) + cos__ * w[i]; + s(j-1,i-1) = temp; + } + + /* store the information necessary to recover the */ + /* givens rotation. */ + w[j] = tau; + } + + /* test for zero diagonal elements in the output s. */ + if (s(j-1,j-1) == 0.) { + *sing = true; + } } - - /* store the information necessary to recover the */ - /* givens rotation. */ - - w[j] = tau; -L120: - - /* test for zero diagonal elements in the output s. */ - - if (s(j-1,j-1) == 0.) { - *sing = true; - } - } -L140: /* move w back into the last column of the output s. */ s(n-1,n-1) = w[n]; -- cgit v1.2.3