GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 11:49:02 +08:00
Code Issues Packages Projects Releases Wiki Activity

fix indentation (and only that)

Browse Source
This commit is contained in:
Thomas Capricelli 2009-09-14 23:47:44 +02:00
parent ab88ba6f7f
commit 6d8baa757e
5 changed files with 377 additions and 377 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

88
unsupported/Eigen/src/NonLinear/qform.h
View File

@ -1,7 +1,7 @@
template <typename Scalar> template <typename Scalar>
void ei_qform(int m, int n, Scalar *q, int void ei_qform(int m, int n, Scalar *q, int
ldq, Scalar *wa) ldq, Scalar *wa)
{ {
/* System generated locals */ /* System generated locals */
int q_dim1, q_offset; int q_dim1, q_offset;
@ -19,71 +19,71 @@ void ei_qform(int m, int n, Scalar *q, int
/* Function Body */ /* 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); minmn = std::min(m,n);
if (minmn < 2) { if (minmn < 2) {
goto L30; goto L30;
} }
for (j = 2; j <= minmn; ++j) { for (j = 2; j <= minmn; ++j) {
jm1 = j - 1; jm1 = j - 1;
for (i = 1; i <= jm1; ++i) { for (i = 1; i <= jm1; ++i) {
q[i + j * q_dim1] = 0.; q[i + j * q_dim1] = 0.;
/* L10: */ /* L10: */
} }
/* L20: */ /* L20: */
} }
L30: L30:
/* initialize remaining columns to those of the identity matrix. */ /* initialize remaining columns to those of the identity matrix. */
np1 = n + 1; np1 = n + 1;
if (m < np1) { if (m < np1) {
goto L60; goto L60;
} }
for (j = np1; j <= m; ++j) { for (j = np1; j <= m; ++j) {
for (i = 1; i <= m; ++i) { for (i = 1; i <= m; ++i) {
q[i + j * q_dim1] = 0.; q[i + j * q_dim1] = 0.;
/* L40: */ /* L40: */
} }
q[j + j * q_dim1] = 1.; q[j + j * q_dim1] = 1.;
/* L50: */ /* L50: */
} }
L60: L60:
/* accumulate q from its factored form. */ /* accumulate q from its factored form. */
for (l = 1; l <= minmn; ++l) { for (l = 1; l <= minmn; ++l) {
k = minmn - l + 1; k = minmn - l + 1;
for (i = k; i <= m; ++i) { for (i = k; i <= m; ++i) {
wa[i] = q[i + k * q_dim1]; wa[i] = q[i + k * q_dim1];
q[i + k * q_dim1] = 0.; q[i + k * q_dim1] = 0.;
/* L70: */ /* L70: */
} }
q[k + k * q_dim1] = 1.; q[k + k * q_dim1] = 1.;
if (wa[k] == 0.) { if (wa[k] == 0.) {
goto L110; goto L110;
} }
for (j = k; j <= m; ++j) { for (j = k; j <= m; ++j) {
sum = 0.; sum = 0.;
for (i = k; i <= m; ++i) { for (i = k; i <= m; ++i) {
sum += q[i + j * q_dim1] * wa[i]; sum += q[i + j * q_dim1] * wa[i];
/* L80: */ /* L80: */
} }
temp = sum / wa[k]; temp = sum / wa[k];
for (i = k; i <= m; ++i) { for (i = k; i <= m; ++i) {
q[i + j * q_dim1] -= temp * wa[i]; q[i + j * q_dim1] -= temp * wa[i];
/* L90: */ /* L90: */
} }
/* L100: */ /* L100: */
} }
L110: L110:
/* L120: */ /* L120: */
; ;
} }
return; return;
/* last card of subroutine qform. */ /* last card of subroutine qform. */
} /* qform_ */ } /* qform_ */

172
unsupported/Eigen/src/NonLinear/qrfac.h
View File

@ -1,8 +1,8 @@
template <typename Scalar> template <typename Scalar>
void ei_qrfac(int m, int n, Scalar *a, int void ei_qrfac(int m, int n, Scalar *a, int
lda, int pivot, int *ipvt, int /* lipvt */, Scalar *rdiag, lda, int pivot, int *ipvt, int /* lipvt */, Scalar *rdiag,
Scalar *acnorm) Scalar *acnorm)
{ {
/* System generated locals */ /* System generated locals */
int a_dim1, a_offset; int a_dim1, a_offset;
@ -28,109 +28,109 @@ void ei_qrfac(int m, int n, Scalar *a, int
/* Function Body */ /* Function Body */
const Scalar epsmch = epsilon<Scalar>(); const Scalar epsmch = epsilon<Scalar>();
/* compute the initial column norms and initialize several arrays. */ /* compute the initial column norms and initialize several arrays. */
for (j = 1; j <= n; ++j) { for (j = 1; j <= n; ++j) {
acnorm[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm(); acnorm[j] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j * a_dim1 + 1],m).blueNorm();
rdiag[j] = acnorm[j]; rdiag[j] = acnorm[j];
wa[j] = rdiag[j]; wa[j] = rdiag[j];
if (pivot) { if (pivot) {
ipvt[j] = j; ipvt[j] = j;
} }
/* L10: */ /* L10: */
} }
/* reduce a to r with householder transformations. */ /* reduce a to r with householder transformations. */
minmn = std::min(m,n); minmn = std::min(m,n);
for (j = 1; j <= minmn; ++j) { for (j = 1; j <= minmn; ++j) {
if (! (pivot)) { if (! (pivot)) {
goto L40; goto L40;
} }
/* bring the column of largest norm into the pivot position. */ /* bring the column of largest norm into the pivot position. */
kmax = j; kmax = j;
for (k = j; k <= n; ++k) { for (k = j; k <= n; ++k) {
if (rdiag[k] > rdiag[kmax]) { if (rdiag[k] > rdiag[kmax]) {
kmax = k; kmax = k;
} }
/* L20: */ /* L20: */
} }
if (kmax == j) { if (kmax == j) {
goto L40; goto L40;
} }
for (i = 1; i <= m; ++i) { for (i = 1; i <= m; ++i) {
temp = a[i + j * a_dim1]; temp = a[i + j * a_dim1];
a[i + j * a_dim1] = a[i + kmax * a_dim1]; a[i + j * a_dim1] = a[i + kmax * a_dim1];
a[i + kmax * a_dim1] = temp; a[i + kmax * a_dim1] = temp;
/* L30: */ /* L30: */
} }
rdiag[kmax] = rdiag[j]; rdiag[kmax] = rdiag[j];
wa[kmax] = wa[j]; wa[kmax] = wa[j];
k = ipvt[j]; k = ipvt[j];
ipvt[j] = ipvt[kmax]; ipvt[j] = ipvt[kmax];
ipvt[kmax] = k; ipvt[kmax] = k;
L40: L40:
/* compute the householder transformation to reduce the */ /* compute the householder transformation to reduce the */
/* j-th column of a to a multiple of the j-th unit vector. */ /* 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(); ajnorm = Map< Matrix< Scalar, Dynamic, 1 > >(&a[j + j * a_dim1],m-j+1).blueNorm();
if (ajnorm == 0.) { if (ajnorm == 0.) {
goto L100; goto L100;
} }
if (a[j + j * a_dim1] < 0.) { if (a[j + j * a_dim1] < 0.) {
ajnorm = -ajnorm; ajnorm = -ajnorm;
} }
for (i = j; i <= m; ++i) { for (i = j; i <= m; ++i) {
a[i + j * a_dim1] /= ajnorm; a[i + j * a_dim1] /= ajnorm;
/* L50: */ /* L50: */
} }
a[j + j * a_dim1] += 1.; a[j + j * a_dim1] += 1.;
/* apply the transformation to the remaining columns */ /* apply the transformation to the remaining columns */
/* and update the norms. */ /* and update the norms. */
jp1 = j + 1; jp1 = j + 1;
if (n < jp1) { if (n < jp1) {
goto L100; goto L100;
} }
for (k = jp1; k <= n; ++k) { for (k = jp1; k <= n; ++k) {
sum = 0.; sum = 0.;
for (i = j; i <= m; ++i) { for (i = j; i <= m; ++i) {
sum += a[i + j * a_dim1] * a[i + k * a_dim1]; sum += a[i + j * a_dim1] * a[i + k * a_dim1];
/* L60: */ /* L60: */
} }
temp = sum / a[j + j * a_dim1]; temp = sum / a[j + j * a_dim1];
for (i = j; i <= m; ++i) { for (i = j; i <= m; ++i) {
a[i + k * a_dim1] -= temp * a[i + j * a_dim1]; a[i + k * a_dim1] -= temp * a[i + j * a_dim1];
/* L70: */ /* L70: */
} }
if (! (pivot) || rdiag[k] == 0.) { if (! (pivot) || rdiag[k] == 0.) {
goto L80; goto L80;
} }
temp = a[j + k * a_dim1] / rdiag[k]; temp = a[j + k * a_dim1] / rdiag[k];
/* Computing MAX */ /* Computing MAX */
/* Computing 2nd power */ /* Computing 2nd power */
rdiag[k] *= ei_sqrt((std::max(Scalar(0.), Scalar(1.)-ei_abs2(temp)))); rdiag[k] *= ei_sqrt((std::max(Scalar(0.), Scalar(1.)-ei_abs2(temp))));
/* Computing 2nd power */ /* Computing 2nd power */
if (Scalar(.05) * ei_abs2(rdiag[k] / wa[k]) > epsmch) { if (Scalar(.05) * ei_abs2(rdiag[k] / wa[k]) > epsmch) {
goto L80; goto L80;
} }
rdiag[k] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[jp1 + k * a_dim1],m-j).blueNorm(); rdiag[k] = Map< Matrix< Scalar, Dynamic, 1 > >(&a[jp1 + k * a_dim1],m-j).blueNorm();
wa[k] = rdiag[k]; wa[k] = rdiag[k];
L80: L80:
/* L90: */ /* L90: */
; ;
} }
L100: L100:
rdiag[j] = -ajnorm; rdiag[j] = -ajnorm;
/* L110: */ /* L110: */
} }
return; return;
/* last card of subroutine qrfac. */ /* last card of subroutine qrfac. */
} /* qrfac_ */ } /* qrfac_ */

196
unsupported/Eigen/src/NonLinear/qrsolv.h
View File

@ -1,8 +1,8 @@
template <typename Scalar> template <typename Scalar>
void ei_qrsolv(int n, Scalar *r__, int ldr, void ei_qrsolv(int n, Scalar *r__, int ldr,
const int *ipvt, const Scalar *diag, const Scalar *qtb, Scalar *x, const int *ipvt, const Scalar *diag, const Scalar *qtb, Scalar *x,
Scalar *sdiag, Scalar *wa) Scalar *sdiag, Scalar *wa)
{ {
/* System generated locals */ /* System generated locals */
int r_dim1, r_offset; int r_dim1, r_offset;
@ -26,141 +26,141 @@ void ei_qrsolv(int n, Scalar *r__, int ldr,
/* Function Body */ /* Function Body */
/* copy r and (q transpose)*b to preserve input and initialize s. */ /* copy r and (q transpose)*b to preserve input and initialize s. */
/* in particular, save the diagonal elements of r in x. */ /* in particular, save the diagonal elements of r in x. */
for (j = 1; j <= n; ++j) { for (j = 1; j <= n; ++j) {
for (i = j; i <= n; ++i) { for (i = j; i <= n; ++i) {
r__[i + j * r_dim1] = r__[j + i * r_dim1]; r__[i + j * r_dim1] = r__[j + i * r_dim1];
/* L10: */ /* L10: */
} }
x[j] = r__[j + j * r_dim1]; x[j] = r__[j + j * r_dim1];
wa[j] = qtb[j]; wa[j] = qtb[j];
/* L20: */ /* 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) { for (j = 1; j <= n; ++j) {
/* prepare the row of d to be eliminated, locating the */ /* prepare the row of d to be eliminated, locating the */
/* diagonal element using p from the qr factorization. */ /* diagonal element using p from the qr factorization. */
l = ipvt[j]; l = ipvt[j];
if (diag[l] == 0.) { if (diag[l] == 0.) {
goto L90; goto L90;
} }
for (k = j; k <= n; ++k) { for (k = j; k <= n; ++k) {
sdiag[k] = 0.; sdiag[k] = 0.;
/* L30: */ /* L30: */
} }
sdiag[j] = diag[l]; sdiag[j] = diag[l];
/* the transformations to eliminate the row of d */ /* the transformations to eliminate the row of d */
/* modify only a single element of (q transpose)*b */ /* modify only a single element of (q transpose)*b */
/* beyond the first n, which is initially zero. */ /* beyond the first n, which is initially zero. */
qtbpj = 0.; qtbpj = 0.;
for (k = j; k <= n; ++k) { for (k = j; k <= n; ++k) {
/* determine a givens rotation which eliminates the */ /* determine a givens rotation which eliminates the */
/* appropriate element in the current row of d. */ /* appropriate element in the current row of d. */
if (sdiag[k] == 0.) if (sdiag[k] == 0.)
goto L70; goto L70;
if ( ei_abs(r__[k + k * r_dim1]) >= ei_abs(sdiag[k])) if ( ei_abs(r__[k + k * r_dim1]) >= ei_abs(sdiag[k]))
goto L40; goto L40;
cotan = r__[k + k * r_dim1] / sdiag[k]; cotan = r__[k + k * r_dim1] / sdiag[k];
/* Computing 2nd power */ /* Computing 2nd power */
sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
cos__ = sin__ * cotan; cos__ = sin__ * cotan;
goto L50; goto L50;
L40: L40:
tan__ = sdiag[k] / r__[k + k * r_dim1]; tan__ = sdiag[k] / r__[k + k * r_dim1];
/* Computing 2nd power */ /* Computing 2nd power */
cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
sin__ = cos__ * tan__; sin__ = cos__ * tan__;
L50: L50:
/* compute the modified diagonal element of r and */ /* compute the modified diagonal element of r and */
/* the modified element of ((q transpose)*b,0). */ /* the modified element of ((q transpose)*b,0). */
r__[k + k * r_dim1] = cos__ * r__[k + k * r_dim1] + sin__ * sdiag[ r__[k + k * r_dim1] = cos__ * r__[k + k * r_dim1] + sin__ * sdiag[
k]; k];
temp = cos__ * wa[k] + sin__ * qtbpj; temp = cos__ * wa[k] + sin__ * qtbpj;
qtbpj = -sin__ * wa[k] + cos__ * qtbpj; qtbpj = -sin__ * wa[k] + cos__ * qtbpj;
wa[k] = temp; wa[k] = temp;
/* accumulate the tranformation in the row of s. */ /* accumulate the tranformation in the row of s. */
kp1 = k + 1; kp1 = k + 1;
if (n < kp1) { if (n < kp1) {
goto L70; goto L70;
} }
for (i = kp1; i <= n; ++i) { for (i = kp1; i <= n; ++i) {
temp = cos__ * r__[i + k * r_dim1] + sin__ * sdiag[i]; temp = cos__ * r__[i + k * r_dim1] + sin__ * sdiag[i];
sdiag[i] = -sin__ * r__[i + k * r_dim1] + cos__ * sdiag[ sdiag[i] = -sin__ * r__[i + k * r_dim1] + cos__ * sdiag[
i]; i];
r__[i + k * r_dim1] = temp; r__[i + k * r_dim1] = temp;
/* L60: */ /* L60: */
} }
L70: L70:
/* L80: */ /* L80: */
; ;
} }
L90: L90:
/* store the diagonal element of s and restore */ /* store the diagonal element of s and restore */
/* the corresponding diagonal element of r. */ /* the corresponding diagonal element of r. */
sdiag[j] = r__[j + j * r_dim1]; sdiag[j] = r__[j + j * r_dim1];
r__[j + j * r_dim1] = x[j]; r__[j + j * r_dim1] = x[j];
/* L100: */ /* L100: */
} }
/* solve the triangular system for z. if the system is */ /* solve the triangular system for z. if the system is */
/* singular, then obtain a least squares solution. */ /* singular, then obtain a least squares solution. */
nsing = n; nsing = n;
for (j = 1; j <= n; ++j) { for (j = 1; j <= n; ++j) {
if (sdiag[j] == 0. && nsing == n) { if (sdiag[j] == 0. && nsing == n) {
nsing = j - 1; nsing = j - 1;
} }
if (nsing < n) { if (nsing < n) {
wa[j] = 0.; wa[j] = 0.;
} }
/* L110: */ /* L110: */
} }
if (nsing < 1) { if (nsing < 1) {
goto L150; goto L150;
} }
for (k = 1; k <= nsing; ++k) { for (k = 1; k <= nsing; ++k) {
j = nsing - k + 1; j = nsing - k + 1;
sum = 0.; sum = 0.;
jp1 = j + 1; jp1 = j + 1;
if (nsing < jp1) { if (nsing < jp1) {
goto L130; goto L130;
} }
for (i = jp1; i <= nsing; ++i) { for (i = jp1; i <= nsing; ++i) {
sum += r__[i + j * r_dim1] * wa[i]; sum += r__[i + j * r_dim1] * wa[i];
/* L120: */ /* L120: */
} }
L130: L130:
wa[j] = (wa[j] - sum) / sdiag[j]; wa[j] = (wa[j] - sum) / sdiag[j];
/* L140: */ /* L140: */
} }
L150: 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) { for (j = 1; j <= n; ++j) {
l = ipvt[j]; l = ipvt[j];
x[l] = wa[j]; x[l] = wa[j];
/* L160: */ /* L160: */
} }
return; return;
/* last card of subroutine qrsolv. */ /* last card of subroutine qrsolv. */
} /* qrsolv_ */ } /* qrsolv_ */

102
unsupported/Eigen/src/NonLinear/r1mpyq.h
View File

@ -1,7 +1,7 @@
template <typename Scalar> template <typename Scalar>
void ei_r1mpyq(int m, int n, Scalar *a, int 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 */ /* System generated locals */
int a_dim1, a_offset; int a_dim1, a_offset;
@ -19,69 +19,69 @@ void ei_r1mpyq(int m, int n, Scalar *a, int
/* Function Body */ /* Function Body */
/* apply the first set of givens rotations to a. */ /* apply the first set of givens rotations to a. */
nm1 = n - 1; nm1 = n - 1;
if (nm1 < 1) { if (nm1 < 1) {
/* goto L50; */ /* goto L50; */
return; return;
} }
for (nmj = 1; nmj <= nm1; ++nmj) { for (nmj = 1; nmj <= nm1; ++nmj) {
j = n - nmj; j = n - nmj;
if (ei_abs(v[j]) > 1.) { if (ei_abs(v[j]) > 1.) {
cos__ = 1. / v[j]; cos__ = 1. / v[j];
} }
if (ei_abs(v[j]) > 1.) { if (ei_abs(v[j]) > 1.) {
/* Computing 2nd power */ /* Computing 2nd power */
sin__ = ei_sqrt(1. - ei_abs2(cos__)); sin__ = ei_sqrt(1. - ei_abs2(cos__));
} }
if (ei_abs(v[j]) <= 1.) { if (ei_abs(v[j]) <= 1.) {
sin__ = v[j]; sin__ = v[j];
} }
if (ei_abs(v[j]) <= 1.) { if (ei_abs(v[j]) <= 1.) {
/* Computing 2nd power */ /* Computing 2nd power */
cos__ = ei_sqrt(1. - ei_abs2(sin__)); cos__ = ei_sqrt(1. - ei_abs2(sin__));
} }
for (i = 1; i <= m; ++i) { for (i = 1; i <= m; ++i) {
temp = cos__ * a[i + j * a_dim1] - sin__ * a[i + n * a_dim1]; 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[ a[i + n * a_dim1] = sin__ * a[i + j * a_dim1] + cos__ * a[
i + n * a_dim1]; i + n * a_dim1];
a[i + j * a_dim1] = temp; a[i + j * a_dim1] = temp;
/* L10: */ /* L10: */
} }
/* L20: */ /* 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) { for (j = 1; j <= nm1; ++j) {
if (ei_abs(w[j]) > 1.) { if (ei_abs(w[j]) > 1.) {
cos__ = 1. / w[j]; cos__ = 1. / w[j];
} }
if (ei_abs(w[j]) > 1.) { if (ei_abs(w[j]) > 1.) {
/* Computing 2nd power */ /* Computing 2nd power */
sin__ = ei_sqrt(1. - ei_abs2(cos__)); sin__ = ei_sqrt(1. - ei_abs2(cos__));
} }
if (ei_abs(w[j]) <= 1.) { if (ei_abs(w[j]) <= 1.) {
sin__ = w[j]; sin__ = w[j];
} }
if (ei_abs(w[j]) <= 1.) { if (ei_abs(w[j]) <= 1.) {
/* Computing 2nd power */ /* Computing 2nd power */
cos__ = ei_sqrt(1. - ei_abs2(sin__)); cos__ = ei_sqrt(1. - ei_abs2(sin__));
} }
for (i = 1; i <= m; ++i) { for (i = 1; i <= m; ++i) {
temp = cos__ * a[i + j * a_dim1] + sin__ * a[i + n * a_dim1]; 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[ a[i + n * a_dim1] = -sin__ * a[i + j * a_dim1] + cos__ * a[
i + n * a_dim1]; i + n * a_dim1];
a[i + j * a_dim1] = temp; a[i + j * a_dim1] = temp;
/* L30: */ /* L30: */
} }
/* L40: */ /* L40: */
} }
/* L50: */ /* L50: */
return; return;
/* last card of subroutine r1mpyq. */ /* last card of subroutine r1mpyq. */
} /* r1mpyq_ */ } /* r1mpyq_ */

196
unsupported/Eigen/src/NonLinear/r1updt.h
View File

@ -1,5 +1,5 @@
template <typename Scalar> template <typename Scalar>
void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v, Scalar *w, bool *sing) void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v, Scalar *w, bool *sing)
{ {
/* Local variables */ /* Local variables */
@ -17,159 +17,159 @@ void ei_r1updt(int m, int n, Scalar *s, int /* ls */, const Scalar *u, Scalar *v
/* Function Body */ /* Function Body */
const Scalar giant = std::numeric_limits<Scalar>::max(); const Scalar giant = std::numeric_limits<Scalar>::max();
/* initialize the diagonal element pointer. */ /* initialize the diagonal element pointer. */
jj = n * ((m << 1) - n + 1) / 2 - (m - n); 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; l = jj;
for (i = n; i <= m; ++i) { for (i = n; i <= m; ++i) {
w[i] = s[l]; w[i] = s[l];
++l; ++l;
/* L10: */ /* L10: */
} }
/* rotate the vector v into a multiple of the n-th unit vector */ /* rotate the vector v into a multiple of the n-th unit vector */
/* in such a way that a spike is introduced into w. */ /* in such a way that a spike is introduced into w. */
nm1 = n - 1; nm1 = n - 1;
if (nm1 < 1) { if (nm1 < 1) {
goto L70; goto L70;
} }
for (nmj = 1; nmj <= nm1; ++nmj) { for (nmj = 1; nmj <= nm1; ++nmj) {
j = n - nmj; j = n - nmj;
jj -= m - j + 1; jj -= m - j + 1;
w[j] = 0.; w[j] = 0.;
if (v[j] == 0.) { if (v[j] == 0.) {
goto L50; goto L50;
} }
/* determine a givens rotation which eliminates the */ /* determine a givens rotation which eliminates the */
/* j-th element of v. */ /* j-th element of v. */
if (ei_abs(v[n]) >= ei_abs(v[j])) if (ei_abs(v[n]) >= ei_abs(v[j]))
goto L20; goto L20;
cotan = v[n] / v[j]; cotan = v[n] / v[j];
/* Computing 2nd power */ /* Computing 2nd power */
sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
cos__ = sin__ * cotan; cos__ = sin__ * cotan;
tau = 1.; tau = 1.;
if (ei_abs(cos__) * giant > 1.) { if (ei_abs(cos__) * giant > 1.) {
tau = 1. / cos__; tau = 1. / cos__;
} }
goto L30; goto L30;
L20: L20:
tan__ = v[j] / v[n]; tan__ = v[j] / v[n];
/* Computing 2nd power */ /* Computing 2nd power */
cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
sin__ = cos__ * tan__; sin__ = cos__ * tan__;
tau = sin__; tau = sin__;
L30: L30:
/* apply the transformation to v and store the information */ /* apply the transformation to v and store the information */
/* necessary to recover the givens rotation. */ /* necessary to recover the givens rotation. */
v[n] = sin__ * v[j] + cos__ * v[n]; v[n] = sin__ * v[j] + cos__ * v[n];
v[j] = tau; 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; l = jj;
for (i = j; i <= m; ++i) { for (i = j; i <= m; ++i) {
temp = cos__ * s[l] - sin__ * w[i]; temp = cos__ * s[l] - sin__ * w[i];
w[i] = sin__ * s[l] + cos__ * w[i]; w[i] = sin__ * s[l] + cos__ * w[i];
s[l] = temp; s[l] = temp;
++l; ++l;
/* L40: */ /* L40: */
} }
L50: L50:
/* L60: */ /* L60: */
; ;
} }
L70: 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) { for (i = 1; i <= m; ++i) {
w[i] += v[n] * u[i]; w[i] += v[n] * u[i];
/* L80: */ /* L80: */
} }
/* eliminate the spike. */ /* eliminate the spike. */
*sing = false; *sing = false;
if (nm1 < 1) { if (nm1 < 1) {
goto L140; goto L140;
} }
for (j = 1; j <= nm1; ++j) { for (j = 1; j <= nm1; ++j) {
if (w[j] == 0.) { if (w[j] == 0.) {
goto L120; goto L120;
} }
/* determine a givens rotation which eliminates the */ /* determine a givens rotation which eliminates the */
/* j-th element of the spike. */ /* j-th element of the spike. */
if (ei_abs(s[jj]) >= ei_abs(w[j])) if (ei_abs(s[jj]) >= ei_abs(w[j]))
goto L90; goto L90;
cotan = s[jj] / w[j]; cotan = s[jj] / w[j];
/* Computing 2nd power */ /* Computing 2nd power */
sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan)); sin__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(cotan));
cos__ = sin__ * cotan; cos__ = sin__ * cotan;
tau = 1.; tau = 1.;
if (ei_abs(cos__) * giant > 1.) { if (ei_abs(cos__) * giant > 1.) {
tau = 1. / cos__; tau = 1. / cos__;
} }
goto L100; goto L100;
L90: L90:
tan__ = w[j] / s[jj]; tan__ = w[j] / s[jj];
/* Computing 2nd power */ /* Computing 2nd power */
cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__)); cos__ = Scalar(.5) / ei_sqrt(Scalar(0.25) + Scalar(0.25) * ei_abs2(tan__));
sin__ = cos__ * tan__; sin__ = cos__ * tan__;
tau = sin__; tau = sin__;
L100: 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; l = jj;
for (i = j; i <= m; ++i) { for (i = j; i <= m; ++i) {
temp = cos__ * s[l] + sin__ * w[i]; temp = cos__ * s[l] + sin__ * w[i];
w[i] = -sin__ * s[l] + cos__ * w[i]; w[i] = -sin__ * s[l] + cos__ * w[i];
s[l] = temp; s[l] = temp;
++l; ++l;
/* L110: */ /* L110: */
} }
/* store the information necessary to recover the */ /* store the information necessary to recover the */
/* givens rotation. */ /* givens rotation. */
w[j] = tau; w[j] = tau;
L120: L120:
/* test for zero diagonal elements in the output s. */ /* test for zero diagonal elements in the output s. */
if (s[jj] == 0.) { if (s[jj] == 0.) {
*sing = true; *sing = true;
} }
jj += m - j + 1; jj += m - j + 1;
/* L130: */ /* L130: */
} }
L140: 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; l = jj;
for (i = n; i <= m; ++i) { for (i = n; i <= m; ++i) {
s[l] = w[i]; s[l] = w[i];
++l; ++l;
/* L150: */ /* L150: */
} }
if (s[jj] == 0.) { if (s[jj] == 0.) {
*sing = true; *sing = true;
} }
return; return;
/* last card of subroutine r1updt. */ /* last card of subroutine r1updt. */
} /* r1updt_ */ } /* r1updt_ */