#include /* flowkm.f -- translated by f2c (version of 3 February 1990 3:36:42). You must link the resulting object file with the libraries: -lF77 -lI77 -lm -lc (in that order) */ #include "f2c.h" #include "autlim.h" /* Table of constant values */ static integer c__1 = 1; static integer c__50 = NAUTO; static doublereal c_b24 = 1e-16; static doublereal c_b30 = 1.; static doublereal c_b32 = 0.; static integer c__0 = 0; static integer c__9 = 9; static doublereal c_b338 = -1.; static integer c__3 = 3; static integer c__5 = 5; /* Routines included in this file: */ /* subroutine flowkm : new routine to compute floquet multipliers */ /* subroutine dhhpr : compute a Householder matrix */ /* subroutine dhhap : appy a Householder matrix */ /* subroutine qzhes : QZ reduction to Hessenberg form (from EISPAC)*/ /* subroutine qzit : QZ reduction to quasi-upper triangular form (from EISPAC)*/ /* subroutine qzval : QZ calculation of eigenvalues (from EISPAC)*/ /* subroutine dscal : scale a vector by a constant (from BLAS-1)*/ /* subroutine dgemm : matrix-matrix multiply (from BLAS-3)*/ /* subroutine xerbla : BLAS error handling routine (from BLAS-2)*/ /* subroutine daxpy : constant times a vector plus a vector (from BLAS-1)*/ /* subroutine drot : apply a plane rotation (from BLAS-1)*/ /* subroutine dswap : swap two vectors (from BLAS-1)*/ /* subroutine dgemc : matrix-matrix copy */ /* subroutines ezsvd, ndrotg, ndsvd, prse, sig22, sigmin, sndrtg : */ /* Demmel-Kahan svd routines */ /* function epslon : machine constant routine (from EISPAC)*/ /* function dnrm2 : compute l2-norm of a vector (from BLAS-1)*/ /* function ddot : dot product of two vectors (from BLAS-1)*/ /* function idamax : find index of element with max abs value (from BLAS-1)*/ /* function lsame : compare character strings (from BLAS-2)*/ /* Subroutine */ int flowkm_(ndim, c0, c1, iid, ir, ic, rwork, ev) integer *ndim; doublereal *c0, *c1; integer *iid, *ir, *ic; doublereal *rwork; doublecomplex *ev; { /* Format strings */ static char fmt_900[] = "(\002 *** SUBROUTINE FLOWKM DOESN'T HAVE ENOUGH\ STORAGE ***\002,/,\002 *** RECOMPILE WITH LARGER VALUE OF MAXDIM **\ *\002)"; static char fmt_101[] = "(\002 UNDEFLATED CIRCUIT PENCIL (C0, C1) \002)"; static char fmt_102[] = "(\002 C0 : \002)"; static char fmt_104[] = "(1x,6e23.16)"; static char fmt_103[] = "(\002 C1 : \002)"; static char fmt_901[] = "(\002 *** WARNING FROM SUBROUTINE FLOWKM ***\ \002,/,\002 *** SVD routine returned with SVDINF = \002,i4,\002 ***\002,/\ ,\002 *** Floquet mulitplier calculations may be wrong ***\002)"; static char fmt_105[] = "(\002 DEFLATED CIRCUIT PENCIL (H2^T)*(C0, C1)*(\ H1) \002)"; static char fmt_106[] = "(\002 (H2^T)*C0*(H1) : \002)"; static char fmt_107[] = "(\002 (H2^T)*C1*(H1) : \002)"; static char fmt_902[] = "(\002 *** WARNING FROM SUBROUTINE FLOWKM ***\ \002,/,\002 *** QZ routine returned with QZIERR = \002,i4,\002 ***\002,/,\ \002 *** Floquet mulitplier calculations may be wrong ***\002)"; static char fmt_903[] = "(\002 *** WARNING FROM SUBROUTINE FLOWKM ***\ \002,/,\002 NOTE: infinite Floquet multiplier represented by \002,/,\002 \ DMPLX( 1.0D+99, 1.0D+99 )\002)"; /* System generated locals */ integer c0_dim1, c0_offset, c1_dim1, c1_offset, rwork_dim1, rwork_offset, i_1, i_2; doublereal d_1, d_2, d_3, d_4; doublecomplex z_1; /* Builtin functions */ integer s_wsfe(), e_wsfe(); /* Subroutine */ int s_stop(); integer do_fio(); /* Local variables */ static doublereal beta, svde[NAUTO], svds[NAUTO+1], svdu[1] /* was [1][1] */, svdv[NxNAUTO] /* was [50][50] */; extern /* Subroutine */ int qzit_(); extern doublereal dnrm2_(); static integer i, j; extern /* Subroutine */ int dgemc_(), dhhap_(); static doublereal v[NAUTO], x[NAUTO]; extern /* Subroutine */ int dgemm_(), dhhpr_(); static logical infev; static doublereal const_; extern /* Subroutine */ int qzhes_(), ezsvd_(), qzval_(); static integer ndimm1; static doublereal nrmc0x, nrmc1x, qzalfi[NAUTO], qzbeta[NAUTO]; static integer svdinf; static doublereal qzalfr[NAUTO]; static integer qzierr; static doublereal svdwrk[NAUTO], qzz[1] /* was [1][1] */; /* Fortran I/O blocks */ static cilist io__1 = { 0, 6, 0, fmt_900, 0 }; static cilist io__2 = { 0, 9, 0, fmt_101, 0 }; static cilist io__3 = { 0, 9, 0, fmt_102, 0 }; static cilist io__5 = { 0, 9, 0, fmt_104, 0 }; static cilist io__7 = { 0, 9, 0, fmt_103, 0 }; static cilist io__8 = { 0, 9, 0, fmt_104, 0 }; static cilist io__15 = { 0, 9, 0, fmt_901, 0 }; static cilist io__22 = { 0, 9, 0, fmt_105, 0 }; static cilist io__23 = { 0, 9, 0, fmt_106, 0 }; static cilist io__24 = { 0, 9, 0, fmt_104, 0 }; static cilist io__25 = { 0, 9, 0, fmt_107, 0 }; static cilist io__26 = { 0, 9, 0, fmt_104, 0 }; static cilist io__30 = { 0, 6, 0, fmt_902, 0 }; static cilist io__35 = { 0, 6, 0, fmt_903, 0 }; /* Parameter adjustments */ c0_dim1 = *ndim; c0_offset = c0_dim1 + 1; c0 -= c0_offset; c1_dim1 = *ndim; c1_offset = c1_dim1 + 1; c1 -= c1_offset; --ir; --ic; rwork_dim1 = *ndim; rwork_offset = rwork_dim1 + 1; rwork -= rwork_offset; --ev; /* Function Body */ /* IMPLICIT UNDEFINED (A-Z,a-z) */ /* Subroutine to compute Floquet multipliers via the */ /* "deflated circuit pencil" method, which is described in */ /* T.F. Fairgrieve and A.D. Jepson, "O.K. Floquet Multipliers". */ /* Available from Tom Fairgrieve, Department of Computer Science, */ /* University of Toronto, Toronto, Ontario, CANADA M5S 1A4 */ /* Please inform Tom Fairgrieve (tff@na.utoronto.ca) of any */ /* modifications to or errors in these routines. */ /* This routine is called by a modified version of */ /* the AUTO'86 function FNSPBV. */ /* Modified by Bob Olsen (January 1991) */ /* Three unnecessary assignment statements (as indicated by the */ /* Apollo FORTRAN compiler) are commented out with c01/91 in */ /* columns 1 through 6. */ /* Parameter declarations: */ /* Local declarations: */ /* Maximum order of problem. **** (To avoid changing too many AUTO'86 */ /* routines, I need to declare some local array storage here --- MAXDIM */ /* will need to be increased if someone has increased the "20"'s in */ /* AUTO'86s DFINIT subroutine.) **** */ /* storage for SVD computations */ /* compute right singular vectors only */ /* storage for generalized eigenvalue computations */ /* don't want to accumulate the transforms --- vectors not needed */ /* BLAS routines */ /* routines from EISPAC */ /* own routines */ /* Jim Demmel's svd routine (demmel@nyu.edu) */ /* builtin F77 functions */ /* Make sure that you have enough local storage. */ if (*ndim > NAUTO) { s_wsfe(&io__1); e_wsfe(); s_stop("", 0L); } /* Print the undeflated circuit pencil (C0, C1). */ if (*iid > 2) { s_wsfe(&io__2); e_wsfe(); s_wsfe(&io__3); e_wsfe(); i_1 = *ndim; for (i = 1; i <= i_1; ++i) { s_wsfe(&io__5); i_2 = *ndim; for (j = 1; j <= i_2; ++j) { do_fio(&c__1, (char *)&c0[i + j * c0_dim1], (ftnlen)sizeof( doublereal)); } e_wsfe(); /* L1: */ } s_wsfe(&io__7); e_wsfe(); i_1 = *ndim; for (i = 1; i <= i_1; ++i) { s_wsfe(&io__8); i_2 = *ndim; for (j = 1; j <= i_2; ++j) { do_fio(&c__1, (char *)&c1[i + j * c1_dim1], (ftnlen)sizeof( doublereal)); } e_wsfe(); /* L2: */ } } /* PART I: */ /* ======= */ /* Deflate the Floquet multiplier at +1.0 so that the deflated */ /* circuit pencil is not defective at periodic branch turning points. */ /* The matrix (C0 - C1) should be (nearly) singular. Find an approximation*/ /* to the right null vector (call it X). This will be our approximation */ /* to the eigenvector corresponding to the fixed multiplier at +1.0. */ /* There are many ways to get this approximation. We could use */ /* 1) p'(0) = f(p(0)) */ /* 2) AUTO'86 routine NLVC applied to C0-C1 */ /* 3) the right singular vector corresponding to the smallest */ /* singular value of C0-C1 */ /* I 've chosen option 3) because it should introduce as little roundof f */ /* error as possible. Although it is more expensive, this is insignificant*/ /* relative to the rest of the AUTO computations. Also, the SVD does give a*/ /* version of the Householder matrix which we would have to compute */ /* anyways. But note that it gives V = ( X perp | X ) and not (X | Xperp) ,*/ /* which the Householder routine would give. This will permute the deflated*/ /* circuit pencil, so that the part to be deflated is in the last column, */ /* not it the first column, as was shown in the paper. */ i_1 = *ndim; for (j = 1; j <= i_1; ++j) { i_2 = *ndim; for (i = 1; i <= i_2; ++i) { rwork[i + j * rwork_dim1] = c0[i + j * c0_dim1] - c1[i + j * c1_dim1]; /* L9: */ } /* L10: */ } ezsvd_(&rwork[rwork_offset], ndim, ndim, ndim, svds, svde, svdu, &c__1, svdv, &c__50, svdwrk, &c__1, &svdinf, &c_b24); if (svdinf != 0) { s_wsfe(&io__15); do_fio(&c__1, (char *)&svdinf, (ftnlen)sizeof(integer)); e_wsfe(); } /* Apply a Householder matrix (call it H1) based on the null vector */ /* to (C0, C1) from the right. H1 = SVDV = ( Xperp | X ), where X */ /* is the null vector. */ dgemm_("n", "n", ndim, ndim, ndim, &c_b30, &c0[c0_offset], ndim, svdv, & c__50, &c_b32, &rwork[rwork_offset], ndim, 1L, 1L); dgemc_(ndim, ndim, &rwork[rwork_offset], ndim, &c0[c0_offset], ndim, & c__0); dgemm_("n", "n", ndim, ndim, ndim, &c_b30, &c1[c1_offset], ndim, svdv, & c__50, &c_b32, &rwork[rwork_offset], ndim, 1L, 1L); dgemc_(ndim, ndim, &rwork[rwork_offset], ndim, &c1[c1_offset], ndim, & c__0); /* Apply a Householder matrix (call it H2) based on */ /* (C0*X/||C0*X|| + C1*X/||C1*X||) / 2 */ /* to (C0*H1, C1*H1) from the left. */ nrmc0x = dnrm2_(ndim, &c0[*ndim * c0_dim1 + 1], &c__1); nrmc1x = dnrm2_(ndim, &c1[*ndim * c1_dim1 + 1], &c__1); i_1 = *ndim; for (i = 1; i <= i_1; ++i) { x[i - 1] = (c0[i + *ndim * c0_dim1] / nrmc0x + c1[i + *ndim * c1_dim1] / nrmc1x) / 2.; /* L20: */ } dhhpr_(&c__1, ndim, ndim, x, &c__1, &beta, v); dhhap_(&c__1, ndim, ndim, ndim, &beta, v, &c__1, &c0[c0_offset], ndim); dhhap_(&c__1, ndim, ndim, ndim, &beta, v, &c__1, &c1[c1_offset], ndim); /* Rescale so that (H2^T)*C0*(H1)(1,NDIM) ~= (H2^T)*C1*(H1)(1,NDIM) ~= 1.0*/ /* Computing MAX */ d_3 = (d_1 = c0[*ndim * c0_dim1 + 1], abs(d_1)), d_4 = (d_2 = c1[*ndim * c1_dim1 + 1], abs(d_2)); const_ = max(d_4,d_3); i_1 = *ndim; for (j = 1; j <= i_1; ++j) { i_2 = *ndim; for (i = 1; i <= i_2; ++i) { c0[i + j * c0_dim1] /= const_; c1[i + j * c1_dim1] /= const_; /* L29: */ } /* L30: */ } /* Finished the deflation process! Print the deflated circuit pencil. */ if (*iid > 2) { s_wsfe(&io__22); e_wsfe(); s_wsfe(&io__23); e_wsfe(); i_1 = *ndim; for (i = 1; i <= i_1; ++i) { s_wsfe(&io__24); i_2 = *ndim; for (j = 1; j <= i_2; ++j) { do_fio(&c__1, (char *)&c0[i + j * c0_dim1], (ftnlen)sizeof( doublereal)); } e_wsfe(); /* L40: */ } s_wsfe(&io__25); e_wsfe(); i_1 = *ndim; for (i = 1; i <= i_1; ++i) { s_wsfe(&io__26); i_2 = *ndim; for (j = 1; j <= i_2; ++j) { do_fio(&c__1, (char *)&c1[i + j * c1_dim1], (ftnlen)sizeof( doublereal)); } e_wsfe(); /* L41: */ } } /* At this point we have */ /* (C0Bar, C1Bar) */ /* ::= (H2^T)*(C0, C1)*(H1). */ /* (( B0^T | Beta0 ) ( B1^T | Beta1 )) 1 */ /* = (( ----------------- ), ( ----------------- )) */ /* (( C0BarDef | Delta0 ) ( C1BarDef | Delta1 )) NDIM-1 */ /* NDIM-1 1 NDIM-1 1 */ /* and approximations to the Floquet multipliers are */ /* (Beta0/Beta1) union the eigenvalues of the deflated pencil */ /* (C0BarDef, C1BarDef). */ /* PART II: */ /* ======== */ /* Compute the eigenvalues of the deflated circuit pencil */ /* (C0BarDef, C1BarDef) */ /* by using the QZ routines from EISPAC. */ ndimm1 = *ndim - 1; /* reduce the generalized eigenvalue problem to a simpler form */ /* (C0BarDef,C1BarDef) = (upper hessenberg, upper triangular) */ qzhes_(ndim, &ndimm1, &c0[c0_dim1 + 2], &c1[c1_dim1 + 2], &c__0, qzz); /* now reduce to an even simpler form */ /* (C0BarDef,C1BarDef) = (quasi-upper triangular, upper triangular) */ qzit_(ndim, &ndimm1, &c0[c0_dim1 + 2], &c1[c1_dim1 + 2], &c_b32, &c__0, qzz, &qzierr); if (qzierr != 0) { s_wsfe(&io__30); do_fio(&c__1, (char *)&qzierr, (ftnlen)sizeof(integer)); e_wsfe(); } /* compute the generalized eigenvalues */ qzval_(ndim, &ndimm1, &c0[c0_dim1 + 2], &c1[c1_dim1 + 2], qzalfr, qzalfi, qzbeta, &c__0, qzz); /* Pack the eigenvalues into complex form. */ d_1 = c0[*ndim * c0_dim1 + 1] / c1[*ndim * c1_dim1 + 1]; z_1.r = d_1, z_1.i = 0.; ev[1].r = z_1.r, ev[1].i = z_1.i; infev = FALSE_; i_1 = ndimm1; for (j = 1; j <= i_1; ++j) { if (qzbeta[j - 1] != 0.) { i_2 = j + 1; d_1 = qzalfr[j - 1] / qzbeta[j - 1]; d_2 = qzalfi[j - 1] / qzbeta[j - 1]; z_1.r = d_1, z_1.i = d_2; ev[i_2].r = z_1.r, ev[i_2].i = z_1.i; } else { i_2 = j + 1; ev[i_2].r = 1e99, ev[i_2].i = 1e99; infev = TRUE_; } /* L50: */ } if (infev) { s_wsfe(&io__35); e_wsfe(); } /* Done! */ return 0; /* Format statements */ } /* flowkm_ */ /* ************************** */ /* * Householder routines * */ /* ************************** */ /* Subroutines for performing Householder plane rotations. */ /* DHHPR: for computing Householder transformations and */ /* DHHAP: for applying them. */ /* Ref: Golub and van Loan, Matrix Calcualtions, */ /* First Edition, Pages 38-43 */ /* Subroutine */ int dhhpr_(k, j, n, x, incx, beta, v) integer *k, *j, *n; doublereal *x; integer *incx; doublereal *beta, *v; { /* System generated locals */ integer i_1, i_2; doublereal d_1; /* Builtin functions */ integer s_wsle(), do_lio(), e_wsle(); /* Subroutine */ int s_stop(); double d_sign(); /* Local variables */ static integer iend, jmkp1; extern doublereal dnrm2_(); static integer i, l; static doublereal m, alpha; extern integer idamax_(); static integer istart; /* Fortran I/O blocks */ static cilist io__36 = { 0, 6, 0, 0, 0 }; static cilist io__37 = { 0, 6, 0, 0, 0 }; static cilist io__38 = { 0, 6, 0, 0, 0 }; /* Parameter adjustments */ --x; --v; /* Function Body */ /* IMPLICIT UNDEFINED (A-Z,a-z) */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* DHHPR computes a Householder Plane Rotation (G&vL Alg. 3.3-1) */ /* defined by v and beta. */ /* (I - beta v vt) * x is such that x_i = 0 for i=k+1 to j. */ /* Parameters */ /* ========== */ /* K - INTEGER. */ /* On entry, K specifies that the K+1st entry of X */ /* be the first to be zeroed. */ /* K must be at least one. */ /* Unchanged on exit. */ /* J - INTEGER. */ /* On entry, J specifies the last entry of X to be zeroed. */ /* J must be >= K and <= N. */ /* Unchanged on exit. */ /* N - INTEGER. */ /* On entry, N specifies the (logical) length of X. */ /* Unchanged on exit. */ /* X - DOUBLE PRECISION array of DIMENSION at least */ /* ( 1 + ( N - 1 )*abs( INCX ) ). */ /* On entry, X specifies the vector to be (partially) zeroed. */ /* Unchanged on exit. */ /* INCX - INTEGER. */ /* On entry, INCX specifies the increment for the elements of */ /* X. INCX must be > zero. If X represents part of a matrix, */ /* then use INCX = 1 if a column vector is being zeroed and */ /* INCX = NDIM if a row vector is being zeroed. */ /* Unchanged on exit. */ /* BETA - DOUBLE PRECISION. */ /* BETA specifies the scalar beta. (see pg. 40 of G and v.L.) */ /* V - DOUBLE PRECISION array of DIMENSION at least n. */ /* Is updated to be the appropriate Householder vector for */ /* the given problem. (Note: space for the implicit zeroes is */ /* assumed to be present. Will save on time for index translation.)*/ /* -- Written by Tom Fairgrieve, */ /* Department of Computer Science, */ /* University of Toronto, */ /* Toronto, Ontario CANADA M5S 1A4 */ /* .. Local Scalars .. */ /* .. External Functions from BLAS .. */ /* .. External Subroutines from BLAS .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Test the input parameters. */ if (*k < 1 || *k > *j) { s_wsle(&io__36); do_lio(&c__9, &c__1, "Domain error for K in DHHPR", 27L); e_wsle(); s_stop("", 0L); } if (*j > *n) { s_wsle(&io__37); do_lio(&c__9, &c__1, "Domain error for J in DHHPR", 27L); e_wsle(); s_stop("", 0L); } if (*incx < 1) { s_wsle(&io__38); do_lio(&c__9, &c__1, "Domain error for INCX in DHHPR", 30L); e_wsle(); s_stop("", 0L); } /* Number of potential non-zero elements in V. */ jmkp1 = *j - *k + 1; /* Find M := max{ |x_k|, ... , |x_j| } */ m = (d_1 = x[idamax_(&jmkp1, &x[*k], incx)], abs(d_1)); /* alpha := 0 */ /* For i = k to j */ /* v_i = x_i / m */ /* alpha := alpha + v_i^2 (i.e. alpha = vtv) */ /* End For */ /* alpha := sqrt( alpha ) */ /* Copy X(K)/M, ... , X(J)/M to V(K), ... , V(J) */ if (*incx == 1) { i_1 = *j; for (i = *k; i <= i_1; ++i) { v[i] = x[i] / m; /* L10: */ } } else { iend = jmkp1 * *incx; istart = (*k - 1) * *incx + 1; l = *k; i_1 = iend; i_2 = *incx; for (i = istart; i_2 < 0 ? i >= i_1 : i <= i_1; i += i_2) { v[l] = x[i] / m; ++l; /* L20: */ } } /* Compute alpha */ alpha = dnrm2_(&jmkp1, &v[*k], &c__1); /* beta := 1/(alpha(alpha + |V_k|)) */ *beta = 1. / (alpha * (alpha + (d_1 = v[*k], abs(d_1)))); /* v_k := v_k + sign(v_k)*alpha */ v[*k] += d_sign(&c_b30, &v[*k]) * alpha; /* Done ! */ return 0; /* End of DHHPR. */ } /* dhhpr_ */ /* Subroutine */ int dhhap_(k, j, n, q, beta, v, job, a, lda) integer *k, *j, *n, *q; doublereal *beta, *v; integer *job; doublereal *a; integer *lda; { /* System generated locals */ integer a_dim1, a_offset, i_1, i_2; /* Builtin functions */ integer s_wsle(), do_lio(), e_wsle(); /* Subroutine */ int s_stop(); /* Local variables */ extern doublereal ddot_(); static integer jmkp1; static doublereal s; static integer col, row; /* Fortran I/O blocks */ static cilist io__46 = { 0, 6, 0, 0, 0 }; static cilist io__47 = { 0, 6, 0, 0, 0 }; static cilist io__48 = { 0, 6, 0, 0, 0 }; static cilist io__49 = { 0, 6, 0, 0, 0 }; /* Parameter adjustments */ --v; a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; /* Function Body */ /* .. Scalar Arguments .. */ /* .. Array Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* DHHAP applies a Householder Plane Rotation defined by v and beta */ /* to the matrix A. If JOB = 1 then A := (I - beta*v*vt)A and if */ /* JOB = 2 then A := A(I - beta*v*vt). (See Golub and van Loan */ /* Alg. 3.3-2.) */ /* Parameters */ /* ========== */ /* K - INTEGER. */ /* On entry, K specifies that the V(K) may be the first */ /* non-zero entry of V. */ /* K must be at least one. */ /* Unchanged on exit. */ /* J - INTEGER. */ /* On entry, J specifies the last non-zero entry of V. */ /* J must be >= K and <= N. */ /* Unchanged on exit. */ /* N - INTEGER. */ /* On entry, N specifies the row dimension of A. */ /* Unchanged on exit. */ /* Q - INTEGER. */ /* On entry, Q specifies the column dimension of A. */ /* Unchanged on exit. */ /* BETA - DOUBLE PRECISION. */ /* BETA specifies the scalar beta. (see pg. 40 of G and v.L.) */ /* Unchanged on exit. */ /* V - DOUBLE PRECISION array of DIMENSION at least n. */ /* Householder vector v. */ /* Unchanged on exit. */ /* JOB - INTEGER. */ /* On entry, JOB specifies the order of the Householder application.*/ /* If JOB = 1 then A := (I - beta*v*vt)A and if JOB = 2 then */ /* A := A(I - beta*v*vt) */ /* Unchanged on exit. */ /* A - DOUBLE PRECISION array of DIMENSION at least */ /* ( LDA, Q ). */ /* On entry, A specifies the matrix to be transformed. */ /* On exit, A specifies the transformed matrix. */ /* LDA - INTEGER. */ /* On entry, LDA specifies the declared leading dimension of A. */ /* Unchanged on exit. */ /* -- Written by Tom Fairgrieve, */ /* Department of Computer Science, */ /* University of Toronto, */ /* Toronto, Ontario CANADA M5S 1A4 */ /* .. Local Scalars .. */ /* .. External Functions from BLAS .. */ /* .. Executable Statements .. */ /* Test the input parameters. */ if (*job != 1 && *job != 2) { s_wsle(&io__46); do_lio(&c__9, &c__1, "Domain error for JOB in DHHAP", 29L); e_wsle(); s_stop("", 0L); } if (*k < 1 || *k > *j) { s_wsle(&io__47); do_lio(&c__9, &c__1, "Domain error for K in DHHAP", 27L); e_wsle(); s_stop("", 0L); } if (*job == 1) { if (*j > *n) { s_wsle(&io__48); do_lio(&c__9, &c__1, "Domain error for J in DHHAP", 27L); e_wsle(); s_stop("", 0L); } } else { if (*j > *q) { s_wsle(&io__49); do_lio(&c__9, &c__1, "Domain error for J in DHHAP", 27L); e_wsle(); s_stop("", 0L); } } /* Minimum {row,col} dimension of update. */ jmkp1 = *j - *k + 1; /* If (JOB = 1) then */ /* For p = 1, ... , q */ /* s := beta*(v_k*a_k,p + ... + v_j*a_j,p) */ /* For i = k, ..., j */ /* a_i,p := a_i,p - s*v_i */ /* End For */ /* End For */ /* Else % JOB=2 */ /* For p = 1, ... , n */ /* s := beta*(v_k*a_p,k + ... + v_j*a_p,j) */ /* For i = k, ..., j */ /* a_p,i := a_p,i - s*v_i */ /* End For */ /* End For */ /* End If */ if (*job == 1) { i_1 = *q; for (col = 1; col <= i_1; ++col) { s = *beta * ddot_(&jmkp1, &v[*k], &c__1, &a[*k + col * a_dim1], & c__1); i_2 = *j; for (row = *k; row <= i_2; ++row) { a[row + col * a_dim1] -= s * v[row]; /* L9: */ } /* L10: */ } } else { i_1 = *n; for (row = 1; row <= i_1; ++row) { s = *beta * ddot_(&jmkp1, &v[*k], &c__1, &a[row + *k * a_dim1], lda); i_2 = *j; for (col = *k; col <= i_2; ++col) { a[row + col * a_dim1] -= s * v[col]; /* L19: */ } /* L20: */ } } /* Done ! */ return 0; /* End of DHHAP. */ } /* dhhap_ */ /* ************************** */ /* * Routines from EISPAC * */ /* ************************** */ /* Subroutine */ int qzhes_(nm, n, a, b, matz, z) integer *nm, *n; doublereal *a, *b; logical *matz; doublereal *z; { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i_1, i_2, i_3; doublereal d_1, d_2; /* Builtin functions */ double sqrt(), d_sign(); /* Local variables */ static integer i, j, k, l; static doublereal r, s, t; static integer l1; static doublereal u1, u2, v1, v2; static integer lb, nk1, nm1, nm2; static doublereal rho; /* Parameter adjustments */ a_dim1 = *nm; a_offset = a_dim1 + 1; a -= a_offset; b_dim1 = *nm; b_offset = b_dim1 + 1; b -= b_offset; z_dim1 = *nm; z_offset = z_dim1 + 1; z -= z_offset; /* Function Body */ /* THIS SUBROUTINE IS THE FIRST STEP OF THE QZ ALGORITHM */ /* FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */ /* SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. */ /* THIS SUBROUTINE ACCEPTS A PAIR OF REAL GENERAL MATRICES AND */ /* REDUCES ONE OF THEM TO UPPER HESSENBERG FORM AND THE OTHER */ /* TO UPPER TRIANGULAR FORM USING ORTHOGONAL TRANSFORMATIONS. */ /* IT IS USUALLY FOLLOWED BY QZIT, QZVAL AND, POSSIBLY, QZVEC. */ /* ON INPUT */ /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ /* DIMENSION STATEMENT. */ /* N IS THE ORDER OF THE MATRICES. */ /* A CONTAINS A REAL GENERAL MATRIX. */ /* B CONTAINS A REAL GENERAL MATRIX. */ /* MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS */ /* ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */ /* EIGENVECTORS, AND TO .FALSE. OTHERWISE. */ /* ON OUTPUT */ /* A HAS BEEN REDUCED TO UPPER HESSENBERG FORM. THE ELEMENTS */ /* BELOW THE FIRST SUBDIAGONAL HAVE BEEN SET TO ZERO. */ /* B HAS BEEN REDUCED TO UPPER TRIANGULAR FORM. THE ELEMENTS */ /* BELOW THE MAIN DIAGONAL HAVE BEEN SET TO ZERO. */ /* Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS IF */ /* MATZ HAS BEEN SET TO .TRUE. OTHERWISE, Z IS NOT REFERENCED. */ /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY */ /* THIS VERSION DATED AUGUST 1983. */ /* ------------------------------------------------------------------ */ /* .......... INITIALIZE Z .......... */ if (! (*matz)) { goto L10; } i_1 = *n; for (j = 1; j <= i_1; ++j) { i_2 = *n; for (i = 1; i <= i_2; ++i) { z[i + j * z_dim1] = 0.; /* L2: */ } z[j + j * z_dim1] = 1.; /* L3: */ } /* .......... REDUCE B TO UPPER TRIANGULAR FORM .......... */ L10: if (*n <= 1) { goto L170; } nm1 = *n - 1; i_1 = nm1; for (l = 1; l <= i_1; ++l) { l1 = l + 1; s = 0.; i_2 = *n; for (i = l1; i <= i_2; ++i) { s += (d_1 = b[i + l * b_dim1], abs(d_1)); /* L20: */ } if (s == 0.) { goto L100; } s += (d_1 = b[l + l * b_dim1], abs(d_1)); r = 0.; i_2 = *n; for (i = l; i <= i_2; ++i) { b[i + l * b_dim1] /= s; /* Computing 2nd power */ d_1 = b[i + l * b_dim1]; r += d_1 * d_1; /* L25: */ } d_1 = sqrt(r); r = d_sign(&d_1, &b[l + l * b_dim1]); b[l + l * b_dim1] += r; rho = r * b[l + l * b_dim1]; i_2 = *n; for (j = l1; j <= i_2; ++j) { t = 0.; i_3 = *n; for (i = l; i <= i_3; ++i) { t += b[i + l * b_dim1] * b[i + j * b_dim1]; /* L30: */ } t = -t / rho; i_3 = *n; for (i = l; i <= i_3; ++i) { b[i + j * b_dim1] += t * b[i + l * b_dim1]; /* L40: */ } /* L50: */ } i_2 = *n; for (j = 1; j <= i_2; ++j) { t = 0.; i_3 = *n; for (i = l; i <= i_3; ++i) { t += b[i + l * b_dim1] * a[i + j * a_dim1]; /* L60: */ } t = -t / rho; i_3 = *n; for (i = l; i <= i_3; ++i) { a[i + j * a_dim1] += t * b[i + l * b_dim1]; /* L70: */ } /* L80: */ } b[l + l * b_dim1] = -s * r; i_2 = *n; for (i = l1; i <= i_2; ++i) { b[i + l * b_dim1] = 0.; /* L90: */ } L100: ;} /* .......... REDUCE A TO UPPER HESSENBERG FORM, WHILE */ /* KEEPING B TRIANGULAR .......... */ if (*n == 2) { goto L170; } nm2 = *n - 2; i_1 = nm2; for (k = 1; k <= i_1; ++k) { nk1 = nm1 - k; /* .......... FOR L=N-1 STEP -1 UNTIL K+1 DO -- .......... */ i_2 = nk1; for (lb = 1; lb <= i_2; ++lb) { l = *n - lb; l1 = l + 1; /* .......... ZERO A(L+1,K) .......... */ s = (d_1 = a[l + k * a_dim1], abs(d_1)) + (d_2 = a[l1 + k * a_dim1], abs(d_2)); if (s == 0.) { goto L150; } u1 = a[l + k * a_dim1] / s; u2 = a[l1 + k * a_dim1] / s; d_1 = sqrt(u1 * u1 + u2 * u2); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; u2 = v2 / v1; i_3 = *n; for (j = k; j <= i_3; ++j) { t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1]; a[l + j * a_dim1] += t * v1; a[l1 + j * a_dim1] += t * v2; /* L110: */ } a[l1 + k * a_dim1] = 0.; i_3 = *n; for (j = l; j <= i_3; ++j) { t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1]; b[l + j * b_dim1] += t * v1; b[l1 + j * b_dim1] += t * v2; /* L120: */ } /* .......... ZERO B(L+1,L) .......... */ s = (d_1 = b[l1 + l1 * b_dim1], abs(d_1)) + (d_2 = b[l1 + l * b_dim1], abs(d_2)); if (s == 0.) { goto L150; } u1 = b[l1 + l1 * b_dim1] / s; u2 = b[l1 + l * b_dim1] / s; d_1 = sqrt(u1 * u1 + u2 * u2); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; u2 = v2 / v1; i_3 = l1; for (i = 1; i <= i_3; ++i) { t = b[i + l1 * b_dim1] + u2 * b[i + l * b_dim1]; b[i + l1 * b_dim1] += t * v1; b[i + l * b_dim1] += t * v2; /* L130: */ } b[l1 + l * b_dim1] = 0.; i_3 = *n; for (i = 1; i <= i_3; ++i) { t = a[i + l1 * a_dim1] + u2 * a[i + l * a_dim1]; a[i + l1 * a_dim1] += t * v1; a[i + l * a_dim1] += t * v2; /* L140: */ } if (! (*matz)) { goto L150; } i_3 = *n; for (i = 1; i <= i_3; ++i) { t = z[i + l1 * z_dim1] + u2 * z[i + l * z_dim1]; z[i + l1 * z_dim1] += t * v1; z[i + l * z_dim1] += t * v2; /* L145: */ } L150: ;} /* L160: */ } L170: return 0; } /* qzhes_ */ /* Subroutine */ int qzit_(nm, n, a, b, eps1, matz, z, ierr) integer *nm, *n; doublereal *a, *b, *eps1; logical *matz; doublereal *z; integer *ierr; { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i_1, i_2, i_3; doublereal d_1, d_2, d_3; /* Builtin functions */ double sqrt(), d_sign(); /* Local variables */ static doublereal epsa, epsb; static integer i, j, k, l; static doublereal r, s, t, anorm, bnorm; static integer enorn; static doublereal a1, a2, a3; static integer k1, k2, l1; static doublereal u1, u2, u3, v1, v2, v3, a11, a12, a21, a22, a33, a34, a43, a44, b11, b12, b22, b33; static integer na, ld; static doublereal b34, b44; static integer en; static doublereal ep; static integer ll; static doublereal sh; extern doublereal epslon_(); static logical notlas; static integer km1, lm1; static doublereal ani, bni; static integer ish, itn, its, enm2, lor1; /* Parameter adjustments */ a_dim1 = *nm; a_offset = a_dim1 + 1; a -= a_offset; b_dim1 = *nm; b_offset = b_dim1 + 1; b -= b_offset; z_dim1 = *nm; z_offset = z_dim1 + 1; z -= z_offset; /* Function Body */ /* THIS SUBROUTINE IS THE SECOND STEP OF THE QZ ALGORITHM */ /* FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */ /* SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART, */ /* AS MODIFIED IN TECHNICAL NOTE NASA TN D-7305(1973) BY WARD. */ /* THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM */ /* IN UPPER HESSENBERG FORM AND THE OTHER IN UPPER TRIANGULAR FORM. */ /* IT REDUCES THE HESSENBERG MATRIX TO QUASI-TRIANGULAR FORM USING */ /* ORTHOGONAL TRANSFORMATIONS WHILE MAINTAINING THE TRIANGULAR FORM */ /* OF THE OTHER MATRIX. IT IS USUALLY PRECEDED BY QZHES AND */ /* FOLLOWED BY QZVAL AND, POSSIBLY, QZVEC. */ /* ON INPUT */ /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ /* DIMENSION STATEMENT. */ /* N IS THE ORDER OF THE MATRICES. */ /* A CONTAINS A REAL UPPER HESSENBERG MATRIX. */ /* B CONTAINS A REAL UPPER TRIANGULAR MATRIX. */ /* EPS1 IS A TOLERANCE USED TO DETERMINE NEGLIGIBLE ELEMENTS. */ /* EPS1 = 0.0 (OR NEGATIVE) MAY BE INPUT, IN WHICH CASE AN */ /* ELEMENT WILL BE NEGLECTED ONLY IF IT IS LESS THAN ROUNDOFF */ /* ERROR TIMES THE NORM OF ITS MATRIX. IF THE INPUT EPS1 IS */ /* POSITIVE, THEN AN ELEMENT WILL BE CONSIDERED NEGLIGIBLE */ /* IF IT IS LESS THAN EPS1 TIMES THE NORM OF ITS MATRIX. A */ /* POSITIVE VALUE OF EPS1 MAY RESULT IN FASTER EXECUTION, */ /* BUT LESS ACCURATE RESULTS. */ /* MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS */ /* ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */ /* EIGENVECTORS, AND TO .FALSE. OTHERWISE. */ /* Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE */ /* TRANSFORMATION MATRIX PRODUCED IN THE REDUCTION */ /* BY QZHES, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. */ /* IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. */ /* ON OUTPUT */ /* A HAS BEEN REDUCED TO QUASI-TRIANGULAR FORM. THE ELEMENTS */ /* BELOW THE FIRST SUBDIAGONAL ARE STILL ZERO AND NO TWO */ /* CONSECUTIVE SUBDIAGONAL ELEMENTS ARE NONZERO. */ /* B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS */ /* HAVE BEEN ALTERED. THE LOCATION B(N,1) IS USED TO STORE */ /* EPS1 TIMES THE NORM OF B FOR LATER USE BY QZVAL AND QZVEC. */ /* Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS */ /* (FOR BOTH STEPS) IF MATZ HAS BEEN SET TO .TRUE.. */ /* IERR IS SET TO */ /* ZERO FOR NORMAL RETURN, */ /* J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED */ /* WHILE THE J-TH EIGENVALUE IS BEING SOUGHT. */ /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY */ /* THIS VERSION DATED AUGUST 1983. */ /* ------------------------------------------------------------------ */ *ierr = 0; /* .......... COMPUTE EPSA,EPSB .......... */ anorm = 0.; bnorm = 0.; i_1 = *n; for (i = 1; i <= i_1; ++i) { ani = 0.; if (i != 1) { ani = (d_1 = a[i + (i - 1) * a_dim1], abs(d_1)); } bni = 0.; i_2 = *n; for (j = i; j <= i_2; ++j) { ani += (d_1 = a[i + j * a_dim1], abs(d_1)); bni += (d_1 = b[i + j * b_dim1], abs(d_1)); /* L20: */ } if (ani > anorm) { anorm = ani; } if (bni > bnorm) { bnorm = bni; } /* L30: */ } if (anorm == 0.) { anorm = 1.; } if (bnorm == 0.) { bnorm = 1.; } ep = *eps1; if (ep > 0.) { goto L50; } /* .......... USE ROUNDOFF LEVEL IF EPS1 IS ZERO .......... */ ep = epslon_(&c_b30); L50: epsa = ep * anorm; epsb = ep * bnorm; /* .......... REDUCE A TO QUASI-TRIANGULAR FORM, WHILE */ /* KEEPING B TRIANGULAR .......... */ lor1 = 1; enorn = *n; en = *n; itn = *n * 30; /* .......... BEGIN QZ STEP .......... */ L60: if (en <= 2) { goto L1001; } if (! (*matz)) { enorn = en; } its = 0; na = en - 1; enm2 = na - 1; L70: ish = 2; /* .......... CHECK FOR CONVERGENCE OR REDUCIBILITY. */ /* FOR L=EN STEP -1 UNTIL 1 DO -- .......... */ i_1 = en; for (ll = 1; ll <= i_1; ++ll) { lm1 = en - ll; l = lm1 + 1; if (l == 1) { goto L95; } if ((d_1 = a[l + lm1 * a_dim1], abs(d_1)) <= epsa) { goto L90; } /* L80: */ } L90: a[l + lm1 * a_dim1] = 0.; if (l < na) { goto L95; } /* .......... 1-BY-1 OR 2-BY-2 BLOCK ISOLATED .......... */ en = lm1; goto L60; /* .......... CHECK FOR SMALL TOP OF B .......... */ L95: ld = l; L100: l1 = l + 1; b11 = b[l + l * b_dim1]; if (abs(b11) > epsb) { goto L120; } b[l + l * b_dim1] = 0.; s = (d_1 = a[l + l * a_dim1], abs(d_1)) + (d_2 = a[l1 + l * a_dim1], abs( d_2)); u1 = a[l + l * a_dim1] / s; u2 = a[l1 + l * a_dim1] / s; d_1 = sqrt(u1 * u1 + u2 * u2); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; u2 = v2 / v1; i_1 = enorn; for (j = l; j <= i_1; ++j) { t = a[l + j * a_dim1] + u2 * a[l1 + j * a_dim1]; a[l + j * a_dim1] += t * v1; a[l1 + j * a_dim1] += t * v2; t = b[l + j * b_dim1] + u2 * b[l1 + j * b_dim1]; b[l + j * b_dim1] += t * v1; b[l1 + j * b_dim1] += t * v2; /* L110: */ } if (l != 1) { a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1]; } lm1 = l; l = l1; goto L90; L120: a11 = a[l + l * a_dim1] / b11; a21 = a[l1 + l * a_dim1] / b11; if (ish == 1) { goto L140; } /* .......... ITERATION STRATEGY .......... */ if (itn == 0) { goto L1000; } if (its == 10) { goto L155; } /* .......... DETERMINE TYPE OF SHIFT .......... */ b22 = b[l1 + l1 * b_dim1]; if (abs(b22) < epsb) { b22 = epsb; } b33 = b[na + na * b_dim1]; if (abs(b33) < epsb) { b33 = epsb; } b44 = b[en + en * b_dim1]; if (abs(b44) < epsb) { b44 = epsb; } a33 = a[na + na * a_dim1] / b33; a34 = a[na + en * a_dim1] / b44; a43 = a[en + na * a_dim1] / b33; a44 = a[en + en * a_dim1] / b44; b34 = b[na + en * b_dim1] / b44; t = (a43 * b34 - a33 - a44) * .5; r = t * t + a34 * a43 - a33 * a44; if (r < 0.) { goto L150; } /* .......... DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A .......... */ ish = 1; r = sqrt(r); sh = -t + r; s = -t - r; if ((d_1 = s - a44, abs(d_1)) < (d_2 = sh - a44, abs(d_2))) { sh = s; } /* .......... LOOK FOR TWO CONSECUTIVE SMALL */ /* SUB-DIAGONAL ELEMENTS OF A. */ /* FOR L=EN-2 STEP -1 UNTIL LD DO -- .......... */ i_1 = enm2; for (ll = ld; ll <= i_1; ++ll) { l = enm2 + ld - ll; if (l == ld) { goto L140; } lm1 = l - 1; l1 = l + 1; t = a[l + l * a_dim1]; if ((d_1 = b[l + l * b_dim1], abs(d_1)) > epsb) { t -= sh * b[l + l * b_dim1]; } if ((d_1 = a[l + lm1 * a_dim1], abs(d_1)) <= (d_2 = t / a[l1 + l * a_dim1], abs(d_2)) * epsa) { goto L100; } /* L130: */ } L140: a1 = a11 - sh; a2 = a21; if (l != ld) { a[l + lm1 * a_dim1] = -a[l + lm1 * a_dim1]; } goto L160; /* .......... DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A .......... */ L150: a12 = a[l + l1 * a_dim1] / b22; a22 = a[l1 + l1 * a_dim1] / b22; b12 = b[l + l1 * b_dim1] / b22; a1 = ((a33 - a11) * (a44 - a11) - a34 * a43 + a43 * b34 * a11) / a21 + a12 - a11 * b12; a2 = a22 - a11 - a21 * b12 - (a33 - a11) - (a44 - a11) + a43 * b34; a3 = a[l1 + 1 + l1 * a_dim1] / b22; goto L160; /* .......... AD HOC SHIFT .......... */ L155: a1 = 0.; a2 = 1.; a3 = 1.1605; L160: ++its; --itn; if (! (*matz)) { lor1 = ld; } /* .......... MAIN LOOP .......... */ i_1 = na; for (k = l; k <= i_1; ++k) { notlas = k != na && ish == 2; k1 = k + 1; k2 = k + 2; /* Computing MAX */ i_2 = k - 1; km1 = max(l,i_2); /* Computing MAX */ i_2 = en, i_3 = k1 + ish; ll = min(i_3,i_2); if (notlas) { goto L190; } /* .......... ZERO A(K+1,K-1) .......... */ if (k == l) { goto L170; } a1 = a[k + km1 * a_dim1]; a2 = a[k1 + km1 * a_dim1]; L170: s = abs(a1) + abs(a2); if (s == 0.) { goto L70; } u1 = a1 / s; u2 = a2 / s; d_1 = sqrt(u1 * u1 + u2 * u2); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; u2 = v2 / v1; i_2 = enorn; for (j = km1; j <= i_2; ++j) { t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1]; a[k + j * a_dim1] += t * v1; a[k1 + j * a_dim1] += t * v2; t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1]; b[k + j * b_dim1] += t * v1; b[k1 + j * b_dim1] += t * v2; /* L180: */ } if (k != l) { a[k1 + km1 * a_dim1] = 0.; } goto L240; /* .......... ZERO A(K+1,K-1) AND A(K+2,K-1) .......... */ L190: if (k == l) { goto L200; } a1 = a[k + km1 * a_dim1]; a2 = a[k1 + km1 * a_dim1]; a3 = a[k2 + km1 * a_dim1]; L200: s = abs(a1) + abs(a2) + abs(a3); if (s == 0.) { goto L260; } u1 = a1 / s; u2 = a2 / s; u3 = a3 / s; d_1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; v3 = -u3 / r; u2 = v2 / v1; u3 = v3 / v1; i_2 = enorn; for (j = km1; j <= i_2; ++j) { t = a[k + j * a_dim1] + u2 * a[k1 + j * a_dim1] + u3 * a[k2 + j * a_dim1]; a[k + j * a_dim1] += t * v1; a[k1 + j * a_dim1] += t * v2; a[k2 + j * a_dim1] += t * v3; t = b[k + j * b_dim1] + u2 * b[k1 + j * b_dim1] + u3 * b[k2 + j * b_dim1]; b[k + j * b_dim1] += t * v1; b[k1 + j * b_dim1] += t * v2; b[k2 + j * b_dim1] += t * v3; /* L210: */ } if (k == l) { goto L220; } a[k1 + km1 * a_dim1] = 0.; a[k2 + km1 * a_dim1] = 0.; /* .......... ZERO B(K+2,K+1) AND B(K+2,K) .......... */ L220: s = (d_1 = b[k2 + k2 * b_dim1], abs(d_1)) + (d_2 = b[k2 + k1 * b_dim1] , abs(d_2)) + (d_3 = b[k2 + k * b_dim1], abs(d_3)); if (s == 0.) { goto L240; } u1 = b[k2 + k2 * b_dim1] / s; u2 = b[k2 + k1 * b_dim1] / s; u3 = b[k2 + k * b_dim1] / s; d_1 = sqrt(u1 * u1 + u2 * u2 + u3 * u3); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; v3 = -u3 / r; u2 = v2 / v1; u3 = v3 / v1; i_2 = ll; for (i = lor1; i <= i_2; ++i) { t = a[i + k2 * a_dim1] + u2 * a[i + k1 * a_dim1] + u3 * a[i + k * a_dim1]; a[i + k2 * a_dim1] += t * v1; a[i + k1 * a_dim1] += t * v2; a[i + k * a_dim1] += t * v3; t = b[i + k2 * b_dim1] + u2 * b[i + k1 * b_dim1] + u3 * b[i + k * b_dim1]; b[i + k2 * b_dim1] += t * v1; b[i + k1 * b_dim1] += t * v2; b[i + k * b_dim1] += t * v3; /* L230: */ } b[k2 + k * b_dim1] = 0.; b[k2 + k1 * b_dim1] = 0.; if (! (*matz)) { goto L240; } i_2 = *n; for (i = 1; i <= i_2; ++i) { t = z[i + k2 * z_dim1] + u2 * z[i + k1 * z_dim1] + u3 * z[i + k * z_dim1]; z[i + k2 * z_dim1] += t * v1; z[i + k1 * z_dim1] += t * v2; z[i + k * z_dim1] += t * v3; /* L235: */ } /* .......... ZERO B(K+1,K) .......... */ L240: s = (d_1 = b[k1 + k1 * b_dim1], abs(d_1)) + (d_2 = b[k1 + k * b_dim1], abs(d_2)); if (s == 0.) { goto L260; } u1 = b[k1 + k1 * b_dim1] / s; u2 = b[k1 + k * b_dim1] / s; d_1 = sqrt(u1 * u1 + u2 * u2); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; u2 = v2 / v1; i_2 = ll; for (i = lor1; i <= i_2; ++i) { t = a[i + k1 * a_dim1] + u2 * a[i + k * a_dim1]; a[i + k1 * a_dim1] += t * v1; a[i + k * a_dim1] += t * v2; t = b[i + k1 * b_dim1] + u2 * b[i + k * b_dim1]; b[i + k1 * b_dim1] += t * v1; b[i + k * b_dim1] += t * v2; /* L250: */ } b[k1 + k * b_dim1] = 0.; if (! (*matz)) { goto L260; } i_2 = *n; for (i = 1; i <= i_2; ++i) { t = z[i + k1 * z_dim1] + u2 * z[i + k * z_dim1]; z[i + k1 * z_dim1] += t * v1; z[i + k * z_dim1] += t * v2; /* L255: */ } L260: ;} /* .......... END QZ STEP .......... */ goto L70; /* .......... SET ERROR -- ALL EIGENVALUES HAVE NOT */ /* CONVERGED AFTER 30*N ITERATIONS .......... */ L1000: *ierr = en; /* .......... SAVE EPSB FOR USE BY QZVAL AND QZVEC .......... */ L1001: if (*n > 1) { b[*n + b_dim1] = epsb; } return 0; } /* qzit_ */ /* Subroutine */ int qzval_(nm, n, a, b, alfr, alfi, beta, matz, z) integer *nm, *n; doublereal *a, *b, *alfr, *alfi, *beta; logical *matz; doublereal *z; { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, z_dim1, z_offset, i_1, i_2; doublereal d_1, d_2, d_3, d_4; /* Builtin functions */ double sqrt(), d_sign(); /* Local variables */ static doublereal epsb, c, d, e; static integer i, j; static doublereal r, s, t, a1, a2, u1, u2, v1, v2, a11, a12, a21, a22, b11, b12, b22, di, ei; static integer na; static doublereal an, bn; static integer en; static doublereal cq, dr; static integer nn; static doublereal cz, ti, tr, a1i, a2i, a11i, a12i, a22i, a11r, a12r, a22r, sqi, ssi; static integer isw; static doublereal sqr, szi, ssr, szr; /* Parameter adjustments */ a_dim1 = *nm; a_offset = a_dim1 + 1; a -= a_offset; b_dim1 = *nm; b_offset = b_dim1 + 1; b -= b_offset; --alfr; --alfi; --beta; z_dim1 = *nm; z_offset = z_dim1 + 1; z -= z_offset; /* Function Body */ /* THIS SUBROUTINE IS THE THIRD STEP OF THE QZ ALGORITHM */ /* FOR SOLVING GENERALIZED MATRIX EIGENVALUE PROBLEMS, */ /* SIAM J. NUMER. ANAL. 10, 241-256(1973) BY MOLER AND STEWART. */ /* THIS SUBROUTINE ACCEPTS A PAIR OF REAL MATRICES, ONE OF THEM */ /* IN QUASI-TRIANGULAR FORM AND THE OTHER IN UPPER TRIANGULAR FORM. */ /* IT REDUCES THE QUASI-TRIANGULAR MATRIX FURTHER, SO THAT ANY */ /* REMAINING 2-BY-2 BLOCKS CORRESPOND TO PAIRS OF COMPLEX */ /* EIGENVALUES, AND RETURNS QUANTITIES WHOSE RATIOS GIVE THE */ /* GENERALIZED EIGENVALUES. IT IS USUALLY PRECEDED BY QZHES */ /* AND QZIT AND MAY BE FOLLOWED BY QZVEC. */ /* ON INPUT */ /* NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL */ /* ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM */ /* DIMENSION STATEMENT. */ /* N IS THE ORDER OF THE MATRICES. */ /* A CONTAINS A REAL UPPER QUASI-TRIANGULAR MATRIX. */ /* B CONTAINS A REAL UPPER TRIANGULAR MATRIX. IN ADDITION, */ /* LOCATION B(N,1) CONTAINS THE TOLERANCE QUANTITY (EPSB) */ /* COMPUTED AND SAVED IN QZIT. */ /* MATZ SHOULD BE SET TO .TRUE. IF THE RIGHT HAND TRANSFORMATIONS */ /* ARE TO BE ACCUMULATED FOR LATER USE IN COMPUTING */ /* EIGENVECTORS, AND TO .FALSE. OTHERWISE. */ /* Z CONTAINS, IF MATZ HAS BEEN SET TO .TRUE., THE */ /* TRANSFORMATION MATRIX PRODUCED IN THE REDUCTIONS BY QZHES */ /* AND QZIT, IF PERFORMED, OR ELSE THE IDENTITY MATRIX. */ /* IF MATZ HAS BEEN SET TO .FALSE., Z IS NOT REFERENCED. */ /* ON OUTPUT */ /* A HAS BEEN REDUCED FURTHER TO A QUASI-TRIANGULAR MATRIX */ /* IN WHICH ALL NONZERO SUBDIAGONAL ELEMENTS CORRESPOND TO */ /* PAIRS OF COMPLEX EIGENVALUES. */ /* B IS STILL IN UPPER TRIANGULAR FORM, ALTHOUGH ITS ELEMENTS */ /* HAVE BEEN ALTERED. B(N,1) IS UNALTERED. */ /* ALFR AND ALFI CONTAIN THE REAL AND IMAGINARY PARTS OF THE */ /* DIAGONAL ELEMENTS OF THE TRIANGULAR MATRIX THAT WOULD BE */ /* OBTAINED IF A WERE REDUCED COMPLETELY TO TRIANGULAR FORM */ /* BY UNITARY TRANSFORMATIONS. NON-ZERO VALUES OF ALFI OCCUR */ /* IN PAIRS, THE FIRST MEMBER POSITIVE AND THE SECOND NEGATIVE. */ /* BETA CONTAINS THE DIAGONAL ELEMENTS OF THE CORRESPONDING B, */ /* NORMALIZED TO BE REAL AND NON-NEGATIVE. THE GENERALIZED */ /* EIGENVALUES ARE THEN THE RATIOS ((ALFR+I*ALFI)/BETA). */ /* Z CONTAINS THE PRODUCT OF THE RIGHT HAND TRANSFORMATIONS */ /* (FOR ALL THREE STEPS) IF MATZ HAS BEEN SET TO .TRUE. */ /* QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, */ /* MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY */ /* THIS VERSION DATED AUGUST 1983. */ /* ------------------------------------------------------------------ */ epsb = b[*n + b_dim1]; isw = 1; /* .......... FIND EIGENVALUES OF QUASI-TRIANGULAR MATRICES. */ /* FOR EN=N STEP -1 UNTIL 1 DO -- .......... */ i_1 = *n; for (nn = 1; nn <= i_1; ++nn) { en = *n + 1 - nn; na = en - 1; if (isw == 2) { goto L505; } if (en == 1) { goto L410; } if (a[en + na * a_dim1] != 0.) { goto L420; } /* .......... 1-BY-1 BLOCK, ONE REAL ROOT .......... */ L410: alfr[en] = a[en + en * a_dim1]; if (b[en + en * b_dim1] < 0.) { alfr[en] = -alfr[en]; } beta[en] = (d_1 = b[en + en * b_dim1], abs(d_1)); alfi[en] = 0.; goto L510; /* .......... 2-BY-2 BLOCK .......... */ L420: if ((d_1 = b[na + na * b_dim1], abs(d_1)) <= epsb) { goto L455; } if ((d_1 = b[en + en * b_dim1], abs(d_1)) > epsb) { goto L430; } a1 = a[en + en * a_dim1]; a2 = a[en + na * a_dim1]; bn = 0.; goto L435; L430: an = (d_1 = a[na + na * a_dim1], abs(d_1)) + (d_2 = a[na + en * a_dim1], abs(d_2)) + (d_3 = a[en + na * a_dim1], abs(d_3)) + ( d_4 = a[en + en * a_dim1], abs(d_4)); bn = (d_1 = b[na + na * b_dim1], abs(d_1)) + (d_2 = b[na + en * b_dim1], abs(d_2)) + (d_3 = b[en + en * b_dim1], abs(d_3)); a11 = a[na + na * a_dim1] / an; a12 = a[na + en * a_dim1] / an; a21 = a[en + na * a_dim1] / an; a22 = a[en + en * a_dim1] / an; b11 = b[na + na * b_dim1] / bn; b12 = b[na + en * b_dim1] / bn; b22 = b[en + en * b_dim1] / bn; e = a11 / b11; ei = a22 / b22; s = a21 / (b11 * b22); t = (a22 - e * b22) / b22; if (abs(e) <= abs(ei)) { goto L431; } e = ei; t = (a11 - e * b11) / b11; L431: c = (t - s * b12) * .5; d = c * c + s * (a12 - e * b12); if (d < 0.) { goto L480; } /* .......... TWO REAL ROOTS. */ /* ZERO BOTH A(EN,NA) AND B(EN,NA) .......... */ d_1 = sqrt(d); e += c + d_sign(&d_1, &c); a11 -= e * b11; a12 -= e * b12; a22 -= e * b22; if (abs(a11) + abs(a12) < abs(a21) + abs(a22)) { goto L432; } a1 = a12; a2 = a11; goto L435; L432: a1 = a22; a2 = a21; /* .......... CHOOSE AND APPLY REAL Z .......... */ L435: s = abs(a1) + abs(a2); u1 = a1 / s; u2 = a2 / s; d_1 = sqrt(u1 * u1 + u2 * u2); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; u2 = v2 / v1; i_2 = en; for (i = 1; i <= i_2; ++i) { t = a[i + en * a_dim1] + u2 * a[i + na * a_dim1]; a[i + en * a_dim1] += t * v1; a[i + na * a_dim1] += t * v2; t = b[i + en * b_dim1] + u2 * b[i + na * b_dim1]; b[i + en * b_dim1] += t * v1; b[i + na * b_dim1] += t * v2; /* L440: */ } if (! (*matz)) { goto L450; } i_2 = *n; for (i = 1; i <= i_2; ++i) { t = z[i + en * z_dim1] + u2 * z[i + na * z_dim1]; z[i + en * z_dim1] += t * v1; z[i + na * z_dim1] += t * v2; /* L445: */ } L450: if (bn == 0.) { goto L475; } if (an < abs(e) * bn) { goto L455; } a1 = b[na + na * b_dim1]; a2 = b[en + na * b_dim1]; goto L460; L455: a1 = a[na + na * a_dim1]; a2 = a[en + na * a_dim1]; /* .......... CHOOSE AND APPLY REAL Q .......... */ L460: s = abs(a1) + abs(a2); if (s == 0.) { goto L475; } u1 = a1 / s; u2 = a2 / s; d_1 = sqrt(u1 * u1 + u2 * u2); r = d_sign(&d_1, &u1); v1 = -(u1 + r) / r; v2 = -u2 / r; u2 = v2 / v1; i_2 = *n; for (j = na; j <= i_2; ++j) { t = a[na + j * a_dim1] + u2 * a[en + j * a_dim1]; a[na + j * a_dim1] += t * v1; a[en + j * a_dim1] += t * v2; t = b[na + j * b_dim1] + u2 * b[en + j * b_dim1]; b[na + j * b_dim1] += t * v1; b[en + j * b_dim1] += t * v2; /* L470: */ } L475: a[en + na * a_dim1] = 0.; b[en + na * b_dim1] = 0.; alfr[na] = a[na + na * a_dim1]; alfr[en] = a[en + en * a_dim1]; if (b[na + na * b_dim1] < 0.) { alfr[na] = -alfr[na]; } if (b[en + en * b_dim1] < 0.) { alfr[en] = -alfr[en]; } beta[na] = (d_1 = b[na + na * b_dim1], abs(d_1)); beta[en] = (d_1 = b[en + en * b_dim1], abs(d_1)); alfi[en] = 0.; alfi[na] = 0.; goto L505; /* .......... TWO COMPLEX ROOTS .......... */ L480: e += c; ei = sqrt(-d); a11r = a11 - e * b11; a11i = ei * b11; a12r = a12 - e * b12; a12i = ei * b12; a22r = a22 - e * b22; a22i = ei * b22; if (abs(a11r) + abs(a11i) + abs(a12r) + abs(a12i) < abs(a21) + abs( a22r) + abs(a22i)) { goto L482; } a1 = a12r; a1i = a12i; a2 = -a11r; a2i = -a11i; goto L485; L482: a1 = a22r; a1i = a22i; a2 = -a21; a2i = 0.; /* .......... CHOOSE COMPLEX Z .......... */ L485: cz = sqrt(a1 * a1 + a1i * a1i); if (cz == 0.) { goto L487; } szr = (a1 * a2 + a1i * a2i) / cz; szi = (a1 * a2i - a1i * a2) / cz; r = sqrt(cz * cz + szr * szr + szi * szi); cz /= r; szr /= r; szi /= r; goto L490; L487: szr = 1.; szi = 0.; L490: if (an < (abs(e) + ei) * bn) { goto L492; } a1 = cz * b11 + szr * b12; a1i = szi * b12; a2 = szr * b22; a2i = szi * b22; goto L495; L492: a1 = cz * a11 + szr * a12; a1i = szi * a12; a2 = cz * a21 + szr * a22; a2i = szi * a22; /* .......... CHOOSE COMPLEX Q .......... */ L495: cq = sqrt(a1 * a1 + a1i * a1i); if (cq == 0.) { goto L497; } sqr = (a1 * a2 + a1i * a2i) / cq; sqi = (a1 * a2i - a1i * a2) / cq; r = sqrt(cq * cq + sqr * sqr + sqi * sqi); cq /= r; sqr /= r; sqi /= r; goto L500; L497: sqr = 1.; sqi = 0.; /* .......... COMPUTE DIAGONAL ELEMENTS THAT WOULD RESULT */ /* IF TRANSFORMATIONS WERE APPLIED .......... */ L500: ssr = sqr * szr + sqi * szi; ssi = sqr * szi - sqi * szr; i = 1; tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21 + ssr * a22; ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22; dr = cq * cz * b11 + cq * szr * b12 + ssr * b22; di = cq * szi * b12 + ssi * b22; goto L503; L502: i = 2; tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21 + cq * cz * a22; ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21; dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22; di = -ssi * b11 - sqi * cz * b12; L503: t = ti * dr - tr * di; j = na; if (t < 0.) { j = en; } r = sqrt(dr * dr + di * di); beta[j] = bn * r; alfr[j] = an * (tr * dr + ti * di) / r; alfi[j] = an * t / r; if (i == 1) { goto L502; } L505: isw = 3 - isw; L510: ;} b[*n + b_dim1] = epsb; return 0; } /* qzval_ */ doublereal epslon_(x) doublereal *x; { /* System generated locals */ doublereal ret_val, d_1; /* Local variables */ static doublereal a, b, c, eps; /* ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. */ /* THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS */ /* SATISFYING THE FOLLOWING TWO ASSUMPTIONS, */ /* 1. THE BASE USED IN REPRESENTING FLOATING POINT */ /* NUMBERS IS NOT A POWER OF THREE. */ /* 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO */ /* THE ACCURACY USED IN FLOATING POINT VARIABLES */ /* THAT ARE STORED IN MEMORY. */ /* THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO */ /* FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING */ /* ASSUMPTION 2. */ /* UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, */ /* A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, */ /* B HAS A ZERO FOR ITS LAST BIT OR DIGIT, */ /* C IS NOT EXACTLY EQUAL TO ONE, */ /* EPS MEASURES THE SEPARATION OF 1.0 FROM */ /* THE NEXT LARGER FLOATING POINT NUMBER. */ /* THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED */ /* ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD. */ /* THIS VERSION DATED 4/6/83. */ a = 1.3333333333333333; L10: b = a - 1.; c = b + b + b; eps = (d_1 = c - 1., abs(d_1)); if (eps == 0.) { goto L10; } ret_val = eps * abs(*x); return ret_val; } /* epslon_ */ /* ************************ */ /* * Routines from BLAS * */ /* ************************ */ doublereal ddot_(n, dx, incx, dy, incy) integer *n; doublereal *dx; integer *incx; doublereal *dy; integer *incy; { /* System generated locals */ integer i_1; doublereal ret_val; /* Local variables */ static integer i, m; static doublereal dtemp; static integer ix, iy, mp1; /* Parameter adjustments */ --dx; --dy; /* Function Body */ /* FORMS THE DOT PRODUCT OF TWO VECTORS. */ /* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */ /* JACK DONGARRA, LINPACK, 3/11/78. */ ret_val = 0.; dtemp = 0.; if (*n <= 0) { return ret_val; } if (*incx == 1 && *incy == 1) { goto L20; } /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */ /* NOT EQUAL TO 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i_1 = *n; for (i = 1; i <= i_1; ++i) { dtemp += dx[ix] * dy[iy]; ix += *incx; iy += *incy; /* L10: */ } ret_val = dtemp; return ret_val; /* CODE FOR BOTH INCREMENTS EQUAL TO 1 */ /* CLEAN-UP LOOP */ L20: m = *n % 5; if (m == 0) { goto L40; } i_1 = m; for (i = 1; i <= i_1; ++i) { dtemp += dx[i] * dy[i]; /* L30: */ } if (*n < 5) { goto L60; } L40: mp1 = m + 1; i_1 = *n; for (i = mp1; i <= i_1; i += 5) { dtemp = dtemp + dx[i] * dy[i] + dx[i + 1] * dy[i + 1] + dx[i + 2] * dy[i + 2] + dx[i + 3] * dy[i + 3] + dx[i + 4] * dy[i + 4]; /* L50: */ } L60: ret_val = dtemp; return ret_val; } /* ddot_ */ integer idamax_(n, dx, incx) integer *n; doublereal *dx; integer *incx; { /* System generated locals */ integer ret_val, i_1; doublereal d_1; /* Local variables */ static doublereal dmax_; static integer i, ix; /* Parameter adjustments */ --dx; /* Function Body */ /* FINDS THE INDEX OF ELEMENT HAVING MAX. ABSOLUTE VALUE. */ /* JACK DONGARRA, LINPACK, 3/11/78. */ ret_val = 0; if (*n < 1) { return ret_val; } ret_val = 1; if (*n == 1) { return ret_val; } if (*incx == 1) { goto L20; } /* CODE FOR INCREMENT NOT EQUAL TO 1 */ ix = 1; dmax_ = abs(dx[1]); ix += *incx; i_1 = *n; for (i = 2; i <= i_1; ++i) { if ((d_1 = dx[ix], abs(d_1)) <= dmax_) { goto L5; } ret_val = i; dmax_ = (d_1 = dx[ix], abs(d_1)); L5: ix += *incx; /* L10: */ } return ret_val; /* CODE FOR INCREMENT EQUAL TO 1 */ L20: dmax_ = abs(dx[1]); i_1 = *n; for (i = 2; i <= i_1; ++i) { if ((d_1 = dx[i], abs(d_1)) <= dmax_) { goto L30; } ret_val = i; dmax_ = (d_1 = dx[i], abs(d_1)); L30: ;} return ret_val; } /* idamax_ */ logical lsame_(ca, cb, ca_len, cb_len) char *ca, *cb; ftnlen ca_len; ftnlen cb_len; { /* System generated locals */ logical ret_val; /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* LSAME tests if CA is the same letter as CB regardless of case. */ /* CB is assumed to be an upper case letter. LSAME returns .TRUE. if */ /* CA is either the same as CB or the equivalent lower case letter. */ /* N.B. This version of the routine is only correct for ASCII code. */ /* Installers must modify the routine for other character-codes. */ /* For EBCDIC systems the constant IOFF must be changed to -64. */ /* For CDC systems using 6-12 bit representations, the system- */ /* specific code in comments must be activated. */ /* Parameters */ /* ========== */ /* CA - CHARACTER*1 */ /* CB - CHARACTER*1 */ /* On entry, CA and CB specify characters to be compared. */ /* Unchanged on exit. */ /* Auxiliary routine for Level 2 Blas. */ /* -- Written on 20-July-1986 */ /* Richard Hanson, Sandia National Labs. */ /* Jeremy Du Croz, Nag Central Office. */ /* .. Parameters .. */ /* .. Intrinsic Functions .. */ /* .. Executable Statements .. */ /* Test if the characters are equal */ ret_val = *ca == *cb; /* Now test for equivalence */ if (! ret_val) { ret_val = *ca - 32 == *cb; } return ret_val; /* The following comments contain code for CDC systems using 6-12 bit */ /* representations. */ /* .. Parameters .. */ /* INTEGER ICIRFX */ /* PARAMETER ( ICIRFX=62 ) */ /* .. Scalar Arguments .. */ /* CHARACTER*1 CB */ /* .. Array Arguments .. */ /* CHARACTER*1 CA(*) */ /* .. Local Scalars .. */ /* INTEGER IVAL */ /* .. Intrinsic Functions .. */ /* INTRINSIC ICHAR, CHAR */ /* .. Executable Statements .. */ /* See if the first character in string CA equals string CB. */ /* LSAME = CA(1) .EQ. CB .AND. CA(1) .NE. CHAR(ICIRFX) */ /* IF (LSAME) RETURN */ /* The characters are not identical. Now check them for equivalence. */ /* Look for the 'escape' character, circumflex, followed by the */ /* letter. */ /* IVAL = ICHAR(CA(2)) */ /* IF (IVAL.GE.ICHAR('A') .AND. IVAL.LE.ICHAR('Z')) THEN */ /* LSAME = CA(1) .EQ. CHAR(ICIRFX) .AND. CA(2) .EQ. CB */ /* END IF */ /* RETURN */ /* End of LSAME. */ } /* lsame_ */ /* Subroutine */ int xerbla_(srname, info, srname_len) char *srname; integer *info; ftnlen srname_len; { /* Format strings */ static char fmt_99999[] = "(\002 ** On entry to \002,a6,\002 parameter n\ umber \002,i2,\002 had an illegal value\002)"; /* Builtin functions */ integer s_wsfe(), do_fio(), e_wsfe(); /* Subroutine */ int s_stop(); /* Fortran I/O blocks */ static cilist io__211 = { 0, 6, 0, fmt_99999, 0 }; /* .. Scalar Arguments .. */ /* .. */ /* Purpose */ /* ======= */ /* XERBLA is an error handler for the Level 2 BLAS routines. */ /* It is called by the Level 2 BLAS routines if an input parameter is */ /* invalid. */ /* Installers should consider modifying the STOP statement in order to */ /* call system-specific exception-handling facilities. */ /* Parameters */ /* ========== */ /* SRNAME - CHARACTER*6. */ /* On entry, SRNAME specifies the name of the routine which */ /* called XERBLA. */ /* INFO - INTEGER. */ /* On entry, INFO specifies the position of the invalid */ /* parameter in the parameter-list of the calling routine. */ /* Auxiliary routine for Level 2 Blas. */ /* Written on 20-July-1986. */ /* .. Executable Statements .. */ s_wsfe(&io__211); do_fio(&c__1, srname, 6L); do_fio(&c__1, (char *)&(*info), (ftnlen)sizeof(integer)); e_wsfe(); s_stop("", 0L); /* End of XERBLA. */ } /* xerbla_ */ /* Subroutine */ int daxpy_(n, da, dx, incx, dy, incy) integer *n; doublereal *da, *dx; integer *incx; doublereal *dy; integer *incy; { /* System generated locals */ integer i_1; /* Local variables */ static integer i, m, ix, iy, mp1; /* Parameter adjustments */ --dx; --dy; /* Function Body */ /* CONSTANT TIMES A VECTOR PLUS A VECTOR. */ /* USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. */ /* JACK DONGARRA, LINPACK, 3/11/78. */ if (*n <= 0) { return 0; } if (*da == 0.) { return 0; } if (*incx == 1 && *incy == 1) { goto L20; } /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS */ /* NOT EQUAL TO 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i_1 = *n; for (i = 1; i <= i_1; ++i) { dy[iy] += *da * dx[ix]; ix += *incx; iy += *incy; /* L10: */ } return 0; /* CODE FOR BOTH INCREMENTS EQUAL TO 1 */ /* CLEAN-UP LOOP */ L20: m = *n % 4; if (m == 0) { goto L40; } i_1 = m; for (i = 1; i <= i_1; ++i) { dy[i] += *da * dx[i]; /* L30: */ } if (*n < 4) { return 0; } L40: mp1 = m + 1; i_1 = *n; for (i = mp1; i <= i_1; i += 4) { dy[i] += *da * dx[i]; dy[i + 1] += *da * dx[i + 1]; dy[i + 2] += *da * dx[i + 2]; dy[i + 3] += *da * dx[i + 3]; /* L50: */ } return 0; } /* daxpy_ */ /* Subroutine */ int drot_(n, dx, incx, dy, incy, c, s) integer *n; doublereal *dx; integer *incx; doublereal *dy; integer *incy; doublereal *c, *s; { /* System generated locals */ integer i_1; /* Local variables */ static integer i; static doublereal dtemp; static integer ix, iy; /* Parameter adjustments */ --dx; --dy; /* Function Body */ /* APPLIES A PLANE ROTATION. */ /* JACK DONGARRA, LINPACK, 3/11/78. */ if (*n <= 0) { return 0; } if (*incx == 1 && *incy == 1) { goto L20; } /* CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS NOT EQUAL */ /* TO 1 */ ix = 1; iy = 1; if (*incx < 0) { ix = (-(*n) + 1) * *incx + 1; } if (*incy < 0) { iy = (-(*n) + 1) * *incy + 1; } i_1 = *n; for (i = 1; i <= i_1; ++i) { dtemp = *c * dx[ix] + *s * dy[iy]; dy[iy] = *c * dy[iy] - *s * dx[ix]; dx[ix] = dtemp; ix += *incx; iy += *incy; /* L10: */ } return 0; /* CODE FOR BOTH INCREMENTS EQUAL TO 1 */ L20: i_1 = *n; for (i = 1; i <= i_1; ++i) { dtemp = *c * dx[i] + *s * dy[i]; dy[i] = *c * dy[i] - *s * dx[i]; dx[i] = dtemp; /* L30: */ } return 0; } /* drot_ */ /* *************************** */ /* * Other useful routines * */ /* *************************** */ /* Subroutine */ int dgemc_(m, n, a, lda, b, ldb, trans) integer *m, *n; doublereal *a; integer *lda; doublereal *b; integer *ldb; logical *trans; { /* System generated locals */ integer a_dim1, a_offset, b_dim1, b_offset, i_1, i_2; /* Local variables */ static integer i, j, mm, mmp1; /* Parameter adjustments */ a_dim1 = *lda; a_offset = a_dim1 + 1; a -= a_offset; b_dim1 = *ldb; b_offset = b_dim1 + 1; b -= b_offset; /* Function Body */ /* This subroutine copies a double precision real */ /* M by N matrix stored in A to double precision real B. */ /* If TRANS is true, B is assigned A transpose. */ if (*trans) { i_1 = *n; for (j = 1; j <= i_1; ++j) { /* USES UNROLLED LOOPS */ /* from JACK DONGARRA, LINPACK, 3/11/78. */ mm = *m % 7; if (mm == 0) { goto L80; } i_2 = mm; for (i = 1; i <= i_2; ++i) { b[j + i * b_dim1] = a[i + j * a_dim1]; /* L70: */ } if (*m < 7) { goto L99; } L80: mmp1 = mm + 1; i_2 = *m; for (i = mmp1; i <= i_2; i += 7) { b[j + i * b_dim1] = a[i + j * a_dim1]; b[j + (i + 1) * b_dim1] = a[i + 1 + j * a_dim1]; b[j + (i + 2) * b_dim1] = a[i + 2 + j * a_dim1]; b[j + (i + 3) * b_dim1] = a[i + 3 + j * a_dim1]; b[j + (i + 4) * b_dim1] = a[i + 4 + j * a_dim1]; b[j + (i + 5) * b_dim1] = a[i + 5 + j * a_dim1]; b[j + (i + 6) * b_dim1] = a[i + 6 + j * a_dim1]; /* L90: */ } L99: /* L100: */ ;} } else { i_1 = *n; for (j = 1; j <= i_1; ++j) { /* USES UNROLLED LOOPS */ /* from JACK DONGARRA, LINPACK, 3/11/78. */ mm = *m % 7; if (mm == 0) { goto L180; } i_2 = mm; for (i = 1; i <= i_2; ++i) { b[i + j * b_dim1] = a[i + j * a_dim1]; /* L170: */ } if (*m < 7) { goto L199; } L180: mmp1 = mm + 1; i_2 = *m; for (i = mmp1; i <= i_2; i += 7) { b[i + j * b_dim1] = a[i + j * a_dim1]; b[i + 1 + j * b_dim1] = a[i + 1 + j * a_dim1]; b[i + 2 + j * b_dim1] = a[i + 2 + j * a_dim1]; b[i + 3 + j * b_dim1] = a[i + 3 + j * a_dim1]; b[i + 4 + j * b_dim1] = a[i + 4 + j * a_dim1]; b[i + 5 + j * b_dim1] = a[i + 5 + j * a_dim1]; b[i + 6 + j * b_dim1] = a[i + 6 + j * a_dim1]; /* L190: */ } L199: /* L200: */ ;} } return 0; } /* dgemc_ */ /* Subroutine */ int ezsvd_(x, ldx, n, p, s, e, u, ldu, v, ldv, work, job, info, tol) doublereal *x; integer *ldx, *n, *p; doublereal *s, *e, *u; integer *ldu; doublereal *v; integer *ldv; doublereal *work; integer *job, *info; doublereal *tol; { /* System generated locals */ integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset; /* Local variables */ static integer idbg, skip, iidir, ifull; extern /* Subroutine */ int ndsvd_(); static integer kount, kount1, kount2, limshf; static doublereal maxsin; static integer maxitr; /* Parameter adjustments */ x_dim1 = *ldx; x_offset = x_dim1 + 1; x -= x_offset; --s; --e; u_dim1 = *ldu; u_offset = u_dim1 + 1; u -= u_offset; v_dim1 = *ldv; v_offset = v_dim1 + 1; v -= v_offset; --work; /* Function Body */ /* new svd by J. Demmel, W. Kahan */ /* finds singular values of bidiagonal matrices with guaranteed high */ /* relative precision */ /* easy to use version of ndsvd ("hard to use" version, below) */ /* with defaults for some ndsvd parameters */ /* all parameters same as linpack dsvdc except for tol: */ /* tol = if positive, desired relative precision in singular values */ /* if negative, desired absolute precision in singular values */ /* (expressed as abs(tol) * sigma-max) */ /* (in both cases, abs(tol) should be less than 1 and */ /* greater than macheps) */ /* I have tested this software on a SUN 3 in double precision */ /* IEEE arithmetic with macheps about 2.2e-16 and tol=1e-14; */ /* In general I recommend tol 10-100 times larger than macheps. */ /* On the average it appears to be as fast or faster than dsvdc. */ /* I have seen it go 3.5 times faster and 2 times slower at the */ /* extremes. */ /* defaults for ndsvd parameters (see ndsvd for more description of */ /* these parameters) are: */ /* set to no debug output */ idbg = 0; /* use zero-shift normally */ ifull = 0; /* use normal bidiagonalization code */ skip = 0; /* choose chase direction normally */ iidir = 0; /* maximum 30 QR sweeps per singular value */ maxitr = 30; ndsvd_(&x[x_offset], ldx, n, p, &s[1], &e[1], &u[u_offset], ldu, &v[ v_offset], ldv, &work[1], job, info, &maxitr, tol, &idbg, &ifull, &kount, &kount1, &kount2, &skip, &limshf, &maxsin, &iidir); return 0; } /* ezsvd_ */ /* Subroutine */ int ndrotg_(f, g, cs, sn) doublereal *f, *g, *cs, *sn; { /* Builtin functions */ double sqrt(); /* Local variables */ static doublereal t, tt; /* faster version of drotg, except g unchanged on return */ /* cs, sn returned so that -sn*f+cs*g = 0 */ /* and returned f = cs*f + sn*g */ /* if g=0, then cs=1 and sn=0 (in case svd adds extra zero row */ /* to bidiagonal, this makes sure last row rotation is trivial) */ /* if f=0 and g.ne.0, then cs=0 and sn=1 without floating point work */ /* (in case s(i)=0 in svd so that bidiagonal deflates, this */ /* computes rotation without any floating point operations) */ if (*f == 0.) { if (*g == 0.) { /* this case needed in case extra zero row added in svd, so */ /* bottom rotation always trivial */ *cs = 1.; *sn = 0.; } else { /* this case needed for s(i)=0 in svd to compute rotation */ /* cheaply */ *cs = 0.; *sn = 1.; *f = *g; } } else { if (abs(*f) > abs(*g)) { t = *g / *f; tt = sqrt(t * t + 1.); *cs = 1. / tt; *sn = t * *cs; *f *= tt; } else { t = *f / *g; tt = sqrt(t * t + 1.); *sn = 1. / tt; *cs = t * *sn; *f = *g * tt; } } return 0; } /* ndrotg_ */ /* Subroutine */ int ndsvd_(x, ldx, n, p, s, e, u, ldu, v, ldv, work, job, info, maxitr, tol, idbg, ifull, kount, kount1, kount2, skip, limshf, maxsin, iidir) doublereal *x; integer *ldx, *n, *p; doublereal *s, *e, *u; integer *ldu; doublereal *v; integer *ldv; doublereal *work; integer *job, *info, *maxitr; doublereal *tol; integer *idbg, *ifull, *kount, *kount1, *kount2, *skip, *limshf; doublereal *maxsin; integer *iidir; { /* System generated locals */ integer x_dim1, x_offset, u_dim1, u_offset, v_dim1, v_offset, i_1, i_2, i_3; doublereal d_1, d_2, d_3, d_4; /* Builtin functions */ double d_sign(); integer s_wsle(), do_lio(), e_wsle(); /* Local variables */ static doublereal abse; extern /* Subroutine */ int sig22_(); static integer idir; static doublereal abss; extern doublereal ddot_(); static integer oldm, jobu; static doublereal cosl; static integer iter; static doublereal temp, smin, smax, cosr, sinl, sinr; extern /* Subroutine */ int prse_(), drot_(); static doublereal test; extern doublereal dnrm2_(); static integer nctp1, nrtp1; static doublereal f, g; static integer i, j, k, l, m; static doublereal t; extern /* Subroutine */ int dscal_(); static doublereal oldcs; static integer oldll, iisub; static doublereal shift, oldsn, sigmn; static integer minnp, maxit; static doublereal sminl; extern /* Subroutine */ int dswap_(), daxpy_(); static doublereal sigmx; static logical wantu, wantv; static doublereal gg, lambda; static integer oldacc; static doublereal cs; static integer ll, mm; static doublereal sm; static integer lu; static doublereal sn, mu; extern doublereal sigmin_(); static doublereal thresh; extern /* Subroutine */ int ndrotg_(); static integer lm1, lp1, lll, nct, ncu; static doublereal sll; static integer nrt; static doublereal emm1, smm1; /* Fortran I/O blocks */ static cilist io__269 = { 0, 6, 0, 0, 0 }; static cilist io__270 = { 0, 6, 0, 0, 0 }; static cilist io__272 = { 0, 6, 0, 0, 0 }; static cilist io__276 = { 0, 6, 0, 0, 0 }; static cilist io__277 = { 0, 6, 0, 0, 0 }; static cilist io__278 = { 0, 6, 0, 0, 0 }; static cilist io__279 = { 0, 6, 0, 0, 0 }; static cilist io__287 = { 0, 6, 0, 0, 0 }; static cilist io__290 = { 0, 6, 0, 0, 0 }; static cilist io__292 = { 0, 6, 0, 0, 0 }; static cilist io__295 = { 0, 6, 0, 0, 0 }; static cilist io__300 = { 0, 6, 0, 0, 0 }; static cilist io__301 = { 0, 6, 0, 0, 0 }; static cilist io__302 = { 0, 6, 0, 0, 0 }; static cilist io__303 = { 0, 6, 0, 0, 0 }; static cilist io__305 = { 0, 6, 0, 0, 0 }; static cilist io__306 = { 0, 6, 0, 0, 0 }; static cilist io__307 = { 0, 6, 0, 0, 0 }; static cilist io__308 = { 0, 6, 0, 0, 0 }; static cilist io__309 = { 0, 6, 0, 0, 0 }; static cilist io__318 = { 0, 6, 0, 0, 0 }; static cilist io__319 = { 0, 6, 0, 0, 0 }; static cilist io__320 = { 0, 6, 0, 0, 0 }; static cilist io__321 = { 0, 6, 0, 0, 0 }; static cilist io__322 = { 0, 6, 0, 0, 0 }; static cilist io__323 = { 0, 6, 0, 0, 0 }; static cilist io__324 = { 0, 6, 0, 0, 0 }; static cilist io__325 = { 0, 6, 0, 0, 0 }; static cilist io__326 = { 0, 6, 0, 0, 0 }; static cilist io__327 = { 0, 6, 0, 0, 0 }; static cilist io__328 = { 0, 6, 0, 0, 0 }; static cilist io__329 = { 0, 6, 0, 0, 0 }; static cilist io__330 = { 0, 6, 0, 0, 0 }; static cilist io__331 = { 0, 6, 0, 0, 0 }; static cilist io__332 = { 0, 6, 0, 0, 0 }; static cilist io__333 = { 0, 6, 0, 0, 0 }; /* Parameter adjustments */ x_dim1 = *ldx; x_offset = x_dim1 + 1; x -= x_offset; --s; --e; u_dim1 = *ldu; u_offset = u_dim1 + 1; u -= u_offset; v_dim1 = *ldv; v_offset = v_dim1 + 1; v -= v_offset; --work; /* Function Body */ /* LINPACK SVD modified by: */ /* James Demmel W. Kahan */ /* Courant Institute Computer Science Dept. */ /* demmel@acf8.nyu.edu U.C. Berkeley */ /* modified version designed to guarantee relative accuracy of */ /* all singular values of intermediate bidiagonal form */ /* extra input/output parameters in addition to those from LINPACK SVD: */ /* extra input paramters: */ /* tol = if positive, desired relative precision in singular values */ /* if negative, desired absolute precision in singular values */ /* (expressed as abs(tol) * sigma-max) */ /* (abs(tol) should be less than 1 and greater than macheps) */ /* idbg = 0 for no debug output (normal setting) */ /* = 1 convergence, shift decisions (written to standard output) */ /* = 2 for above plus before, after qr */ /* ifull= 0 if decision to use zero-shift set normally (normal setting) */ /* = 1 if always set to nonzero-shift */ /* = 2 if always set to zero-shift */ /* skip =-1 means standard code but do all work of bidiagonalization */ /* (even if input bidiagonal) */ /* 0 means standard code (normal setting) */ /* 1 means assume x is bidiagonal, and skip bidiagonalization */ /* entirely */ /* (skip used for timing tests) */ /* iidir = 0 if idir (chase direction) chosen normally */ /* 1 if idir=1 (chase top to bottom) always */ /* 2 if idir=2 (chase bottom to top) always */ /* extra output parameters: */ /* kount =number of qr sweeps taken */ /* kount1=number of passes through inner loop of full qr */ /* kount2=number of passes through inner loop of zero-shift qr */ /* limshf = number of times the shift was greater than its threshold */ /* (nct*smin) and had to be decreased */ /* maxsin = maximum sin in inner loop of zero-shift */ /* new version designed to be robust with respect to over/underflow */ /* have fast inner loop when shift is zero, */ /* guarantee relative accuracy of all singular values */ /* dsvdc is a subroutine to reduce a double precision nxp matrix x */ /* by orthogonal transformations u and v to diagonal form. the */ /* diagonal elements s(i) are the singular values of x. the */ /* columns of u are the corresponding left singular vectors, */ /* and the columns of v the right singular vectors. */ /* on entry */ /* x double precision(ldx,p), where ldx.ge.n. */ /* x contains the matrix whose singular value */ /* decomposition is to be computed. x is */ /* destroyed by dsvdc. */ /* ldx integer. */ /* ldx is the leading dimension of the array x. */ /* n integer. */ /* n is the number of rows of the matrix x. */ /* p integer. */ /* p is the number of columns of the matrix x. */ /* ldu integer. */ /* ldu is the leading dimension of the array u. */ /* (see below). */ /* ldv integer. */ /* ldv is the leading dimension of the array v. */ /* (see below). */ /* work double precision(n). */ /* work is a scratch array. */ /* job integer. */ /* job controls the computation of the singular */ /* vectors. it has the decimal expansion ab */ /* with the following meaning */ /* a.eq.0 do not compute the left singular */ /* vectors. */ /* a.eq.1 return the n left singular vectors */ /* in u. */ /* a.ge.2 return the first min(n,p) singular */ /* vectors in u. */ /* external drot */ /* blas daxpy,ddot,dscal,dswap,dnrm2,drotg */ /* fortran dabs,dmax1,max0,min0,mod,dsqrt,dsign */ /* prse,ndrotg,sig22,sndrtg,sigmin */ /* internal variables */ /* new variables */ /* double precision sg1,sg2 */ /* set the maximum number of iterations. */ /* maxit = 30 */ *kount = 0; *kount1 = 0; *kount2 = 0; *limshf = 0; *maxsin = (float)0.; /* determine what is to be computed. */ wantu = FALSE_; wantv = FALSE_; jobu = *job % 100 / 10; ncu = *n; if (jobu > 1) { ncu = min(*n,*p); } if (jobu != 0) { wantu = TRUE_; } if (*job % 10 != 0) { wantv = TRUE_; } /* reduce x to bidiagonal form, storing the diagonal elements */ /* in s and the super-diagonal elements in e. */ *info = 0; /* Computing MAX */ i_1 = *n - 1; nct = min(*p,i_1); /* Computing MAX */ /* Computing MAX */ i_3 = *p - 2; i_1 = 0, i_2 = min(*n,i_3); nrt = max(i_2,i_1); lu = max(nct,nrt); if (*skip <= 0) { if (lu < 1) { goto L170; } i_1 = lu; for (l = 1; l <= i_1; ++l) { lp1 = l + 1; if (l > nct) { goto L20; } /* compute the transformation for the l-th column and */ /* place the l-th diagonal in s(l). */ i_2 = *n - l + 1; s[l] = dnrm2_(&i_2, &x[l + l * x_dim1], &c__1); if (s[l] == 0. && *skip == 0) { goto L10; } if (x[l + l * x_dim1] != 0.) { s[l] = d_sign(&s[l], &x[l + l * x_dim1]); } i_2 = *n - l + 1; d_1 = 1. / s[l]; dscal_(&i_2, &d_1, &x[l + l * x_dim1], &c__1); x[l + l * x_dim1] += 1.; L10: s[l] = -s[l]; L20: if (*p < lp1) { goto L50; } i_2 = *p; for (j = lp1; j <= i_2; ++j) { if (l > nct) { goto L30; } if (s[l] == 0. && *skip == 0) { goto L30; } /* apply the transformation. */ i_3 = *n - l + 1; t = -ddot_(&i_3, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1] , &c__1) / x[l + l * x_dim1]; i_3 = *n - l + 1; daxpy_(&i_3, &t, &x[l + l * x_dim1], &c__1, &x[l + j * x_dim1] , &c__1); L30: /* place the l-th row of x into e for the */ /* subsequent calculation of the row transformation. */ e[j] = x[l + j * x_dim1]; /* L40: */ } L50: if (! wantu || l > nct) { goto L70; } /* place the transformation in u for subsequent back */ /* multiplication. */ i_2 = *n; for (i = l; i <= i_2; ++i) { u[i + l * u_dim1] = x[i + l * x_dim1]; /* L60: */ } L70: if (l > nrt) { goto L150; } /* compute the l-th row transformation and place the */ /* l-th super-diagonal in e(l). */ i_2 = *p - l; e[l] = dnrm2_(&i_2, &e[lp1], &c__1); if (e[l] == 0. && *skip == 0) { goto L80; } if (e[lp1] != 0.) { e[l] = d_sign(&e[l], &e[lp1]); } i_2 = *p - l; d_1 = 1. / e[l]; dscal_(&i_2, &d_1, &e[lp1], &c__1); e[lp1] += 1.; L80: e[l] = -e[l]; if (lp1 > *n || e[l] == 0. && *skip == 0) { goto L120; } /* apply the transformation. */ i_2 = *n; for (i = lp1; i <= i_2; ++i) { work[i] = 0.; /* L90: */ } i_2 = *p; for (j = lp1; j <= i_2; ++j) { i_3 = *n - l; daxpy_(&i_3, &e[j], &x[lp1 + j * x_dim1], &c__1, &work[lp1], & c__1); /* L100: */ } i_2 = *p; for (j = lp1; j <= i_2; ++j) { i_3 = *n - l; d_1 = -e[j] / e[lp1]; daxpy_(&i_3, &d_1, &work[lp1], &c__1, &x[lp1 + j * x_dim1], & c__1); /* L110: */ } L120: if (! wantv) { goto L140; } /* place the transformation in v for subsequent */ /* back multiplication. */ i_2 = *p; for (i = lp1; i <= i_2; ++i) { v[i + l * v_dim1] = e[i]; /* L130: */ } L140: L150: /* L160: */ ;} L170: ;} /* set up the final bidiagonal matrix or order m. */ /* Computing MAX */ i_1 = *p, i_2 = *n + 1; m = min(i_2,i_1); nctp1 = nct + 1; nrtp1 = nrt + 1; if (*skip <= 0) { if (nct < *p) { s[nctp1] = x[nctp1 + nctp1 * x_dim1]; } if (*n < m) { s[m] = 0.; } if (nrtp1 < m) { e[nrtp1] = x[nrtp1 + m * x_dim1]; } e[m] = 0.; /* if required, generate u. */ if (! wantu) { goto L300; } if (ncu < nctp1) { goto L200; } i_1 = ncu; for (j = nctp1; j <= i_1; ++j) { i_2 = *n; for (i = 1; i <= i_2; ++i) { u[i + j * u_dim1] = 0.; /* L180: */ } u[j + j * u_dim1] = 1.; /* L190: */ } L200: if (nct < 1) { goto L290; } i_1 = nct; for (ll = 1; ll <= i_1; ++ll) { l = nct - ll + 1; if (s[l] == 0.) { goto L250; } lp1 = l + 1; if (ncu < lp1) { goto L220; } i_2 = ncu; for (j = lp1; j <= i_2; ++j) { i_3 = *n - l + 1; t = -ddot_(&i_3, &u[l + l * u_dim1], &c__1, &u[l + j * u_dim1] , &c__1) / u[l + l * u_dim1]; i_3 = *n - l + 1; daxpy_(&i_3, &t, &u[l + l * u_dim1], &c__1, &u[l + j * u_dim1] , &c__1); /* L210: */ } L220: i_2 = *n - l + 1; dscal_(&i_2, &c_b338, &u[l + l * u_dim1], &c__1); u[l + l * u_dim1] += 1.; lm1 = l - 1; if (lm1 < 1) { goto L240; } i_2 = lm1; for (i = 1; i <= i_2; ++i) { u[i + l * u_dim1] = 0.; /* L230: */ } L240: goto L270; L250: i_2 = *n; for (i = 1; i <= i_2; ++i) { u[i + l * u_dim1] = 0.; /* L260: */ } u[l + l * u_dim1] = 1.; L270: /* L280: */ ;} L290: L300: /* if it is required, generate v. */ if (! wantv) { goto L350; } i_1 = *p; for (ll = 1; ll <= i_1; ++ll) { l = *p - ll + 1; lp1 = l + 1; if (l > nrt) { goto L320; } if (e[l] == 0.) { goto L320; } i_2 = *p; for (j = lp1; j <= i_2; ++j) { i_3 = *p - l; t = -ddot_(&i_3, &v[lp1 + l * v_dim1], &c__1, &v[lp1 + j * v_dim1], &c__1) / v[lp1 + l * v_dim1]; i_3 = *p - l; daxpy_(&i_3, &t, &v[lp1 + l * v_dim1], &c__1, &v[lp1 + j * v_dim1], &c__1); /* L310: */ } L320: i_2 = *p; for (i = 1; i <= i_2; ++i) { v[i + l * v_dim1] = 0.; /* L330: */ } v[l + l * v_dim1] = 1.; /* L340: */ } L350: ;} if (*skip == 1) { /* set up s,e,u,v assuming x bidiagonal on input */ minnp = min(*n,*p); i_1 = minnp; for (i = 1; i <= i_1; ++i) { s[i] = x[i + i * x_dim1]; if (i < *p) { e[i] = x[i + (i + 1) * x_dim1]; } /* L351: */ } if (*n < *p) { s[*n + 1] = (float)0.; } e[m] = 0.; if (wantu) { i_1 = ncu; for (j = 1; j <= i_1; ++j) { i_2 = *n; for (i = 1; i <= i_2; ++i) { u[i + j * u_dim1] = (float)0.; /* L353: */ } u[j + j * u_dim1] = 1.; /* L352: */ } } if (wantv) { i_1 = *p; for (j = 1; j <= i_1; ++j) { i_2 = *p; for (i = 1; i <= i_2; ++i) { v[i + j * v_dim1] = (float)0.; /* L355: */ } v[j + j * v_dim1] = 1.; /* L354: */ } } } /* main iteration loop for the singular values. */ /* convert maxit to bound on total number of passes through */ /* inner loops of qr iteration (half number of rotations) */ maxit = *maxitr * m * m / 2; iter = 0; oldll = -1; oldm = -1; oldacc = -1; if (*tol > 0.) { /* relative accuracy desired */ thresh = 0.; } else { /* absolute accuracy desired */ smax = (d_1 = s[m], abs(d_1)); i_1 = m - 1; for (i = 1; i <= i_1; ++i) { /* Computing MAX */ d_3 = smax, d_4 = (d_1 = s[i], abs(d_1)), d_3 = max(d_4,d_3), d_4 = (d_2 = e[i], abs(d_2)); smax = max(d_4,d_3); /* L1111: */ } thresh = abs(*tol) * smax; } mm = m; /* begin loop */ L999: if (*idbg > 0) { s_wsle(&io__269); do_lio(&c__9, &c__1, "top of loop", 11L); e_wsle(); s_wsle(&io__270); do_lio(&c__9, &c__1, "oldll,oldm,oldacc,m,iter,maxit,ifull,thresh=", 44L); do_lio(&c__3, &c__1, (char *)&oldll, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&oldm, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&oldacc, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&m, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&iter, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&maxit, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&(*ifull), (ftnlen)sizeof(integer)); do_lio(&c__5, &c__1, (char *)&thresh, (ftnlen)sizeof(doublereal)); e_wsle(); prse_(&c__1, &mm, n, p, &s[1], &e[1]); } /* check for being done */ if (m == 1) { goto L998; } /* check number of iterations */ if (iter >= maxit) { goto L997; } /* compute minimum s(i) and max of all s(i),e(i) */ if (*tol <= 0. && (d_1 = s[m], abs(d_1)) <= thresh) { s[m] = (float)0.; } smax = (d_1 = s[m], abs(d_1)); smin = smax; /* reset convergence threshold if starting new part of matrix */ if (m <= oldll && *tol > 0.) { thresh = 0.; } if (*idbg > 0) { s_wsle(&io__272); do_lio(&c__9, &c__1, "thresh=", 7L); do_lio(&c__5, &c__1, (char *)&thresh, (ftnlen)sizeof(doublereal)); e_wsle(); } i_1 = m; for (lll = 1; lll <= i_1; ++lll) { ll = m - lll; if (ll == 0) { goto L1003; } if (*tol <= 0. && (d_1 = s[ll], abs(d_1)) <= thresh) { s[ll] = (float)0.; } if ((d_1 = e[ll], abs(d_1)) <= thresh) { goto L1002; } abss = (d_1 = s[ll], abs(d_1)); abse = (d_1 = e[ll], abs(d_1)); smin = min(smin,abss); /* Computing MAX */ d_1 = smax, d_1 = max(abss,d_1); smax = max(abse,d_1); /* L1001: */ } L1002: e[ll] = 0.; /* matrix splits since e(ll)=0 */ if (ll == m - 1) { /* convergence of bottom singular values */ --m; if (*idbg > 0) { s_wsle(&io__276); do_lio(&c__9, &c__1, "convergence", 11L); e_wsle(); } goto L999; } L1003: ++ll; /* e(ll) ... e(m-1) are nonzero */ if (*idbg > 0) { s_wsle(&io__277); do_lio(&c__9, &c__1, "work on block ll,m=", 19L); do_lio(&c__3, &c__1, (char *)&ll, (ftnlen)sizeof(integer)); do_lio(&c__3, &c__1, (char *)&m, (ftnlen)sizeof(integer)); e_wsle(); s_wsle(&io__278); do_lio(&c__9, &c__1, "smin=", 5L); do_lio(&c__5, &c__1, (char *)&smin, (ftnlen)sizeof(doublereal)); e_wsle(); s_wsle(&io__279); do_lio(&c__9, &c__1, "smax=", 5L); do_lio(&c__5, &c__1, (char *)&smax, (ftnlen)sizeof(doublereal)); e_wsle(); } /* 2 by 2 block - handle specially to guarantee convergence */ if (ll == m - 1) { /* after one step */ ++(*kount1); /* shift = sigmin(s(m-1),e(m-1),s(m)) */ /* rotate, setting e(m-1)=0 and s(m)=+-shift */ /* if (s(ll).eq.0.0d0) then */ /* f = 0.0d0 */ /* else */ /* f = (abs(s(ll)) - shift)*(dsign(1.0d0,s(ll))+shift/s(ll)) */ /* endif */ /* g=e(ll) */ /* call ndrotg(f,g,cs,sn) */ /* sg1=dsign(1.0d0,s(m)) */ /* sg2=dsign(1.0d0,cs) */ /* f = cs*s(ll) + sn*e(ll) */ /* g = sn*s(m) */ /* if (idbg.gt.0) then */ /* abss = cs*s(m) */ /* abse = -sn*s(ll) + cs*e(ll) */ /* endif */ /* if (wantv) call drot(p,v(1,ll),1,v(1,m),1,cs,sn) */ /* call ndrotg(f,g,cs,sn) */ /* s(ll)=f */ /* if (wantu.and.ll.lt.n) call drot(n,u(1,ll),1,u(1,m),1,cs,sn) */ /* e(ll) = 0.0d0 */ /* s(m) = shift * dsign(1.0d0,cs) * sg1 * sg2 */ /* if (idbg.gt.0) then */ /* print *,'2 by 2 block' */ /* print *,'shift=',shift */ /* print *,'check shift=',-sn*abse+cs*abss */ /* print *,'check zero=',cs*abse+sn*abss */ /* endif */ sig22_(&s[m - 1], &e[m - 1], &s[m], &sigmn, &sigmx, &sinr, &cosr, & sinl, &cosl); s[m - 1] = sigmx; e[m - 1] = (float)0.; s[m] = sigmn; if (wantv) { drot_(p, &v[ll * v_dim1 + 1], &c__1, &v[m * v_dim1 + 1], &c__1, & cosr, &sinr); } /* if wantu and ll.eq.n, then rotation trivial */ if (wantu && ll < *n) { drot_(n, &u[ll * u_dim1 + 1], &c__1, &u[m * u_dim1 + 1], &c__1, & cosl, &sinl); } goto L999; } /* choose shift direction if new submatrix */ /* if (ll.ne.oldll .or. m.ne.oldm) then */ /* choose shift direction if working on entirely new submatrix */ if (ll > oldm || m < oldll) { if ((d_1 = s[ll], abs(d_1)) >= (d_2 = s[m], abs(d_2)) && *iidir == 0 || *iidir == 1) { /* chase bulge from top (big end) to bottom (small end) */ /* if m=n+1, chase from top to bottom even if s(ll)=0 */ idir = 1; } else { /* chase bulge from bottom (big end) to top (small end) */ idir = 2; } } if (*idbg > 0) { s_wsle(&io__287); do_lio(&c__9, &c__1, "idir=", 5L); do_lio(&c__3, &c__1, (char *)&idir, (ftnlen)sizeof(integer)); e_wsle(); } /* compute lower bound on smallest singular value */ /* if old lower bound still good, do not recompute it */ /* if (ll.ne.oldll .or. m.ne.oldm .or. oldacc.ne.1) then */ /* compute lower bound */ /* sminl = smin */ /* oldacc = 1 */ /* if (sminl.gt.0.0d0) then */ /* if (idir.eq.1) then */ /* do 1004 lll=ll,m-1 */ /* abse = abs(e(lll)) */ /* abss = abs(s(lll)) */ /* if (abss.lt.abse) then */ /* sminl = sminl * (abss/abse) */ /* oldacc = -1 */ /* endif */ /* L1004: */ /* else */ /* do 1005 lll=ll,m-1 */ /* abse = abs(e(lll)) */ /* abss = abs(s(lll+1)) */ /* if (abss.lt.abse) then */ /* sminl = sminl * (abss/abse) */ /* oldacc = -1 */ /* endif */ /* L1005: */ /* endif */ /* endif */ /* oldll = ll */ /* oldm = m */ /* sminl is lower bound on smallest singular value */ /* within a factor of sqrt(m*(m+1)/2) */ /* if oldacc = 1 as well, sminl is also upper bound */ /* note that smin is always an upper bound */ /* compute convergence threshold */ /* thresh = tol*sminl */ /* endif */ /* if (idbg.gt.0) then */ /* print *,'oldll,oldm,oldacc=',oldll,oldm,oldacc */ /* print *,'sminl=',sminl */ /* print *,'thresh=',thresh */ /* endif */ /* test again for convergence using new thresh */ /* iconv = 0 */ /* do 1014 lll=ll,m-1 */ /* if (dabs(e(lll)).le.thresh) then */ /* e(lll) = 0.0d0 */ /* iconv = 1 */ /* endif */ /* L1014: */ /* if (iconv.eq.1) goto 999 */ /* Kahan's convergence test */ sminl = (float)0.; if (*tol > 0.) { if (idir == 1) { /* forward direction */ /* apply test on bottom 2 by 2 only */ if ((d_1 = e[m - 1], abs(d_1)) <= *tol * (d_2 = s[m], abs(d_2))) { /* convergence of bottom element */ e[m - 1] = 0.; goto L999; } /* apply test in forward direction */ mu = (d_1 = s[ll], abs(d_1)); sminl = mu; i_1 = m - 1; for (lll = ll; lll <= i_1; ++lll) { if ((d_1 = e[lll], abs(d_1)) <= *tol * mu) { /* test for negligibility satisfied */ if (*idbg >= 1) { s_wsle(&io__290); do_lio(&c__9, &c__1, "knew: e(lll),mu=", 16L); do_lio(&c__5, &c__1, (char *)&e[lll], (ftnlen)sizeof( doublereal)); do_lio(&c__5, &c__1, (char *)&mu, (ftnlen)sizeof( doublereal)); e_wsle(); } e[lll] = 0.; goto L999; } else { mu = (d_1 = s[lll + 1], abs(d_1)) * (mu / (mu + (d_2 = e[ lll], abs(d_2)))); } sminl = min(sminl,mu); /* L3330: */ } } else { /* idir=2, backwards direction */ /* apply test on top 2 by 2 only */ if ((d_1 = e[ll], abs(d_1)) <= *tol * (d_2 = s[ll], abs(d_2))) { /* convergence of top element */ e[ll] = 0.; goto L999; } /* apply test in backward direction */ lambda = (d_1 = s[m], abs(d_1)); sminl = lambda; i_1 = ll; for (lll = m - 1; lll >= i_1; --lll) { if ((d_1 = e[lll], abs(d_1)) <= *tol * lambda) { /* test for negligibility satisfied */ if (*idbg >= 1) { s_wsle(&io__292); do_lio(&c__9, &c__1, "knew: e(lll),lambda=", 20L); do_lio(&c__5, &c__1, (char *)&e[lll], (ftnlen)sizeof( doublereal)); do_lio(&c__5, &c__1, (char *)&lambda, (ftnlen)sizeof( doublereal)); e_wsle(); } e[lll] = 0.; goto L999; } else { lambda = (d_1 = s[lll], abs(d_1)) * (lambda / (lambda + ( d_2 = e[lll], abs(d_2)))); } sminl = min(sminl,lambda); /* L3331: */ } } thresh = *tol * sminl; /* thresh = 0 */ } oldll = ll; oldm = m; /* test for zero shift */ test = nct * *tol * (sminl / smax) + 1.; if (test == 1. && *ifull != 1 && *tol > 0. || *ifull == 2) { /* do a zero shift so that roundoff does not contaminate */ /* smallest singular value */ shift = 0.; if (*idbg > 0) { s_wsle(&io__295); do_lio(&c__9, &c__1, "sminl test for shift is zero", 28L); e_wsle(); } } else { /* compute shift from 2 by 2 block at end of matrix */ if (idir == 1) { smm1 = s[m - 1]; emm1 = e[m - 1]; sm = s[m]; sll = s[ll]; } else { smm1 = s[ll + 1]; emm1 = e[ll]; sm = s[ll]; sll = s[m]; } if (*idbg > 0) { s_wsle(&io__300); do_lio(&c__9, &c__1, "smm1,emm1,sm=", 13L); do_lio(&c__5, &c__1, (char *)&smm1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&emm1, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&sm, (ftnlen)sizeof(doublereal)); e_wsle(); } shift = sigmin_(&smm1, &emm1, &sm); if (*idbg > 0) { s_wsle(&io__301); do_lio(&c__9, &c__1, "sigma-min of 2 by 2 corner=", 27L); do_lio(&c__5, &c__1, (char *)&shift, (ftnlen)sizeof(doublereal)); e_wsle(); } if (*tol > 0.) { if (shift > nct * smin) { ++(*limshf); shift = nct * smin; if (*idbg > 0) { s_wsle(&io__302); do_lio(&c__9, &c__1, "shift limited", 13L); e_wsle(); } } if (*idbg > 0) { s_wsle(&io__303); do_lio(&c__9, &c__1, "shift=", 6L); do_lio(&c__5, &c__1, (char *)&shift, (ftnlen)sizeof( doublereal)); e_wsle(); } temp = shift / sll; if (*idbg > 0) { s_wsle(&io__305); do_lio(&c__9, &c__1, "temp=", 5L); do_lio(&c__5, &c__1, (char *)&temp, (ftnlen)sizeof(doublereal) ); e_wsle(); } /* Computing 2nd power */ d_1 = temp; test = 1. - d_1 * d_1; /* test to see if shift negligible */ if (*ifull != 1 && test == 1.) { shift = 0.; } } else { /* if shift much larger than s(ll), first rotation could be 0, */ /* leading to infinite loop; avoid by doing 0 shift in this case */ if (shift > (d_1 = s[ll], abs(d_1))) { test = s[ll] / shift; if (test + 1. == 1.) { ++(*limshf); if (*ifull != 1) { shift = (float)0.; } if (*idbg > 0 && *ifull != 1) { s_wsle(&io__306); do_lio(&c__9, &c__1, "shift limited", 13L); e_wsle(); } } } test = smax + smin; if (test == smax && *ifull != 1) { shift = (float)0.; } } if (*idbg > 0) { s_wsle(&io__307); do_lio(&c__9, &c__1, "test,shift=", 11L); do_lio(&c__5, &c__1, (char *)&test, (ftnlen)sizeof(doublereal)); do_lio(&c__5, &c__1, (char *)&shift, (ftnlen)sizeof(doublereal)); e_wsle(); } } /* increment iteration counter */ iter = iter + m - ll; ++(*kount); if (*idbg > 1) { s_wsle(&io__308); do_lio(&c__9, &c__1, "s,e before qr", 13L); e_wsle(); prse_(&ll, &m, n, p, &s[1], &e[1]); } /* if shift = 0, do simplified qr iteration */ if (shift == 0.) { *kount2 = *kount2 + m - ll; /* if idir=1, chase bulge from top to bottom */ if (idir == 1) { if (*idbg > 2) { s_wsle(&io__309); do_lio(&c__9, &c__1, "qr with zero shift, top to bottom", 33L) ; e_wsle(); } oldcs = 1.; f = s[ll]; g = e[ll]; i_1 = m - 1; for (k = ll; k <= i_1; ++k) { /* if (idbg.gt.2) print *,'qr inner loop, k=',k */ /* if (idbg.gt.3) print *,'f,g=',f,g */ ndrotg_(&f, &g, &cs, &sn); /* Computing MAX */ d_1 = *maxsin, d_2 = abs(sn); *maxsin = max(d_2,d_1); /* if (idbg.gt.3) print *,'f,cs,sn=',f,cs,sn */ if (wantv) { drot_(p, &v[k * v_dim1 + 1], &c__1, &v[(k + 1) * v_dim1 + 1], &c__1, &cs, &sn); } if (k != ll) { e[k - 1] = oldsn * f; } /* if (k.ne.ll .and. idbg.gt.3) print *,'e(k-1)=',e( k-1) */ f = oldcs * f; /* if (idbg.gt.3) print *,'f=',f */ temp = s[k + 1]; /* if (idbg.gt.3) print *,'temp=',temp */ g = temp * sn; /* if (idbg.gt.3) print *,'g=',g */ gg = temp * cs; /* if (idbg.gt.3) print *,'gg=',gg */ ndrotg_(&f, &g, &cs, &sn); /* Computing MAX */ d_1 = *maxsin, d_2 = abs(sn); *maxsin = max(d_2,d_1); /* if (idbg.gt.3) print *,'f,cs,sn=',f,cs,sn */ /* if wantu and k.eq.n, then s(k+1)=0 so g=0 so cs= 1 and sn=0 */ if (wantu && k < *n) { drot_(n, &u[k * u_dim1 + 1], &c__1, &u[(k + 1) * u_dim1 + 1], &c__1, &cs, &sn); } s[k] = f; /* if (idbg.gt.3) print *,'s(k)=',s(k) */ f = gg; /* if (idbg.gt.3) print *,'f=',f */ g = e[k + 1]; /* if (idbg.gt.3) print *,'g=',g */ oldcs = cs; /* if (idbg.gt.3) print *,'oldcs=',oldcs */ oldsn = sn; /* if (idbg.gt.3) print *,'oldsn=',oldsn */ /* if (idbg.gt.2) call prse(ll,m,n,p,s,e) */ /* L1006: */ } e[m - 1] = gg * sn; /* if (idbg.gt.3) print *,'e(m-1)=',e(m-1) */ s[m] = gg * cs; /* if (idbg.gt.3) print *,'s(m)=',s(m) */ /* test convergence */ if (*idbg > 0) { s_wsle(&io__318); do_lio(&c__9, &c__1, "convergence decision for zero shift to\ p to bottom", 49L); e_wsle(); s_wsle(&io__319); do_lio(&c__9, &c__1, "e(m-1), threshold=", 18L); do_lio(&c__5, &c__1, (char *)&e[m - 1], (ftnlen)sizeof( doublereal)); do_lio(&c__5, &c__1, (char *)&thresh, (ftnlen)sizeof( doublereal)); e_wsle(); if ((d_1 = e[m - 1], abs(d_1)) <= thresh) { s_wsle(&io__320); do_lio(&c__9, &c__1, "***converged***", 15L); e_wsle(); } } if ((d_1 = e[m - 1], abs(d_1)) <= thresh) { e[m - 1] = 0.; } } else { /* (idir=2, so chase bulge from bottom to top) */ if (*idbg > 2) { s_wsle(&io__321); do_lio(&c__9, &c__1, "qr with zero shift, bottom to top", 33L) ; e_wsle(); } oldcs = 1.; f = s[m]; g = e[m - 1]; i_1 = ll + 1; for (k = m; k >= i_1; --k) { /* if (idbg.gt.2) print *,'qr inner loop, k=',k */ /* if (idbg.gt.3) print *,'f,g=',f,g */ ndrotg_(&f, &g, &cs, &sn); /* Computing MAX */ d_1 = *maxsin, d_2 = abs(sn); *maxsin = max(d_2,d_1); /* if (idbg.gt.3) print *,'f,cs,sn=',f,cs,sn */ /* if m=n+1, always chase from top to bottom so no test for */ /* k.lt.n necessary */ if (wantu) { d_1 = -sn; drot_(n, &u[(k - 1) * u_dim1 + 1], &c__1, &u[k * u_dim1 + 1], &c__1, &cs, &d_1); } if (k != m) { e[k] = oldsn * f; } /* if (k.ne.m .and. idbg.gt.3) print *,'e(k)=',e(k) */ f = oldcs * f; /* if (idbg.gt.3) print *,'f=',f */ temp = s[k - 1]; /* if (idbg.gt.3) print *,'temp=',temp */ g = sn * temp; /* if (idbg.gt.3) print *,'g=',g */ gg = cs * temp; /* if (idbg.gt.3) print *,'gg=',gg */ ndrotg_(&f, &g, &cs, &sn); /* Computing MAX */ d_1 = *maxsin, d_2 = abs(sn); *maxsin = max(d_2,d_1); /* if (idbg.gt.3) print *,'f,cs,sn=',f,cs,sn */ if (wantv) { d_1 = -sn; drot_(p, &v[(k - 1) * v_dim1 + 1], &c__1, &v[k * v_dim1 + 1], &c__1, &cs, &d_1); } s[k] = f; /* if (idbg.gt.3) print *,'s(k)=',s(k) */ f = gg; /* if (idbg.gt.3) print *,'f=',f */ if (k != ll + 1) { g = e[k - 2]; } /* if (k.ne.ll+1 .and. idbg.gt.3) print *,'g=',g */ oldcs = cs; /* if (idbg.gt.3) print *,'oldcs=',oldcs */ oldsn = sn; /* if (idbg.gt.3) print *,'oldsn=',oldsn */ /* if (idbg.gt.2) call prse(ll,m,n,p,s,e) */ /* L1007: */ } e[ll] = gg * sn; /* if (idbg.gt.3) print *,'e(ll)=',e(ll) */ s[ll] = gg * cs; /* if (idbg.gt.3) print *,'s(ll)=',s(ll) */ /* test convergence */ if (*idbg > 0) { s_wsle(&io__322); do_lio(&c__9, &c__1, "convergence decision for zero shift bo\ ttom to top", 49L); e_wsle(); s_wsle(&io__323); do_lio(&c__9, &c__1, "e(ll), threshold=", 17L); do_lio(&c__5, &c__1, (char *)&e[ll], (ftnlen)sizeof( doublereal)); do_lio(&c__5, &c__1, (char *)&thresh, (ftnlen)sizeof( doublereal)); e_wsle(); if ((d_1 = e[ll], abs(d_1)) <= thresh) { s_wsle(&io__324); do_lio(&c__9, &c__1, "***converged***", 15L); e_wsle(); } } if ((d_1 = e[ll], abs(d_1)) <= thresh) { e[ll] = 0.; } } } else { /* (shift.ne.0, so do standard qr iteration) */ *kount1 = *kount1 + m - ll; /* if idir=1, chase bulge from top to bottom */ if (idir == 1) { if (*idbg > 2) { s_wsle(&io__325); do_lio(&c__9, &c__1, "qr with nonzero shift, top to bottom", 36L); e_wsle(); } f = ((d_1 = s[ll], abs(d_1)) - shift) * (d_sign(&c_b30, &s[ll]) + shift / s[ll]); g = e[ll]; i_1 = m - 1; for (k = ll; k <= i_1; ++k) { /* if (idbg.gt.2) print *,'qr inner loop, k=',k */ /* if (idbg.gt.3) print *,'f,g=',f,g */ ndrotg_(&f, &g, &cs, &sn); /* if (idbg.gt.3) print *,'f,cs,sn=',f,cs,sn */ if (k != ll) { e[k - 1] = f; } /* if (k.ne.ll .and. idbg.gt.3) print *,'e(k-1)=',e( k-1) */ f = cs * s[k] + sn * e[k]; /* if (idbg.gt.3) print *,'f=',f */ e[k] = cs * e[k] - sn * s[k]; /* if (idbg.gt.3) print *,'e(k)=',e(k) */ g = sn * s[k + 1]; /* if (idbg.gt.3) print *,'g=',g */ s[k + 1] = cs * s[k + 1]; /* if (idbg.gt.3) print *,'s(k+1)=',s(k+1) */ if (wantv) { drot_(p, &v[k * v_dim1 + 1], &c__1, &v[(k + 1) * v_dim1 + 1], &c__1, &cs, &sn); } ndrotg_(&f, &g, &cs, &sn); /* if (idbg.gt.3) print *,'f,cs,sn=',f,cs,sn */ s[k] = f; /* if (idbg.gt.3) print *,'s(k)=',s(k) */ f = cs * e[k] + sn * s[k + 1]; /* if (idbg.gt.3) print *,'f=',f */ s[k + 1] = -sn * e[k] + cs * s[k + 1]; /* if (idbg.gt.3) print *,'s(k+1)=',s(k+1) */ g = sn * e[k + 1]; /* if (idbg.gt.3) print *,'g=',g */ e[k + 1] = cs * e[k + 1]; /* if (idbg.gt.3) print *,'e(k+1)=',e(k+1) */ /* test for k.lt.n seems unnecessary since k=n causes zero */ /* shift, so test removed from original code */ if (wantu) { drot_(n, &u[k * u_dim1 + 1], &c__1, &u[(k + 1) * u_dim1 + 1], &c__1, &cs, &sn); } /* if (idbg.gt.2) call prse(ll,m,n,p,s,e) */ /* L1008: */ } e[m - 1] = f; /* if (idbg.gt.3) print *,'e(m-1)=',e(m-1) */ /* check convergence */ if (*idbg > 0) { s_wsle(&io__326); do_lio(&c__9, &c__1, "convergence decision for shift top to \ bottom", 44L); e_wsle(); s_wsle(&io__327); do_lio(&c__9, &c__1, "e(m-1), threshold=", 18L); do_lio(&c__5, &c__1, (char *)&e[m - 1], (ftnlen)sizeof( doublereal)); do_lio(&c__5, &c__1, (char *)&thresh, (ftnlen)sizeof( doublereal)); e_wsle(); if ((d_1 = e[m - 1], abs(d_1)) <= thresh) { s_wsle(&io__328); do_lio(&c__9, &c__1, "***converged***", 15L); e_wsle(); } } if ((d_1 = e[m - 1], abs(d_1)) <= thresh) { e[m - 1] = 0.; } } else { /* (idir=2, so chase bulge from bottom to top) */ if (*idbg > 2) { s_wsle(&io__329); do_lio(&c__9, &c__1, "qr with nonzero shift, bottom to top", 36L); e_wsle(); } f = ((d_1 = s[m], abs(d_1)) - shift) * (d_sign(&c_b30, &s[m]) + shift / s[m]); g = e[m - 1]; i_1 = ll + 1; for (k = m; k >= i_1; --k) { /* if (idbg.gt.2) print *,'qr inner loop, k=',k */ /* if (idbg.gt.3) print *,'f,g=',f,g */ ndrotg_(&f, &g, &cs, &sn); /* if (idbg.gt.3) print *,'f,cs,sn=',f,cs,sn */ if (k != m) { e[k] = f; } /* if (k.ne.m .and. idbg.gt.3) print *,'e(k)=',e(k) */ f = cs * s[k] + sn * e[k - 1]; /* if (idbg.gt.3) print *,'f=',f */ e[k - 1] = -sn * s[k] + cs * e[k - 1]; /* if (idbg.gt.3) print *,'e(k-1)=',e(k-1) */ g = sn * s[k - 1]; /* if (idbg.gt.3) print *,'g=',g */ s[k - 1] = cs * s[k - 1]; /* if (idbg.gt.3) print *,'s(k-1)=',s(k-1) */ if (wantu && k <= *n) { d_1 = -sn; drot_(n, &u[(k - 1) * u_dim1 + 1], &c__1, &u[k * u_dim1 + 1], &c__1, &cs, &d_1); } ndrotg_(&f, &g, &cs, &sn); /* if (idbg.gt.3) print *,'f,cs,sn=',f,cs,sn */ if (wantv) { d_1 = -sn; drot_(p, &v[(k - 1) * v_dim1 + 1], &c__1, &v[k * v_dim1 + 1], &c__1, &cs, &d_1); } s[k] = f; /* if (idbg.gt.3) print *,'s(k)=',s(k) */ f = sn * s[k - 1] + cs * e[k - 1]; /* if (idbg.gt.3) print *,'f=',f */ s[k - 1] = cs * s[k - 1] - sn * e[k - 1]; /* if (idbg.gt.3) print *,'s(k-1)=',s(k-1) */ if (k != ll + 1) { g = sn * e[k - 2]; /* if (idbg.gt.3) print *,'g=',g */ e[k - 2] = cs * e[k - 2]; /* if (idbg.gt.3) print *,'e(k-2)=',e(k-2) */ } /* if (idbg.gt.2) call prse(ll,m,n,p,s,e) */ /* L1009: */ } e[ll] = f; /* if (idbg.gt.3) print *,'e(ll)=',e(ll) */ /* test convergence */ if (*idbg > 0) { s_wsle(&io__330); do_lio(&c__9, &c__1, "convergence decision for shift bottom \ to top", 44L); e_wsle(); s_wsle(&io__331); do_lio(&c__9, &c__1, "e(ll), threshold=", 17L); do_lio(&c__5, &c__1, (char *)&e[ll], (ftnlen)sizeof( doublereal)); do_lio(&c__5, &c__1, (char *)&thresh, (ftnlen)sizeof( doublereal)); e_wsle(); if ((d_1 = e[ll], abs(d_1)) <= thresh) { s_wsle(&io__332); do_lio(&c__9, &c__1, "***converged***", 15L); e_wsle(); } } if ((d_1 = e[ll], abs(d_1)) <= thresh) { e[ll] = 0.; } } } if (*idbg > 1) { s_wsle(&io__333); do_lio(&c__9, &c__1, "s,e after qr", 12L); e_wsle(); prse_(&ll, &m, n, p, &s[1], &e[1]); } /* qr iteration finished, go back to check convergence */ goto L999; L998: /* make singular values positive */ m = min(*n,*p); i_1 = m; for (i = 1; i <= i_1; ++i) { if (s[i] < 0.) { s[i] = -s[i]; if (wantv) { dscal_(p, &c_b338, &v[i * v_dim1 + 1], &c__1); } } /* L1010: */ } /* sort singular values from largest at top to smallest */ /* at bottom (use insertion sort) */ i_1 = m - 1; for (i = 1; i <= i_1; ++i) { /* scan for smallest s(j) */ iisub = 1; smin = s[1]; i_2 = m + 1 - i; for (j = 2; j <= i_2; ++j) { if (s[j] < smin) { iisub = j; smin = s[j]; } /* L1012: */ } if (iisub != m + 1 - i) { /* swap singular values, vectors */ temp = s[m + 1 - i]; s[m + 1 - i] = s[iisub]; s[iisub] = temp; if (wantv) { dswap_(p, &v[(m + 1 - i) * v_dim1 + 1], &c__1, &v[iisub * v_dim1 + 1], &c__1); } if (wantu) { dswap_(n, &u[(m + 1 - i) * u_dim1 + 1], &c__1, &u[iisub * u_dim1 + 1], &c__1); } } /* L1011: */ } /* finished, return */ return 0; L997: /* maximum number of iterations exceeded */ for (i = m - 1; i >= 1; --i) { *info = i + 1; if (e[i] != 0.) { goto L996; } /* L1013: */ } L996: return 0; } /* ndsvd_ */ /* Subroutine */ int prse_(ll, m, nrow, ncol, s, e) integer *ll, *m; integer *nrow; integer *ncol; doublereal *s, *e; { /* Format strings */ static char fmt_100[] = "(22x,\002s(.)\002,22x,\002e(.) for ll,m=\002,2i\ 3)"; static char fmt_101[] = "(2d26.17)"; /* System generated locals */ integer i_1; /* Builtin functions */ integer s_wsfe(), do_fio(), e_wsfe(); /* Local variables */ static integer i; /* Fortran I/O blocks */ static cilist io__335 = { 0, 6, 0, fmt_100, 0 }; static cilist io__337 = { 0, 6, 0, fmt_101, 0 }; static cilist io__338 = { 0, 6, 0, fmt_101, 0 }; static cilist io__339 = { 0, 6, 0, fmt_101, 0 }; /* Parameter adjustments */ --s; --e; /* Function Body */ /* debug routine to print s,e */ s_wsfe(&io__335); do_fio(&c__1, (char *)&(*ll), (ftnlen)sizeof(integer)); do_fio(&c__1, (char *)&(*m), (ftnlen)sizeof(integer)); e_wsfe(); i_1 = *m - 1; for (i = *ll; i <= i_1; ++i) { s_wsfe(&io__337); do_fio(&c__1, (char *)&s[i], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&e[i], (ftnlen)sizeof(doublereal)); e_wsfe(); /* L1: */ } if (*m >= *ncol) { s_wsfe(&io__338); do_fio(&c__1, (char *)&s[*m], (ftnlen)sizeof(doublereal)); e_wsfe(); } if (*m < *ncol) { s_wsfe(&io__339); do_fio(&c__1, (char *)&s[*m], (ftnlen)sizeof(doublereal)); do_fio(&c__1, (char *)&e[*m], (ftnlen)sizeof(doublereal)); e_wsfe(); } return 0; } /* prse_ */ /* Subroutine */ int sig22_(a, b, c, sigmin, sigmax, snr, csr, snl, csl) doublereal *a, *b, *c, *sigmin, *sigmax, *snr, *csr, *snl, *csl; { /* System generated locals */ doublereal d_1, d_2; /* Builtin functions */ double d_sign(), sqrt(); /* Local variables */ static doublereal absa, absb, absc, acmn, acmx, sgna, sgnb, sgnc, cosl, sinl, cosr, temp, sinr, temp1, temp2, temp3, sgnmn, sgnmx, ac, ca; static integer ia; static doublereal absbac; static integer ib; static doublereal as, at, au; extern /* Subroutine */ int sndrtg_(); static doublereal bac; /* compute the singular value decomposition of the 2 by 2 */ /* upper triangular matrix [[a,b];[0,c]] */ /* inputs - */ /* a,b,c - real*8 - matrix entries */ /* outputs - */ /* sigmin - real*8 - +-smaller singular value */ /* sigmax - real*8 - +-larger singular value */ /* snr, csr - real*8 - sin and cos of right rotation (see below) */ /* snl, csl - real*8 - sin and cos of left rotation (see below) */ /* [ csl snl ] * [ a b ] * [ csr -snr ] = [ sigmax 0 ] */ /* [ -snl csl ] [ 0 c ] [ snr csr ] [ 0 sigmin ] */ /* barring over/underflow all output quantities are correct to */ /* within a few units in their last places */ /* let UF denote the underflow and OF the overflow threshold */ /* let eps denote the machine precision */ /* overflow is impossible unless the true value of sigmax exceeds */ /* OF (or does so within a few units in its last place) */ /* underflow cannot adversely effect sigmin,sigmax unless they are */ /* less than UF/eps for conventional underflow */ /* underflow cannot adversely effect sigmin,sigmax if underflow is */ /* gradual and results normalized */ /* overflow is impossible in computing sinr, cosr, sinl, cosl */ /* underflow can adversely effect results only if sigmax 0.) { /* 0 <= absbac <= 1 */ absbac = absb / acmx; /* 1 <= temp3 <= 1+sqrt(5), underflow harmless */ temp3 = absbac + sqrt(au + 4.); /* 1/3 <= temp3 <= (1+sqrt(10))/2 */ temp3 /= absbac * temp3 + 2.; sinr = d_sign(&c_b30, b); cosr = d_sign(&temp3, a); sinl = d_sign(&temp3, c); cosl = sinr; sgnmn = sgna * sgnb * sgnc; sgnmx = sgnb; } else { /* matrix diagonal */ sinr = 0.; cosr = 1.; sinl = 0.; cosl = 1.; sgnmn = sgnc; sgnmx = sgna; } } else { /* at .ne. 0, so eps/2 <= at <= 1 */ /* eps/2 <= temp3 <= 1 + sqrt(10) */ temp3 = au + temp1 * temp2; if (absa < absc) { /* abs(ac) <= 1 */ ac = *a / *c; /* eps <= sinr <= sqrt(13)+3 */ /* Computing 2nd power */ d_1 = as * at + au; sinr = sqrt(d_1 * d_1 + ac * 4. * ac * au) + as * at + au; /* abs(cosr) <= 2; if underflow, true cosr abs(*g)) { t = *g / *f; tt = sqrt(t * t + 1.); d_1 = 1. / tt; *cs = d_sign(&d_1, f); d_1 = t * *cs; *sn = d_sign(&d_1, g); } else { t = *f / *g; tt = sqrt(t * t + 1.); d_1 = 1. / tt; *sn = d_sign(&d_1, g); d_1 = t * *sn; *cs = d_sign(&d_1, f); } } return 0; } /* sndrtg_ */