| |
| 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_ */ |
| |
