aboutsummaryrefslogtreecommitdiffhomepage
path: root/unsupported/Eigen/src/NonLinearOptimization/r1updt.h
diff options
context:
space:
mode:
authorGravatar Thomas Capricelli <orzel@freehackers.org>2009-11-09 03:27:15 +0100
committerGravatar Thomas Capricelli <orzel@freehackers.org>2009-11-09 03:27:15 +0100
commit71a3e96b49f7f770df9e44909e14379a9a75f2a1 (patch)
tree57e9115dc622ef1ee6dfcc47866d62b8f55a8061 /unsupported/Eigen/src/NonLinearOptimization/r1updt.h
parent09cb27c587c6ee240c9a24a10c131d27fb0a3ed8 (diff)
rename NonLinear to NonLinearOptimization
Diffstat (limited to 'unsupported/Eigen/src/NonLinearOptimization/r1updt.h')
-rw-r--r--unsupported/Eigen/src/NonLinearOptimization/r1updt.h175
1 files changed, 175 insertions, 0 deletions
diff --git a/unsupported/Eigen/src/NonLinearOptimization/r1updt.h b/unsupported/Eigen/src/NonLinearOptimization/r1updt.h
new file mode 100644
index 000000000..f2ddb189b
--- /dev/null
+++ b/unsupported/Eigen/src/NonLinearOptimization/r1updt.h
@@ -0,0 +1,175 @@
+
+ template <typename Scalar>
+void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v, Scalar *w, bool *sing)
+{
+ /* Local variables */
+ int i, j, l, jj, nm1;
+ Scalar tan__;
+ int nmj;
+ Scalar cos__, sin__, tau, temp, cotan;
+
+ /* Parameter adjustments */
+ --w;
+ --u;
+ --v;
+ --s;
+
+ /* Function Body */
+ const Scalar giant = std::numeric_limits<Scalar>::max();
+
+ /* initialize the diagonal element pointer. */
+
+ jj = n * ((m << 1) - n + 1) / 2 - (m - n);
+
+ /* move the nontrivial part of the last column of s into w. */
+
+ l = jj;
+ for (i = n; i <= m; ++i) {
+ w[i] = s[l];
+ ++l;
+ /* L10: */
+ }
+
+ /* 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;
+ jj -= m - j + 1;
+ 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. */
+
+ l = jj;
+ for (i = j; i <= m; ++i) {
+ temp = cos__ * s[l] - sin__ * w[i];
+ w[i] = sin__ * s[l] + cos__ * w[i];
+ s[l] = temp;
+ ++l;
+ /* L40: */
+ }
+L50:
+ /* L60: */
+ ;
+ }
+L70:
+
+ /* add the spike from the rank 1 update to w. */
+
+ 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[jj]) >= ei_abs(w[j]))
+ goto L90;
+ cotan = s[jj] / 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[jj];
+ /* 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. */
+
+ l = jj;
+ for (i = j; i <= m; ++i) {
+ temp = cos__ * s[l] + sin__ * w[i];
+ w[i] = -sin__ * s[l] + cos__ * w[i];
+ s[l] = temp;
+ ++l;
+ /* L110: */
+ }
+
+ /* store the information necessary to recover the */
+ /* givens rotation. */
+
+ w[j] = tau;
+L120:
+
+ /* test for zero diagonal elements in the output s. */
+
+ if (s[jj] == 0.) {
+ *sing = true;
+ }
+ jj += m - j + 1;
+ /* L130: */
+ }
+L140:
+
+ /* move w back into the last column of the output s. */
+
+ l = jj;
+ for (i = n; i <= m; ++i) {
+ s[l] = w[i];
+ ++l;
+ /* L150: */
+ }
+ if (s[jj] == 0.) {
+ *sing = true;
+ }
+ return;
+
+ /* last card of subroutine r1updt. */
+
+} /* r1updt_ */
+