From 6d8baa757e7f02e5752538882bc3ce8f578089ce Mon Sep 17 00:00:00 2001 From: Thomas Capricelli Date: Mon, 14 Sep 2009 23:47:44 +0200 Subject: [PATCH] fix indentation (and only that) --- unsupported/Eigen/src/NonLinear/qform.h | 88 +++++----- unsupported/Eigen/src/NonLinear/qrfac.h | 172 ++++++++++---------- unsupported/Eigen/src/NonLinear/qrsolv.h | 196 +++++++++++------------ unsupported/Eigen/src/NonLinear/r1mpyq.h | 102 ++++++------ unsupported/Eigen/src/NonLinear/r1updt.h | 196 +++++++++++------------ 5 files changed, 377 insertions(+), 377 deletions(-) diff --git a/unsupported/Eigen/src/NonLinear/qform.h b/unsupported/Eigen/src/NonLinear/qform.h index 844c8c2f6..47bbd1fbe 100644 --- a/unsupported/Eigen/src/NonLinear/qform.h +++ b/unsupported/Eigen/src/NonLinear/qform.h @@ -1,7 +1,7 @@ template void ei_qform(int m, int n, Scalar *q, int - ldq, Scalar *wa) + ldq, Scalar *wa) { /* System generated locals */ int q_dim1, q_offset; @@ -19,71 +19,71 @@ void ei_qform(int m, int n, Scalar *q, int /* Function Body */ -/* zero out upper triangle of q in the first min(m,n) columns. */ + /* zero out upper triangle of q in the first min(m,n) columns. */ minmn = std::min(m,n); if (minmn < 2) { - goto L30; + goto L30; } for (j = 2; j <= minmn; ++j) { - jm1 = j - 1; - for (i = 1; i <= jm1; ++i) { - q[i + j * q_dim1] = 0.; -/* L10: */ - } -/* L20: */ + jm1 = j - 1; + for (i = 1; i <= jm1; ++i) { + q[i + j * q_dim1] = 0.; + /* L10: */ + } + /* L20: */ } L30: -/* initialize remaining columns to those of the identity matrix. */ + /* initialize remaining columns to those of the identity matrix. */ np1 = n + 1; if (m < np1) { - goto L60; + goto L60; } for (j = np1; j <= m; ++j) { - for (i = 1; i <= m; ++i) { - q[i + j * q_dim1] = 0.; -/* L40: */ - } - q[j + j * q_dim1] = 1.; -/* L50: */ + for (i = 1; i <= m; ++i) { + q[i + j * q_dim1] = 0.; + /* L40: */ + } + q[j + j * q_dim1] = 1.; + /* L50: */ } L60: -/* accumulate q from its factored form. */ + /* accumulate q from its factored form. */ for (l = 1; l <= minmn; ++l) { - k = minmn - l + 1; - for (i = k; i <= m; ++i) { - wa[i] = q[i + k * q_dim1]; - q[i + k * q_dim1] = 0.; -/* L70: */ - } - q[k + k * q_dim1] = 1.; - if (wa[k] == 0.) { - goto L110; - } - for (j = k; j <= m; ++j) { - sum = 0.; - for (i = k; i <= m; ++i) { - sum += q[i + j * q_dim1] * wa[i]; -/* L80: */ - } - temp = sum / wa[k]; - for (i = k; i <= m; ++i) { - q[i + j * q_dim1] -= temp * wa[i]; -/* L90: */ - } -/* L100: */ - } + k = minmn - l + 1; + for (i = k; i <= m; ++i) { + wa[i] = q[i + k * q_dim1]; + q[i + k * q_dim1] = 0.; + /* L70: */ + } + q[k + k * q_dim1] = 1.; + if (wa[k] == 0.) { + goto L110; + } + for (j = k; j <= m; ++j) { + sum = 0.; + for (i = k; i <= m; ++i) { + sum += q[i + j * q_dim1] * wa[i]; + /* L80: */ + } + temp = sum / wa[k]; + for (i = k; i <= m; ++i) { + q[i + j * q_dim1] -= temp * wa[i]; + /* L90: */ + } + /* L100: */ + } L110: -/* L120: */ - ; + /* L120: */ + ; } return; -/* last card of subroutine qform. */ + /* last card of subroutine qform. */ } /* qform_ */ diff --git a/unsupported/Eigen/src/NonLinear/qrfac.h b/unsupported/Eigen/src/NonLinear/qrfac.h index 94e36bb26..c14fd0714 100644 --- a/unsupported/Eigen/src/NonLinear/qrfac.h +++ b/unsupported/Eigen/src/NonLinear/qrfac.h @@ -1,8 +1,8 @@ template void ei_qrfac(int m, int n, Scalar *a, int - lda, int pivot, int *ipvt, int /* lipvt */, Scalar *rdiag, - Scalar *acnorm) + lda, int pivot, int *ipvt, int /* lipvt */, Scalar *rdiag, + Scalar *acnorm) { /* System generated locals */ int a_dim1, a_offset; @@ -28,109 +28,109 @@ void ei_qrfac(int m, int n, Scalar *a, int /* Function Body */ const Scalar epsmch = epsilon(); -/* compute the initial column norms and initialize several arrays. */ + /* compute the initial column norms and initialize several arrays. */ for (j = 1; j <= n; ++j) { - acnorm[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm(); - rdiag[j] = acnorm[j]; - wa[j] = rdiag[j]; - if (pivot) { - ipvt[j] = j; - } -/* L10: */ + acnorm[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm(); + rdiag[j] = acnorm[j]; + wa[j] = rdiag[j]; + if (pivot) { + ipvt[j] = j; + } + /* L10: */ } -/* reduce a to r with householder transformations. */ + /* reduce a to r with householder transformations. */ minmn = std::min(m,n); for (j = 1; j <= minmn; ++j) { - if (! (pivot)) { - goto L40; - } + if (! (pivot)) { + goto L40; + } -/* bring the column of largest norm into the pivot position. */ + /* bring the column of largest norm into the pivot position. */ - kmax = j; - for (k = j; k <= n; ++k) { - if (rdiag[k] > rdiag[kmax]) { - kmax = k; - } -/* L20: */ - } - if (kmax == j) { - goto L40; - } - for (i = 1; i <= m; ++i) { - temp = a[i + j * a_dim1]; - a[i + j * a_dim1] = a[i + kmax * a_dim1]; - a[i + kmax * a_dim1] = temp; -/* L30: */ - } - rdiag[kmax] = rdiag[j]; - wa[kmax] = wa[j]; - k = ipvt[j]; - ipvt[j] = ipvt[kmax]; - ipvt[kmax] = k; + kmax = j; + for (k = j; k <= n; ++k) { + if (rdiag[k] > rdiag[kmax]) { + kmax = k; + } + /* L20: */ + } + if (kmax == j) { + goto L40; + } + for (i = 1; i <= m; ++i) { + temp = a[i + j * a_dim1]; + a[i + j * a_dim1] = a[i + kmax * a_dim1]; + a[i + kmax * a_dim1] = temp; + /* L30: */ + } + rdiag[kmax] = rdiag[j]; + wa[kmax] = wa[j]; + k = ipvt[j]; + ipvt[j] = ipvt[kmax]; + ipvt[kmax] = k; L40: -/* compute the householder transformation to reduce the */ -/* j-th column of a to a multiple of the j-th unit vector. */ + /* compute the householder transformation to reduce the */ + /* j-th column of a to a multiple of the j-th unit vector. */ - ajnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j + j * a_dim1],m-j+1).blueNorm(); - if (ajnorm == 0.) { - goto L100; - } - if (a[j + j * a_dim1] < 0.) { - ajnorm = -ajnorm; - } - for (i = j; i <= m; ++i) { - a[i + j * a_dim1] /= ajnorm; -/* L50: */ - } - a[j + j * a_dim1] += 1.; + ajnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j + j * a_dim1],m-j+1).blueNorm(); + if (ajnorm == 0.) { + goto L100; + } + if (a[j + j * a_dim1] < 0.) { + ajnorm = -ajnorm; + } + for (i = j; i <= m; ++i) { + a[i + j * a_dim1] /= ajnorm; + /* L50: */ + } + a[j + j * a_dim1] += 1.; -/* apply the transformation to the remaining columns */ -/* and update the norms. */ + /* apply the transformation to the remaining columns */ + /* and update the norms. */ - jp1 = j + 1; - if (n < jp1) { - goto L100; - } - for (k = jp1; k <= n; ++k) { - sum = 0.; - for (i = j; i <= m; ++i) { - sum += a[i + j * a_dim1] * a[i + k * a_dim1]; -/* L60: */ - } - temp = sum / a[j + j * a_dim1]; - for (i = j; i <= m; ++i) { - a[i + k * a_dim1] -= temp * a[i + j * a_dim1]; -/* L70: */ - } - if (! (pivot) || rdiag[k] == 0.) { - goto L80; - } - temp = a[j + k * a_dim1] / rdiag[k]; -/* Computing MAX */ -/* Computing 2nd power */ - rdiag[k] *= ei_sqrt((std::max(Scalar(0.), Scalar(1.)-ei_abs2(temp)))); -/* Computing 2nd power */ - if (Scalar(.05) * ei_abs2(rdiag[k] / wa[k]) > epsmch) { - goto L80; - } - rdiag[k] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[jp1 + k * a_dim1],m-j).blueNorm(); - wa[k] = rdiag[k]; + jp1 = j + 1; + if (n < jp1) { + goto L100; + } + for (k = jp1; k <= n; ++k) { + sum = 0.; + for (i = j; i <= m; ++i) { + sum += a[i + j * a_dim1] * a[i + k * a_dim1]; + /* L60: */ + } + temp = sum / a[j + j * a_dim1]; + for (i = j; i <= m; ++i) { + a[i + k * a_dim1] -= temp * a[i + j * a_dim1]; + /* L70: */ + } + if (! (pivot) || rdiag[k] == 0.) { + goto L80; + } + temp = a[j + k * a_dim1] / rdiag[k]; + /* Computing MAX */ + /* Computing 2nd power */ + rdiag[k] *= ei_sqrt((std::max(Scalar(0.), Scalar(1.)-ei_abs2(temp)))); + /* Computing 2nd power */ + if (Scalar(.05) * ei_abs2(rdiag[k] / wa[k]) > epsmch) { + goto L80; + } + rdiag[k] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[jp1 + k * a_dim1],m-j).blueNorm(); + wa[k] = rdiag[k]; L80: -/* L90: */ - ; - } + /* L90: */ + ; + } L100: - rdiag[j] = -ajnorm; -/* L110: */ + rdiag[j] = -ajnorm; + /* L110: */ } return; -/* last card of subroutine qrfac. */ + /* last card of subroutine qrfac. */ } /* qrfac_ */ diff --git a/unsupported/Eigen/src/NonLinear/qrsolv.h b/unsupported/Eigen/src/NonLinear/qrsolv.h index 778c66a82..57884870c 100644 --- a/unsupported/Eigen/src/NonLinear/qrsolv.h +++ b/unsupported/Eigen/src/NonLinear/qrsolv.h @@ -1,8 +1,8 @@ -template + template void ei_qrsolv(int n, Scalar *r__, int ldr, - const int *ipvt, const Scalar *diag, const Scalar *qtb, Scalar *x, - Scalar *sdiag, Scalar *wa) + const int *ipvt, const Scalar *diag, const Scalar *qtb, Scalar *x, + Scalar *sdiag, Scalar *wa) { /* System generated locals */ int r_dim1, r_offset; @@ -26,141 +26,141 @@ void ei_qrsolv(int n, Scalar *r__, int ldr, /* Function Body */ -/* copy r and (q transpose)*b to preserve input and initialize s. */ -/* in particular, save the diagonal elements of r in x. */ + /* copy r and (q transpose)*b to preserve input and initialize s. */ + /* in particular, save the diagonal elements of r in x. */ for (j = 1; j <= n; ++j) { - for (i = j; i <= n; ++i) { - r__[i + j * r_dim1] = r__[j + i * r_dim1]; -/* L10: */ - } - x[j] = r__[j + j * r_dim1]; - wa[j] = qtb[j]; -/* L20: */ + for (i = j; i <= n; ++i) { + r__[i + j * r_dim1] = r__[j + i * r_dim1]; + /* L10: */ + } + x[j] = r__[j + j * r_dim1]; + wa[j] = qtb[j]; + /* L20: */ } -/* eliminate the diagonal matrix d using a givens rotation. */ + /* eliminate the diagonal matrix d using a givens rotation. */ for (j = 1; j <= n; ++j) { -/* prepare the row of d to be eliminated, locating the */ -/* diagonal element using p from the qr factorization. */ + /* prepare the row of d to be eliminated, locating the */ + /* diagonal element using p from the qr factorization. */ - l = ipvt[j]; - if (diag[l] == 0.) { - goto L90; - } - for (k = j; k <= n; ++k) { - sdiag[k] = 0.; -/* L30: */ - } - sdiag[j] = diag[l]; + l = ipvt[j]; + if (diag[l] == 0.) { + goto L90; + } + for (k = j; k <= n; ++k) { + sdiag[k] = 0.; + /* L30: */ + } + sdiag[j] = diag[l]; -/* the transformations to eliminate the row of d */ -/* modify only a single element of (q transpose)*b */ -/* beyond the first n, which is initially zero. */ + /* the transformations to eliminate the row of d */ + /* modify only a single element of (q transpose)*b */ + /* beyond the first n, which is initially zero. */ - qtbpj = 0.; - for (k = j; k <= n; ++k) { + qtbpj = 0.; + for (k = j; k <= n; ++k) { -/* determine a givens rotation which eliminates the */ -/* appropriate element in the current row of d. */ + /* determine a givens rotation which eliminates the */ + /* appropriate element in the current row of d. */ - if (sdiag[k] == 0.) - goto L70; - if ( ei_abs(r__[k + k * r_dim1]) >= ei_abs(sdiag[k])) - goto L40; - cotan = r__[k + k * r_dim1] / sdiag[k]; -/* Computing 2nd power */ - sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); - cos__ = sin__ * cotan; - goto L50; + if (sdiag[k] == 0.) + goto L70; + if ( ei_abs(r__[k + k * r_dim1]) >= ei_abs(sdiag[k])) + goto L40; + cotan = r__[k + k * r_dim1] / sdiag[k]; + /* Computing 2nd power */ + sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); + cos__ = sin__ * cotan; + goto L50; L40: - tan__ = sdiag[k] / r__[k + k * r_dim1]; -/* Computing 2nd power */ - cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); - sin__ = cos__ * tan__; + tan__ = sdiag[k] / r__[k + k * r_dim1]; + /* Computing 2nd power */ + cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); + sin__ = cos__ * tan__; L50: -/* compute the modified diagonal element of r and */ -/* the modified element of ((q transpose)*b,0). */ + /* compute the modified diagonal element of r and */ + /* the modified element of ((q transpose)*b,0). */ - r__[k + k * r_dim1] = cos__ * r__[k + k * r_dim1] + sin__ * sdiag[ - k]; - temp = cos__ * wa[k] + sin__ * qtbpj; - qtbpj = -sin__ * wa[k] + cos__ * qtbpj; - wa[k] = temp; + r__[k + k * r_dim1] = cos__ * r__[k + k * r_dim1] + sin__ * sdiag[ + k]; + temp = cos__ * wa[k] + sin__ * qtbpj; + qtbpj = -sin__ * wa[k] + cos__ * qtbpj; + wa[k] = temp; -/* accumulate the tranformation in the row of s. */ + /* accumulate the tranformation in the row of s. */ - kp1 = k + 1; - if (n < kp1) { - goto L70; - } - for (i = kp1; i <= n; ++i) { - temp = cos__ * r__[i + k * r_dim1] + sin__ * sdiag[i]; - sdiag[i] = -sin__ * r__[i + k * r_dim1] + cos__ * sdiag[ - i]; - r__[i + k * r_dim1] = temp; -/* L60: */ - } + kp1 = k + 1; + if (n < kp1) { + goto L70; + } + for (i = kp1; i <= n; ++i) { + temp = cos__ * r__[i + k * r_dim1] + sin__ * sdiag[i]; + sdiag[i] = -sin__ * r__[i + k * r_dim1] + cos__ * sdiag[ + i]; + r__[i + k * r_dim1] = temp; + /* L60: */ + } L70: -/* L80: */ - ; - } + /* L80: */ + ; + } L90: -/* store the diagonal element of s and restore */ -/* the corresponding diagonal element of r. */ + /* store the diagonal element of s and restore */ + /* the corresponding diagonal element of r. */ - sdiag[j] = r__[j + j * r_dim1]; - r__[j + j * r_dim1] = x[j]; -/* L100: */ + sdiag[j] = r__[j + j * r_dim1]; + r__[j + j * r_dim1] = x[j]; + /* L100: */ } -/* solve the triangular system for z. if the system is */ -/* singular, then obtain a least squares solution. */ + /* solve the triangular system for z. if the system is */ + /* singular, then obtain a least squares solution. */ nsing = n; for (j = 1; j <= n; ++j) { - if (sdiag[j] == 0. && nsing == n) { - nsing = j - 1; - } - if (nsing < n) { - wa[j] = 0.; - } -/* L110: */ + if (sdiag[j] == 0. && nsing == n) { + nsing = j - 1; + } + if (nsing < n) { + wa[j] = 0.; + } + /* L110: */ } if (nsing < 1) { - goto L150; + goto L150; } for (k = 1; k <= nsing; ++k) { - j = nsing - k + 1; - sum = 0.; - jp1 = j + 1; - if (nsing < jp1) { - goto L130; - } - for (i = jp1; i <= nsing; ++i) { - sum += r__[i + j * r_dim1] * wa[i]; -/* L120: */ - } + j = nsing - k + 1; + sum = 0.; + jp1 = j + 1; + if (nsing < jp1) { + goto L130; + } + for (i = jp1; i <= nsing; ++i) { + sum += r__[i + j * r_dim1] * wa[i]; + /* L120: */ + } L130: - wa[j] = (wa[j] - sum) / sdiag[j]; -/* L140: */ + wa[j] = (wa[j] - sum) / sdiag[j]; + /* L140: */ } L150: -/* permute the components of z back to components of x. */ + /* permute the components of z back to components of x. */ for (j = 1; j <= n; ++j) { - l = ipvt[j]; - x[l] = wa[j]; -/* L160: */ + l = ipvt[j]; + x[l] = wa[j]; + /* L160: */ } return; -/* last card of subroutine qrsolv. */ + /* last card of subroutine qrsolv. */ } /* qrsolv_ */ diff --git a/unsupported/Eigen/src/NonLinear/r1mpyq.h b/unsupported/Eigen/src/NonLinear/r1mpyq.h index 0ae7e0bbe..ef91da2da 100644 --- a/unsupported/Eigen/src/NonLinear/r1mpyq.h +++ b/unsupported/Eigen/src/NonLinear/r1mpyq.h @@ -1,7 +1,7 @@ template void ei_r1mpyq(int m, int n, Scalar *a, int - lda, const Scalar *v, const Scalar *w) + lda, const Scalar *v, const Scalar *w) { /* System generated locals */ int a_dim1, a_offset; @@ -19,69 +19,69 @@ void ei_r1mpyq(int m, int n, Scalar *a, int /* Function Body */ -/* apply the first set of givens rotations to a. */ + /* apply the first set of givens rotations to a. */ nm1 = n - 1; if (nm1 < 1) { - /* goto L50; */ + /* goto L50; */ return; } for (nmj = 1; nmj <= nm1; ++nmj) { - j = n - nmj; - if (ei_abs(v[j]) > 1.) { - cos__ = 1. / v[j]; - } - if (ei_abs(v[j]) > 1.) { -/* Computing 2nd power */ - sin__ = ei_sqrt(1. - ei_abs2(cos__)); - } - if (ei_abs(v[j]) <= 1.) { - sin__ = v[j]; - } - if (ei_abs(v[j]) <= 1.) { -/* Computing 2nd power */ - cos__ = ei_sqrt(1. - ei_abs2(sin__)); - } - for (i = 1; i <= m; ++i) { - temp = cos__ * a[i + j * a_dim1] - sin__ * a[i + n * a_dim1]; - a[i + n * a_dim1] = sin__ * a[i + j * a_dim1] + cos__ * a[ - i + n * a_dim1]; - a[i + j * a_dim1] = temp; -/* L10: */ - } -/* L20: */ + j = n - nmj; + if (ei_abs(v[j]) > 1.) { + cos__ = 1. / v[j]; + } + if (ei_abs(v[j]) > 1.) { + /* Computing 2nd power */ + sin__ = ei_sqrt(1. - ei_abs2(cos__)); + } + if (ei_abs(v[j]) <= 1.) { + sin__ = v[j]; + } + if (ei_abs(v[j]) <= 1.) { + /* Computing 2nd power */ + cos__ = ei_sqrt(1. - ei_abs2(sin__)); + } + for (i = 1; i <= m; ++i) { + temp = cos__ * a[i + j * a_dim1] - sin__ * a[i + n * a_dim1]; + a[i + n * a_dim1] = sin__ * a[i + j * a_dim1] + cos__ * a[ + i + n * a_dim1]; + a[i + j * a_dim1] = temp; + /* L10: */ + } + /* L20: */ } -/* apply the second set of givens rotations to a. */ + /* apply the second set of givens rotations to a. */ for (j = 1; j <= nm1; ++j) { - if (ei_abs(w[j]) > 1.) { - cos__ = 1. / w[j]; - } - if (ei_abs(w[j]) > 1.) { -/* Computing 2nd power */ - sin__ = ei_sqrt(1. - ei_abs2(cos__)); - } - if (ei_abs(w[j]) <= 1.) { - sin__ = w[j]; - } - if (ei_abs(w[j]) <= 1.) { -/* Computing 2nd power */ - cos__ = ei_sqrt(1. - ei_abs2(sin__)); - } - for (i = 1; i <= m; ++i) { - temp = cos__ * a[i + j * a_dim1] + sin__ * a[i + n * a_dim1]; - a[i + n * a_dim1] = -sin__ * a[i + j * a_dim1] + cos__ * a[ - i + n * a_dim1]; - a[i + j * a_dim1] = temp; -/* L30: */ - } -/* L40: */ + if (ei_abs(w[j]) > 1.) { + cos__ = 1. / w[j]; + } + if (ei_abs(w[j]) > 1.) { + /* Computing 2nd power */ + sin__ = ei_sqrt(1. - ei_abs2(cos__)); + } + if (ei_abs(w[j]) <= 1.) { + sin__ = w[j]; + } + if (ei_abs(w[j]) <= 1.) { + /* Computing 2nd power */ + cos__ = ei_sqrt(1. - ei_abs2(sin__)); + } + for (i = 1; i <= m; ++i) { + temp = cos__ * a[i + j * a_dim1] + sin__ * a[i + n * a_dim1]; + a[i + n * a_dim1] = -sin__ * a[i + j * a_dim1] + cos__ * a[ + i + n * a_dim1]; + a[i + j * a_dim1] = temp; + /* L30: */ + } + /* L40: */ } -/* L50: */ + /* L50: */ return; -/* last card of subroutine r1mpyq. */ + /* last card of subroutine r1mpyq. */ } /* r1mpyq_ */ diff --git a/unsupported/Eigen/src/NonLinear/r1updt.h b/unsupported/Eigen/src/NonLinear/r1updt.h index b3ae90bd8..f2ddb189b 100644 --- a/unsupported/Eigen/src/NonLinear/r1updt.h +++ b/unsupported/Eigen/src/NonLinear/r1updt.h @@ -1,5 +1,5 @@ -template + template void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v, Scalar *w, bool *sing) { /* Local variables */ @@ -17,159 +17,159 @@ void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v /* Function Body */ const Scalar giant = std::numeric_limits::max(); -/* initialize the diagonal element pointer. */ + /* 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. */ + /* 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: */ + 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. */ + /* 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; + goto L70; } for (nmj = 1; nmj <= nm1; ++nmj) { - j = n - nmj; - jj -= m - j + 1; - w[j] = 0.; - if (v[j] == 0.) { - goto L50; - } + 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. */ + /* 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; + 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__; + 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. */ + /* 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; + v[n] = sin__ * v[j] + cos__ * v[n]; + v[j] = tau; -/* apply the transformation to s and extend the spike in w. */ + /* 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: */ - } + 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: */ - ; + /* L60: */ + ; } L70: -/* add the spike from the rank 1 update to w. */ + /* add the spike from the rank 1 update to w. */ for (i = 1; i <= m; ++i) { - w[i] += v[n] * u[i]; -/* L80: */ + w[i] += v[n] * u[i]; + /* L80: */ } -/* eliminate the spike. */ + /* eliminate the spike. */ *sing = false; if (nm1 < 1) { - goto L140; + goto L140; } for (j = 1; j <= nm1; ++j) { - if (w[j] == 0.) { - goto L120; - } + if (w[j] == 0.) { + goto L120; + } -/* determine a givens rotation which eliminates the */ -/* j-th element of the spike. */ + /* 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; + 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__; + 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. */ + /* 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: */ - } + 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. */ + /* store the information necessary to recover the */ + /* givens rotation. */ - w[j] = tau; + w[j] = tau; L120: -/* test for zero diagonal elements in the output s. */ + /* test for zero diagonal elements in the output s. */ - if (s[jj] == 0.) { - *sing = true; - } - jj += m - j + 1; -/* L130: */ + if (s[jj] == 0.) { + *sing = true; + } + jj += m - j + 1; + /* L130: */ } L140: -/* move w back into the last column of the output s. */ + /* 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: */ + s[l] = w[i]; + ++l; + /* L150: */ } if (s[jj] == 0.) { - *sing = true; + *sing = true; } return; -/* last card of subroutine r1updt. */ + /* last card of subroutine r1updt. */ } /* r1updt_ */