Mercurial > hg > octave-nkf
changeset 3596:edcaebe1b81b
[project @ 2000-02-10 09:26:48 by jwe]
| author | jwe |
|---|---|
| date | Thu, 10 Feb 2000 09:26:51 +0000 |
| parents | fa5817b05b0f |
| children | 26662775f4e9 |
| files | libcruft/lapack/dbdsqr.f libcruft/lapack/dgelss.f libcruft/lapack/dgesvd.f libcruft/lapack/dlasq1.f libcruft/lapack/dlasq2.f libcruft/lapack/dlasq3.f libcruft/lapack/dlasq4.f libcruft/lapack/dlasq5.f libcruft/lapack/dlasq6.f libcruft/lapack/zbdsqr.f libcruft/lapack/zgelss.f libcruft/lapack/zgesvd.f libcruft/lapack/zhetd2.f |
| diffstat | 13 files changed, 387 insertions(+), 266 deletions(-) [+] |
line wrap: on
line diff
--- a/libcruft/lapack/dbdsqr.f +++ b/libcruft/lapack/dbdsqr.f @@ -4,7 +4,7 @@ * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. CHARACTER UPLO @@ -95,9 +95,7 @@ * The leading dimension of the array C. * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. * -* WORK (workspace) DOUBLE PRECISION array, dimension -* 2*N if only singular values wanted (NCVT = NRU = NCC = 0) -* max( 1, 4*N-4 ) otherwise +* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) * * INFO (output) INTEGER * = 0: successful exit
--- a/libcruft/lapack/dgelss.f +++ b/libcruft/lapack/dgelss.f @@ -4,7 +4,7 @@ * -- LAPACK driver routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK @@ -175,7 +175,7 @@ * * Compute workspace needed for DBDSQR * - BDSPAC = MAX( 1, 5*N-4 ) + BDSPAC = MAX( 1, 5*N ) MAXWRK = MAX( MAXWRK, 3*N+( MM+N )* $ ILAENV( 1, 'DGEBRD', ' ', MM, N, -1, -1 ) ) MAXWRK = MAX( MAXWRK, 3*N+NRHS* @@ -191,7 +191,7 @@ * * Compute workspace needed for DBDSQR * - BDSPAC = MAX( 1, 5*M-4 ) + BDSPAC = MAX( 1, 5*M ) MINWRK = MAX( 3*M+NRHS, 3*M+N, BDSPAC ) IF( N.GE.MNTHR ) THEN *
--- a/libcruft/lapack/dgesvd.f +++ b/libcruft/lapack/dgesvd.f @@ -4,7 +4,7 @@ * -- LAPACK driver routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. CHARACTER JOBU, JOBVT @@ -115,7 +115,7 @@ * * LWORK (input) INTEGER * The dimension of the array WORK. LWORK >= 1. -* LWORK >= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)-4). +* LWORK >= MAX(3*MIN(M,N)+MAX(M,N),5*MIN(M,N)). * For good performance, LWORK should generally be larger. * * If LWORK = -1, then a workspace query is assumed; the routine @@ -214,7 +214,7 @@ * * Compute space needed for DBDSQR * - BDSPAC = MAX( 3*N, 5*N-4 ) + BDSPAC = 5*N IF( M.GE.MNTHR ) THEN IF( WNTUN ) THEN * @@ -387,7 +387,7 @@ * * Compute space needed for DBDSQR * - BDSPAC = MAX( 3*M, 5*M-4 ) + BDSPAC = 5*M IF( N.GE.MNTHR ) THEN IF( WNTVN ) THEN *
--- a/libcruft/lapack/dlasq1.f +++ b/libcruft/lapack/dlasq1.f @@ -3,7 +3,7 @@ * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. INTEGER INFO, N @@ -26,13 +26,7 @@ * 1994, * * and the present implementation is described in "An implementation of -* dqds", LAPACK technical report. -* -* Note : DLASQ1 works only on machines which follow ieee-754 -* floating-point standard in their handling of infinities and NaNs. -* Normal execution of DLASQ1 may create NaNs and infinities and hence -* may abort due to a floating point exception in environments which -* do not conform to the ieee standard. +* the dqds Algorithm (Positive Case)", LAPACK Working Note. * * Arguments * ========= @@ -50,17 +44,17 @@ * of the bidiagonal matrix whose SVD is desired. * On exit, E is overwritten. * -* WORK (workspace) DOUBLE PRECISION array, dimension (2*N) +* WORK (workspace) DOUBLE PRECISION array, dimension (4*N) * * INFO (output) INTEGER -* = 0: successful exit -* < 0: if INFO = -i, the i-th argument had an illegal value +* = 0: successful exit +* < 0: if INFO = -i, the i-th argument had an illegal value * > 0: the algorithm failed -* = 1, a split was marked by a positive value in E -* = 2, current block of Z not diagonalized after 30*N -* iterations (in inner while loop) -* = 3, termination criterion of outer while loop not met -* (program created more than N unreduced blocks) +* = 1, a split was marked by a positive value in E +* = 2, current block of Z not diagonalized after 30*N +* iterations (in inner while loop) +* = 3, termination criterion of outer while loop not met +* (program created more than N unreduced blocks) * * ===================================================================== * @@ -70,10 +64,10 @@ * .. * .. Local Scalars .. INTEGER I, IINFO - DOUBLE PRECISION EPS, SCALE, SFMIN, SIGMN, SIGMX + DOUBLE PRECISION EPS, SCALE, SAFMIN, SIGMN, SIGMX * .. * .. External Subroutines .. - EXTERNAL DCOPY, DLAS2, DLASCL, DLASQ2, DLASRT, XERBLA + EXTERNAL DLAS2, DLASQ2, DLASRT, XERBLA * .. * .. External Functions .. DOUBLE PRECISION DLAMCH @@ -114,7 +108,7 @@ * IF( SIGMX.EQ.ZERO ) THEN CALL DLASRT( 'D', N, D, IINFO ) - GO TO 50 + RETURN END IF * DO 20 I = 1, N @@ -125,13 +119,13 @@ * input data makes scaling by a power of the radix pointless). * EPS = DLAMCH( 'Precision' ) - SFMIN = DLAMCH( 'Safe minimum' ) - SCALE = SQRT( EPS / SFMIN ) + SAFMIN = DLAMCH( 'Safe minimum' ) + SCALE = SQRT( EPS / SAFMIN ) CALL DCOPY( N, D, 1, WORK( 1 ), 2 ) CALL DCOPY( N-1, E, 1, WORK( 2 ), 2 ) CALL DLASCL( 'G', 0, 0, SIGMX, SCALE, 2*N-1, 1, WORK, 2*N-1, $ IINFO ) -* +* * Compute the q's and e's. * DO 30 I = 1, 2*N - 1 @@ -148,7 +142,6 @@ CALL DLASCL( 'G', 0, 0, SCALE, SIGMX, N, 1, D, N, IINFO ) END IF * - 50 CONTINUE RETURN * * End of DLASQ1
--- a/libcruft/lapack/dlasq2.f +++ b/libcruft/lapack/dlasq2.f @@ -1,9 +1,9 @@ SUBROUTINE DLASQ2( N, Z, INFO ) * -* -- LAPACK auxiliary routine (version 3.0) -- +* -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. INTEGER INFO, N @@ -15,7 +15,7 @@ * Purpose * ======= * -* DLASQ2 computes all the eigenvalues of the symmetric positive +* DLASQ2 computes all the eigenvalues of the symmetric positive * definite tridiagonal matrix associated with the qd array Z to high * relative accuracy are computed to high relative accuracy, in the * absence of denormalization, underflow and overflow. @@ -26,11 +26,10 @@ * Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the * symmetric tridiagonal to which it is similar. * -* Note : DLASQ2 works only on machines which follow ieee-754 -* floating-point standard in their handling of infinities and NaNs. -* Normal execution of DLASQ2 may create NaNs and infinities and hence -* may abort due to a floating point exception in environments which -* do not conform to the ieee standard. +* Note : DLASQ2 defines a logical variable, IEEE, which is true +* on machines which follow ieee-754 floating-point standard in their +* handling of infinities and NaNs, and false otherwise. This variable +* is passed to DLASQ3. * * Arguments * ========= @@ -41,9 +40,10 @@ * Z (workspace) DOUBLE PRECISION array, dimension ( 4*N ) * On entry Z holds the qd array. On exit, entries 1 to N hold * the eigenvalues in decreasing order, Z( 2*N+1 ) holds the -* trace, Z( 2*N+2 ) holds the sum of the eigenvalues, Z( 2*N+3 ) -* holds the iteration count, Z( 2*N+4 ) holds NDIVS/NIN^2, and -* Z( 2*N+5 ) holds the percentage of shifts that failed. +* trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If +* N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) +* holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of +* shifts that failed. * * INFO (output) INTEGER * = 0: successful exit @@ -55,7 +55,7 @@ * = 1, a split was marked by a positive value in E * = 2, current block of Z not diagonalized after 30*N * iterations (in inner while loop) -* = 3, termination criterion of outer while loop not met +* = 3, termination criterion of outer while loop not met * (program created more than N unreduced blocks) * * Further Details @@ -69,36 +69,39 @@ * .. Parameters .. DOUBLE PRECISION CBIAS PARAMETER ( CBIAS = 1.50D0 ) - DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, TEN, HNDRD + DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, - $ TWO = 2.0D0, FOUR = 4.0D0, TEN = 10.0D0, - $ HNDRD = 100.0D0 ) + $ TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 ) * .. * .. Local Scalars .. - INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K, + LOGICAL IEEE + INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K, $ N0, NBIG, NDIV, NFAIL, PP, SPLT - DOUBLE PRECISION D, DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, E, - $ EMAX, EMIN, EPS, EPS2, OLDEMN, QMAX, QMIN, S, - $ SIGMA, T, TAU, TEMP, TRACE, ZMAX + DOUBLE PRECISION D, DESIG, DMIN, E, EMAX, EMIN, EPS, OLDEMN, + $ QMAX, QMIN, S, SAFMIN, SIGMA, T, TEMP, TOL, + $ TOL2, TRACE, ZMAX * .. * .. External Subroutines .. - EXTERNAL DLASQ3, DLASQ5, DLASRT, XERBLA + EXTERNAL DLASQ3, DLASRT, XERBLA * .. * .. External Functions .. + INTEGER ILAENV DOUBLE PRECISION DLAMCH - EXTERNAL DLAMCH + EXTERNAL DLAMCH, ILAENV * .. * .. Intrinsic Functions .. - INTRINSIC ABS, DBLE, MAX, MIN, SQRT + INTRINSIC DBLE, MAX, MIN, SQRT * .. * .. Executable Statements .. -* +* * Test the input arguments. * (in case DLASQ2 is not called by DLASQ1) * INFO = 0 - EPS = DLAMCH( 'Precision' )*TEN - EPS2 = EPS**2 + EPS = DLAMCH( 'Precision' ) + SAFMIN = DLAMCH( 'Safe minimum' ) + TOL = EPS*HUNDRD + TOL2 = TOL**2 * IF( N.LT.0 ) THEN INFO = -1 @@ -129,8 +132,8 @@ Z( 1 ) = D END IF Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 ) - IF( Z( 2 ).GT.Z( 3 )*EPS2 ) THEN - T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) + IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN + T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) S = Z( 3 )*( Z( 2 ) / T ) IF( S.LE.T ) THEN S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) ) @@ -143,9 +146,6 @@ END IF Z( 2 ) = Z( 3 ) Z( 6 ) = Z( 2 ) + Z( 1 ) - Z( 7 ) = ZERO - Z( 8 ) = ZERO - Z( 9 ) = ZERO RETURN END IF * @@ -154,40 +154,47 @@ Z( 2*N ) = ZERO EMIN = Z( 2 ) QMAX = ZERO + ZMAX = ZERO D = ZERO E = ZERO * - DO 10 K = 1, N + DO 10 K = 1, 2*( N-1 ), 2 IF( Z( K ).LT.ZERO ) THEN INFO = -( 200+K ) CALL XERBLA( 'DLASQ2', 2 ) RETURN - ELSE IF( Z( N+K ).LT.ZERO ) THEN - INFO = -( 200+N+K ) + ELSE IF( Z( K+1 ).LT.ZERO ) THEN + INFO = -( 200+K+1 ) CALL XERBLA( 'DLASQ2', 2 ) RETURN END IF D = D + Z( K ) - E = E + Z( N+K ) + E = E + Z( K+1 ) QMAX = MAX( QMAX, Z( K ) ) + EMIN = MIN( EMIN, Z( K+1 ) ) + ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) ) 10 CONTINUE - ZMAX = QMAX - DO 20 K = 1, N - 1 - EMIN = MIN( EMIN, Z( N+K ) ) - ZMAX = MAX( ZMAX, Z( N+K ) ) - 20 CONTINUE + IF( Z( 2*N-1 ).LT.ZERO ) THEN + INFO = -( 200+2*N-1 ) + CALL XERBLA( 'DLASQ2', 2 ) + RETURN + END IF + D = D + Z( 2*N-1 ) + QMAX = MAX( QMAX, Z( 2*N-1 ) ) + ZMAX = MAX( QMAX, ZMAX ) * * Check for diagonality. * IF( E.EQ.ZERO ) THEN + DO 20 K = 2, N + Z( K ) = Z( 2*K-1 ) + 20 CONTINUE CALL DLASRT( 'D', N, Z, IINFO ) Z( 2*N-1 ) = D RETURN END IF * TRACE = D + E - I0 = 1 - N0 = N * * Check for zero data. * @@ -195,16 +202,24 @@ Z( 2*N-1 ) = ZERO RETURN END IF -* +* +* Check whether the machine is IEEE conformable. +* + IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND. + $ ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 +* * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). * DO 30 K = 2*N, 2, -2 - Z( 2*K ) = ZERO - Z( 2*K-1 ) = Z( K ) - Z( 2*K-2 ) = ZERO - Z( 2*K-3 ) = Z( K-1 ) + Z( 2*K ) = ZERO + Z( 2*K-1 ) = Z( K ) + Z( 2*K-2 ) = ZERO + Z( 2*K-3 ) = Z( K-1 ) 30 CONTINUE * + I0 = 1 + N0 = N +* * Reverse the qd-array, if warranted. * IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN @@ -225,47 +240,39 @@ * DO 80 K = 1, 2 * - IF( EMIN.LE.EPS2*QMAX ) THEN -* -* Li's reverse test. + D = Z( 4*N0+PP-3 ) + DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4 + IF( Z( I4-1 ).LE.TOL2*D ) THEN + Z( I4-1 ) = -ZERO + D = Z( I4-3 ) + ELSE + D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) ) + END IF + 50 CONTINUE * - D = Z( 4*N0+PP-3 ) - DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4 - IF( Z( I4-1 ).LE.EPS2*D ) THEN - Z( I4-1 ) = -ZERO - D = Z( I4-3 ) - ELSE - D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) ) - END IF - 50 CONTINUE -* -* dqd maps Z to ZZ plus Li's test. +* dqd maps Z to ZZ plus Li's test. * - EMIN = Z( 4*I0+PP+1 ) - D = Z( 4*I0+PP-3 ) - DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4 - IF( Z( I4-1 ).LE.EPS2*D ) THEN - Z( I4-1 ) = -ZERO - Z( I4-2*PP-2 ) = D - Z( I4-2*PP ) = ZERO - D = Z( I4+1 ) - EMIN = ZERO - ELSE - Z( I4-2*PP-2 ) = D + Z( I4-1 ) - Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / - $ Z( I4-2*PP-2 ) ) - D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) - EMIN = MIN( EMIN, Z( I4-2*PP ) ) - END IF - 60 CONTINUE - Z( 4*N0-PP-2 ) = D - ELSE - TAU = ZERO - CALL DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN, - $ DN1, DN2 ) -* - EMIN = Z( 4*N0 ) - END IF + EMIN = Z( 4*I0+PP+1 ) + D = Z( 4*I0+PP-3 ) + DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4 + Z( I4-2*PP-2 ) = D + Z( I4-1 ) + IF( Z( I4-1 ).LE.TOL2*D ) THEN + Z( I4-1 ) = -ZERO + Z( I4-2*PP-2 ) = D + Z( I4-2*PP ) = ZERO + D = Z( I4+1 ) + ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND. + $ SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN + TEMP = Z( I4+1 ) / Z( I4-2*PP-2 ) + Z( I4-2*PP ) = Z( I4-1 )*TEMP + D = D*TEMP + ELSE + Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) ) + D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) + END IF + EMIN = MIN( EMIN, Z( I4-2*PP ) ) + 60 CONTINUE + Z( 4*N0-PP-2 ) = D * * Now find qmax. * @@ -284,14 +291,14 @@ NDIV = 2*( N0-I0 ) * DO 140 IWHILA = 1, N + 1 - IF( N0.LT.1 ) + IF( N0.LT.1 ) $ GO TO 150 * -* While array unfinished do +* While array unfinished do * * E(N0) holds the value of SIGMA when submatrix in I0:N0 * splits from the rest of the array, but is negated. -* +* DESIG = ZERO IF( N0.EQ.N ) THEN SIGMA = ZERO @@ -306,8 +313,12 @@ * Find last unreduced submatrix's top index I0, find QMAX and * EMIN. Find Gershgorin-type bound if Q's much greater than E's. * - EMAX = ZERO - EMIN = ABS( Z( 4*N0-5 ) ) + EMAX = ZERO + IF( N0.GT.I0 ) THEN + EMIN = ABS( Z( 4*N0-5 ) ) + ELSE + EMIN = ZERO + END IF QMIN = Z( 4*N0-3 ) QMAX = QMIN DO 90 I4 = 4*N0, 8, -4 @@ -320,7 +331,7 @@ QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) EMIN = MIN( EMIN, Z( I4-5 ) ) 90 CONTINUE - I4 = 4 + I4 = 4 * 100 CONTINUE I0 = I4 / 4 @@ -335,32 +346,32 @@ * * Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong. * - PP = 0 + PP = 0 * NBIG = 30*( N0-I0+1 ) DO 120 IWHILB = 1, NBIG - IF( I0.GT.N0 ) + IF( I0.GT.N0 ) $ GO TO 130 * * While submatrix unfinished take a good dqds step. * CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, - $ ITER, NDIV ) + $ ITER, NDIV, IEEE ) * PP = 1 - PP * * When EMIN is very small check for splits. * IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN - IF( Z( 4*N0 ).LE.EPS2*QMAX .OR. Z( 4*N0-1 ).LE.EPS2* - $ SIGMA ) THEN + IF( Z( 4*N0 ).LE.TOL2*QMAX .OR. + $ Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN SPLT = I0 - 1 QMAX = Z( 4*I0-3 ) EMIN = Z( 4*I0-1 ) OLDEMN = Z( 4*I0 ) DO 110 I4 = 4*I0, 4*( N0-3 ), 4 - IF( Z( I4 ).LE.EPS2*Z( I4-3 ) .OR. Z( I4-1 ).LE. - $ EPS2*SIGMA ) THEN + IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR. + $ Z( I4-1 ).LE.TOL2*SIGMA ) THEN Z( I4-1 ) = -SIGMA SPLT = I4 / 4 QMAX = ZERO @@ -392,16 +403,16 @@ INFO = 3 RETURN * -* end IWHILA +* end IWHILA * 150 CONTINUE -* +* * Move q's to the front. -* +* DO 160 K = 2, N Z( K ) = Z( 4*K-3 ) 160 CONTINUE -* +* * Sort and compute sum of eigenvalues. * CALL DLASRT( 'D', N, Z, IINFO ) @@ -413,11 +424,11 @@ * * Store trace, sum(eigenvalues) and information on performance. * - Z( 2*N+1 ) = TRACE + Z( 2*N+1 ) = TRACE Z( 2*N+2 ) = E Z( 2*N+3 ) = DBLE( ITER ) Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 ) - Z( 2*N+5 ) = HNDRD*NFAIL / DBLE( ITER ) + Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER ) RETURN * * End of DLASQ2
--- a/libcruft/lapack/dlasq3.f +++ b/libcruft/lapack/dlasq3.f @@ -4,7 +4,7 @@ * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* December 22, 1999 +* October 31, 1999 * * .. Scalar Arguments .. LOGICAL IEEE @@ -223,7 +223,7 @@ * Z( 4*N0+PP-1 ) = ZERO GO TO 70 - ELSE IF( Z( 4*N0-PP ).LE.ZERO ) THEN + ELSE IF( Z( 4*N0+PP ).LE.ZERO ) THEN * * Possible unnecessary underflow in the e's. * Call safe dqd. @@ -231,7 +231,7 @@ Z( 4*N0+PP-1 ) = ZERO DMIN = ZERO GO TO 70 - ELSE IF( DMIN.EQ.ZERO .AND. TAU.EQ.ZERO ) THEN + ELSE IF( DMIN.EQ.ZERO .AND. DMIN1.LE.ZERO ) THEN * * Possible unnecessary underflow in the d's. * Call safe dqd.
--- a/libcruft/lapack/dlasq4.f +++ b/libcruft/lapack/dlasq4.f @@ -4,7 +4,7 @@ * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. INTEGER I0, N0, N0IN, PP, TTYPE @@ -16,7 +16,8 @@ * * Purpose * ======= -* DLASQ4 computes an approximation TAU to the smallest eigenvalue +* +* DLASQ4 computes an approximation TAU to the smallest eigenvalue * using values of d from the previous transform. * * I0 (input) INTEGER @@ -68,10 +69,10 @@ DOUBLE PRECISION CNST1, CNST2, CNST3 PARAMETER ( CNST1 = 0.5630D0, CNST2 = 1.010D0, $ CNST3 = 1.050D0 ) - DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HNDRD + DOUBLE PRECISION QURTR, THIRD, HALF, ZERO, ONE, TWO, HUNDRD PARAMETER ( QURTR = 0.250D0, THIRD = 0.3330D0, $ HALF = 0.50D0, ZERO = 0.0D0, ONE = 1.0D0, - $ TWO = 2.0D0, HNDRD = 100.0D0 ) + $ TWO = 2.0D0, HUNDRD = 100.0D0 ) * .. * .. Local Scalars .. INTEGER I4, NN, NP @@ -83,7 +84,7 @@ * .. Save statement .. SAVE G * .. -* .. Data statements .. +* .. Data statement .. DATA G / ZERO / * .. * .. Executable Statements .. @@ -96,7 +97,7 @@ TTYPE = -1 RETURN END IF -* +* NN = 4*N0 + PP IF( N0IN.EQ.N0 ) THEN * @@ -133,30 +134,40 @@ * * Case 4. * + TTYPE = -4 + S = QURTR*DMIN IF( DMIN.EQ.DN ) THEN GAM = DN A2 = ZERO + IF( Z( NN-5 ) .GT. Z( NN-7 ) ) + $ RETURN B2 = Z( NN-5 ) / Z( NN-7 ) NP = NN - 9 ELSE NP = NN - 2*PP B2 = Z( NP-2 ) GAM = DN1 + IF( Z( NP-4 ) .GT. Z( NP-2 ) ) + $ RETURN A2 = Z( NP-4 ) / Z( NP-2 ) + IF( Z( NN-9 ) .GT. Z( NN-11 ) ) + $ RETURN B2 = Z( NN-9 ) / Z( NN-11 ) NP = NN - 13 END IF * * Approximate contribution to norm squared from I < NN-1. * - IF( B2.EQ.ZERO ) - $ GO TO 20 A2 = A2 + B2 DO 10 I4 = NP, 4*I0 - 1 + PP, -4 + IF( B2.EQ.ZERO ) + $ GO TO 20 B1 = B2 + IF( Z( I4 ) .GT. Z( I4-2 ) ) + $ RETURN B2 = B2*( Z( I4 ) / Z( I4-2 ) ) A2 = A2 + B2 - IF( HNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) + IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) $ GO TO 20 10 CONTINUE 20 CONTINUE @@ -164,49 +175,48 @@ * * Rayleigh quotient residual bound. * - IF( A2.LT.CNST1 ) THEN - S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) - ELSE - S = QURTR*GAM - END IF - TTYPE = -4 + IF( A2.LT.CNST1 ) + $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) END IF ELSE IF( DMIN.EQ.DN2 ) THEN * * Case 5. * + TTYPE = -5 + S = QURTR*DMIN +* * Compute contribution to norm squared from I > NN-2. * NP = NN - 2*PP B1 = Z( NP-2 ) B2 = Z( NP-6 ) GAM = DN2 + IF( Z( NP-8 ).GT.B2 .OR. Z( NP-4 ).GT.B1 ) + $ RETURN A2 = ( Z( NP-8 ) / B2 )*( ONE+Z( NP-4 ) / B1 ) * * Approximate contribution to norm squared from I < NN-2. * IF( N0-I0.GT.2 ) THEN B2 = Z( NN-13 ) / Z( NN-15 ) - IF( B2.EQ.ZERO ) - $ GO TO 40 A2 = A2 + B2 DO 30 I4 = NN - 17, 4*I0 - 1 + PP, -4 + IF( B2.EQ.ZERO ) + $ GO TO 40 B1 = B2 + IF( Z( I4 ) .GT. Z( I4-2 ) ) + $ RETURN B2 = B2*( Z( I4 ) / Z( I4-2 ) ) A2 = A2 + B2 - IF( HNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) + IF( HUNDRD*MAX( B2, B1 ).LT.A2 .OR. CNST1.LT.A2 ) $ GO TO 40 30 CONTINUE 40 CONTINUE A2 = CNST3*A2 END IF * - IF( A2.LT.CNST1 ) THEN - S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) - ELSE - S = QURTR*GAM / ( ONE+A2 ) - END IF - TTYPE = -5 + IF( A2.LT.CNST1 ) + $ S = GAM*( ONE-SQRT( A2 ) ) / ( ONE+A2 ) ELSE * * Case 6, no information to guide us. @@ -226,19 +236,25 @@ * * One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. * - IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN + IF( DMIN1.EQ.DN1 .AND. DMIN2.EQ.DN2 ) THEN * * Cases 7 and 8. * + TTYPE = -7 + S = THIRD*DMIN1 + IF( Z( NN-5 ).GT.Z( NN-7 ) ) + $ RETURN B1 = Z( NN-5 ) / Z( NN-7 ) B2 = B1 IF( B2.EQ.ZERO ) $ GO TO 60 DO 50 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 A2 = B1 + IF( Z( I4 ).GT.Z( I4-2 ) ) + $ RETURN B1 = B1*( Z( I4 ) / Z( I4-2 ) ) B2 = B2 + B1 - IF( HNDRD*MAX( B1, A2 ).LT.B2 ) + IF( HUNDRD*MAX( B1, A2 ).LT.B2 ) $ GO TO 60 50 CONTINUE 60 CONTINUE @@ -246,11 +262,9 @@ A2 = DMIN1 / ( ONE+B2**2 ) GAP2 = HALF*DMIN2 - A2 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN - S = MAX( A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ), - $ THIRD*DMIN1 ) - TTYPE = -7 - ELSE - S = MAX( A2*( ONE-CNST2*B2 ), THIRD*DMIN1 ) + S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) + ELSE + S = MAX( S, A2*( ONE-CNST2*B2 ) ) TTYPE = -8 END IF ELSE @@ -269,15 +283,21 @@ * * Cases 10 and 11. * - IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN + IF( DMIN2.EQ.DN2 .AND. TWO*Z( NN-5 ).LT.Z( NN-7 ) ) THEN + TTYPE = -10 + S = THIRD*DMIN2 + IF( Z( NN-5 ).GT.Z( NN-7 ) ) + $ RETURN B1 = Z( NN-5 ) / Z( NN-7 ) B2 = B1 IF( B2.EQ.ZERO ) $ GO TO 80 DO 70 I4 = 4*N0 - 9 + PP, 4*I0 - 1 + PP, -4 + IF( Z( I4 ).GT.Z( I4-2 ) ) + $ RETURN B1 = B1*( Z( I4 ) / Z( I4-2 ) ) B2 = B2 + B1 - IF( HNDRD*B1.LT.B2 ) + IF( HUNDRD*B1.LT.B2 ) $ GO TO 80 70 CONTINUE 80 CONTINUE @@ -286,12 +306,10 @@ GAP2 = Z( NN-7 ) + Z( NN-9 ) - $ SQRT( Z( NN-11 ) )*SQRT( Z( NN-9 ) ) - A2 IF( GAP2.GT.ZERO .AND. GAP2.GT.B2*A2 ) THEN - S = MAX( A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ), - $ THIRD*DMIN2 ) - ELSE - S = MAX( A2*( ONE-CNST2*B2 ), THIRD*DMIN2 ) + S = MAX( S, A2*( ONE-CNST2*A2*( B2 / GAP2 )*B2 ) ) + ELSE + S = MAX( S, A2*( ONE-CNST2*B2 ) ) END IF - TTYPE = -10 ELSE S = QURTR*DMIN2 TTYPE = -11 @@ -300,7 +318,7 @@ * * Case 12, more than two eigenvalues deflated. No information. * - S = ZERO + S = ZERO TTYPE = -12 END IF *
--- a/libcruft/lapack/dlasq5.f +++ b/libcruft/lapack/dlasq5.f @@ -1,12 +1,13 @@ SUBROUTINE DLASQ5( I0, N0, Z, PP, TAU, DMIN, DMIN1, DMIN2, DN, - $ DNM1, DNM2 ) + $ DNM1, DNM2, IEEE ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. + LOGICAL IEEE INTEGER I0, N0, PP DOUBLE PRECISION DMIN, DMIN1, DMIN2, DN, DNM1, DNM2, TAU * .. @@ -16,7 +17,9 @@ * * Purpose * ======= -* DLASQ5 computes one dqds transform in ping-pong form. +* +* DLASQ5 computes one dqds transform in ping-pong form, one +* version for IEEE machines another for non IEEE machines. * * Arguments * ========= @@ -55,8 +58,15 @@ * DNM2 (output) DOUBLE PRECISION * d(N0-2). * +* IEEE (input) LOGICAL +* Flag for IEEE or non IEEE arithmetic. +* * ===================================================================== * +* .. Parameter .. + DOUBLE PRECISION ZERO + PARAMETER ( ZERO = 0.0D0 ) +* .. * .. Local Scalars .. INTEGER J4, J4P2 DOUBLE PRECISION D, EMIN, TEMP @@ -70,49 +80,113 @@ $ RETURN * J4 = 4*I0 + PP - 3 - EMIN = Z( J4+4 ) + EMIN = Z( J4+4 ) D = Z( J4 ) - TAU DMIN = D + DMIN1 = -Z( J4 ) * - IF( PP.EQ.0 ) THEN - DO 10 J4 = 4*I0, 4*( N0-3 ), 4 - Z( J4-2 ) = D + Z( J4-1 ) - TEMP = Z( J4+1 ) / Z( J4-2 ) - D = D*TEMP - TAU - DMIN = MIN( DMIN, D ) - Z( J4 ) = Z( J4-1 )*TEMP - EMIN = MIN( Z( J4 ), EMIN ) - 10 CONTINUE + IF( IEEE ) THEN +* +* Code for IEEE arithmetic. +* + IF( PP.EQ.0 ) THEN + DO 10 J4 = 4*I0, 4*( N0-3 ), 4 + Z( J4-2 ) = D + Z( J4-1 ) + TEMP = Z( J4+1 ) / Z( J4-2 ) + D = D*TEMP - TAU + DMIN = MIN( DMIN, D ) + Z( J4 ) = Z( J4-1 )*TEMP + EMIN = MIN( Z( J4 ), EMIN ) + 10 CONTINUE + ELSE + DO 20 J4 = 4*I0, 4*( N0-3 ), 4 + Z( J4-3 ) = D + Z( J4 ) + TEMP = Z( J4+2 ) / Z( J4-3 ) + D = D*TEMP - TAU + DMIN = MIN( DMIN, D ) + Z( J4-1 ) = Z( J4 )*TEMP + EMIN = MIN( Z( J4-1 ), EMIN ) + 20 CONTINUE + END IF +* +* Unroll last two steps. +* + DNM2 = D + DMIN2 = DMIN + J4 = 4*( N0-2 ) - PP + J4P2 = J4 + 2*PP - 1 + Z( J4-2 ) = DNM2 + Z( J4P2 ) + Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) + DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU + DMIN = MIN( DMIN, DNM1 ) +* + DMIN1 = DMIN + J4 = J4 + 4 + J4P2 = J4 + 2*PP - 1 + Z( J4-2 ) = DNM1 + Z( J4P2 ) + Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) + DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU + DMIN = MIN( DMIN, DN ) +* ELSE - DO 20 J4 = 4*I0, 4*( N0-3 ), 4 - Z( J4-3 ) = D + Z( J4 ) - TEMP = Z( J4+2 ) / Z( J4-3 ) - D = D*TEMP - TAU - DMIN = MIN( DMIN, D ) - Z( J4-1 ) = Z( J4 )*TEMP - EMIN = MIN( Z( J4-1 ), EMIN ) - 20 CONTINUE +* +* Code for non IEEE arithmetic. +* + IF( PP.EQ.0 ) THEN + DO 30 J4 = 4*I0, 4*( N0-3 ), 4 + Z( J4-2 ) = D + Z( J4-1 ) + IF( D.LE.ZERO ) THEN + RETURN + ELSE + Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) + D = Z( J4+1 )*( D / Z( J4-2 ) ) - TAU + END IF + DMIN = MIN( DMIN, D ) + EMIN = MIN( EMIN, Z( J4 ) ) + 30 CONTINUE + ELSE + DO 40 J4 = 4*I0, 4*( N0-3 ), 4 + Z( J4-3 ) = D + Z( J4 ) + IF( D.LE.ZERO ) THEN + RETURN + ELSE + Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) + D = Z( J4+2 )*( D / Z( J4-3 ) ) - TAU + END IF + DMIN = MIN( DMIN, D ) + EMIN = MIN( EMIN, Z( J4-1 ) ) + 40 CONTINUE + END IF +* +* Unroll last two steps. +* + DNM2 = D + DMIN2 = DMIN + J4 = 4*( N0-2 ) - PP + J4P2 = J4 + 2*PP - 1 + Z( J4-2 ) = DNM2 + Z( J4P2 ) + IF( DNM2.LE.ZERO ) THEN + RETURN + ELSE + Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) + DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU + END IF + DMIN = MIN( DMIN, DNM1 ) +* + DMIN1 = DMIN + J4 = J4 + 4 + J4P2 = J4 + 2*PP - 1 + Z( J4-2 ) = DNM1 + Z( J4P2 ) + IF( DNM1.LE.ZERO ) THEN + RETURN + ELSE + Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) + DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU + END IF + DMIN = MIN( DMIN, DN ) +* END IF * -* Unroll last two steps. -* - DNM2 = D - DMIN2 = DMIN - J4 = 4*( N0-2 ) - PP - J4P2 = J4 + 2*PP - 1 - Z( J4-2 ) = DNM2 + Z( J4P2 ) - Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) - DNM1 = Z( J4P2+2 )*( DNM2 / Z( J4-2 ) ) - TAU - DMIN = MIN( DMIN, DNM1 ) -* - DMIN1 = DMIN - J4 = J4 + 4 - J4P2 = J4 + 2*PP - 1 - Z( J4-2 ) = DNM1 + Z( J4P2 ) - Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) - DN = Z( J4P2+2 )*( DNM1 / Z( J4-2 ) ) - TAU - DMIN = MIN( DMIN, DN ) -* Z( J4+2 ) = DN Z( 4*N0-PP ) = EMIN RETURN
--- a/libcruft/lapack/dlasq6.f +++ b/libcruft/lapack/dlasq6.f @@ -1,10 +1,10 @@ - SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, DNM1, - $ DNM2 ) + SUBROUTINE DLASQ6( I0, N0, Z, PP, DMIN, DMIN1, DMIN2, DN, + $ DNM1, DNM2 ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. INTEGER I0, N0, PP @@ -16,7 +16,9 @@ * * Purpose * ======= -* DLASQ6 computes one dqds transform in ping-pong form. +* +* DLASQ6 computes one dqd (shift equal to zero) transform in +* ping-pong form, with protection against underflow and overflow. * * Arguments * ========= @@ -54,15 +56,15 @@ * * ===================================================================== * -* .. Parameters .. +* .. Parameter .. DOUBLE PRECISION ZERO PARAMETER ( ZERO = 0.0D0 ) * .. * .. Local Scalars .. INTEGER J4, J4P2 - DOUBLE PRECISION D, EMIN, SFMIN, TEMP + DOUBLE PRECISION D, EMIN, SAFMIN, TEMP * .. -* .. External Functions .. +* .. External Function .. DOUBLE PRECISION DLAMCH EXTERNAL DLAMCH * .. @@ -74,33 +76,55 @@ IF( ( N0-I0-1 ).LE.0 ) $ RETURN * - SFMIN = DLAMCH( 'Safe minimum' ) + SAFMIN = DLAMCH( 'Safe minimum' ) J4 = 4*I0 + PP - 3 - EMIN = Z( J4+4 ) + EMIN = Z( J4+4 ) D = Z( J4 ) DMIN = D * - DO 10 J4 = 4*I0 - PP, 4*( N0-3 ) - PP, 4 - J4P2 = J4 + 2*PP - 1 - Z( J4-2 ) = D + Z( J4P2 ) - IF( Z( J4-2 ).EQ.ZERO ) THEN - Z( J4 ) = ZERO - D = Z( J4P2+2 ) - DMIN = D - EMIN = ZERO - ELSE IF( SFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) ) THEN - TEMP = Z( J4P2+2 ) / Z( J4-2 ) - Z( J4 ) = Z( J4P2 )*TEMP - D = D*TEMP - ELSE - Z( J4 ) = Z( J4P2+2 )*( Z( J4P2 ) / Z( J4-2 ) ) - D = Z( J4P2+2 )*( D / Z( J4-2 ) ) - END IF - DMIN = MIN( DMIN, D ) - EMIN = MIN( EMIN, Z( J4 ) ) - 10 CONTINUE + IF( PP.EQ.0 ) THEN + DO 10 J4 = 4*I0, 4*( N0-3 ), 4 + Z( J4-2 ) = D + Z( J4-1 ) + IF( Z( J4-2 ).EQ.ZERO ) THEN + Z( J4 ) = ZERO + D = Z( J4+1 ) + DMIN = D + EMIN = ZERO + ELSE IF( SAFMIN*Z( J4+1 ).LT.Z( J4-2 ) .AND. + $ SAFMIN*Z( J4-2 ).LT.Z( J4+1 ) ) THEN + TEMP = Z( J4+1 ) / Z( J4-2 ) + Z( J4 ) = Z( J4-1 )*TEMP + D = D*TEMP + ELSE + Z( J4 ) = Z( J4+1 )*( Z( J4-1 ) / Z( J4-2 ) ) + D = Z( J4+1 )*( D / Z( J4-2 ) ) + END IF + DMIN = MIN( DMIN, D ) + EMIN = MIN( EMIN, Z( J4 ) ) + 10 CONTINUE + ELSE + DO 20 J4 = 4*I0, 4*( N0-3 ), 4 + Z( J4-3 ) = D + Z( J4 ) + IF( Z( J4-3 ).EQ.ZERO ) THEN + Z( J4-1 ) = ZERO + D = Z( J4+2 ) + DMIN = D + EMIN = ZERO + ELSE IF( SAFMIN*Z( J4+2 ).LT.Z( J4-3 ) .AND. + $ SAFMIN*Z( J4-3 ).LT.Z( J4+2 ) ) THEN + TEMP = Z( J4+2 ) / Z( J4-3 ) + Z( J4-1 ) = Z( J4 )*TEMP + D = D*TEMP + ELSE + Z( J4-1 ) = Z( J4+2 )*( Z( J4 ) / Z( J4-3 ) ) + D = Z( J4+2 )*( D / Z( J4-3 ) ) + END IF + DMIN = MIN( DMIN, D ) + EMIN = MIN( EMIN, Z( J4-1 ) ) + 20 CONTINUE + END IF * -* Unroll last two steps. +* Unroll last two steps. * DNM2 = D DMIN2 = DMIN @@ -112,7 +136,8 @@ DNM1 = Z( J4P2+2 ) DMIN = DNM1 EMIN = ZERO - ELSE IF( SFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) ) THEN + ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND. + $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN TEMP = Z( J4P2+2 ) / Z( J4-2 ) Z( J4 ) = Z( J4P2 )*TEMP DNM1 = DNM2*TEMP @@ -131,7 +156,8 @@ DN = Z( J4P2+2 ) DMIN = DN EMIN = ZERO - ELSE IF( SFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) ) THEN + ELSE IF( SAFMIN*Z( J4P2+2 ).LT.Z( J4-2 ) .AND. + $ SAFMIN*Z( J4-2 ).LT.Z( J4P2+2 ) ) THEN TEMP = Z( J4P2+2 ) / Z( J4-2 ) Z( J4 ) = Z( J4P2 )*TEMP DN = DNM1*TEMP
--- a/libcruft/lapack/zbdsqr.f +++ b/libcruft/lapack/zbdsqr.f @@ -4,7 +4,7 @@ * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. CHARACTER UPLO @@ -95,9 +95,7 @@ * The leading dimension of the array C. * LDC >= max(1,N) if NCC > 0; LDC >=1 if NCC = 0. * -* RWORK (workspace) DOUBLE PRECISION array, dimension -* 2*N if only singular values wanted (NCVT = NRU = NCC = 0) -* max( 1, 4*N-4 ) otherwise +* RWORK (workspace) DOUBLE PRECISION array, dimension (4*N) * * INFO (output) INTEGER * = 0: successful exit
--- a/libcruft/lapack/zgelss.f +++ b/libcruft/lapack/zgelss.f @@ -4,7 +4,7 @@ * -- LAPACK driver routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. INTEGER INFO, LDA, LDB, LWORK, M, N, NRHS, RANK @@ -92,7 +92,7 @@ * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * -* RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N)-1) +* RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N)) * * INFO (output) INTEGER * = 0: successful exit @@ -171,7 +171,7 @@ * * Path 1a - overdetermined, with many more rows than columns * -* Space needed for ZBDSQR is BDSPAC = 5*N-1 +* Space needed for ZBDSQR is BDSPAC = 5*N * MM = N MAXWRK = MAX( MAXWRK, N+N*ILAENV( 1, 'ZGEQRF', ' ', M, N, @@ -201,7 +201,7 @@ * Path 2a - underdetermined, with many more columns * than rows * -* Space needed for ZBDSQR is BDSPAC = 5*M-1 +* Space needed for ZBDSQR is BDSPAC = 5*M * MAXWRK = M + M*ILAENV( 1, 'ZGELQF', ' ', M, N, -1, -1 ) MAXWRK = MAX( MAXWRK, 3*M+M*M+2*M* @@ -221,7 +221,7 @@ * * Path 2 - underdetermined * -* Space needed for ZBDSQR is BDSPAC = 5*M-1 +* Space needed for ZBDSQR is BDSPAC = 5*M * MAXWRK = 2*M + ( N+M )*ILAENV( 1, 'ZGEBRD', ' ', M, N, $ -1, -1 )
--- a/libcruft/lapack/zgesvd.f +++ b/libcruft/lapack/zgesvd.f @@ -4,7 +4,7 @@ * -- LAPACK driver routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* June 30, 1999 +* October 31, 1999 * * .. Scalar Arguments .. CHARACTER JOBU, JOBVT @@ -119,8 +119,7 @@ * this value as the first entry of the WORK array, and no error * message related to LWORK is issued by XERBLA. * -* RWORK (workspace) DOUBLE PRECISION array, dimension -* (max(3*min(M,N),5*min(M,N)-4)) +* RWORK (workspace) DOUBLE PRECISION array, dimension (5*min(M,N)) * On exit, if INFO > 0, RWORK(1:MIN(M,N)-1) contains the * unconverged superdiagonal elements of an upper bidiagonal * matrix B whose diagonal is in S (not necessarily sorted). @@ -221,7 +220,7 @@ $ N.GT.0 ) THEN IF( M.GE.N ) THEN * -* Space needed for ZBDSQR is BDSPAC = MAX( 3*N, 5*N-4 ) +* Space needed for ZBDSQR is BDSPAC = 5*N * IF( M.GE.MNTHR ) THEN IF( WNTUN ) THEN @@ -383,7 +382,7 @@ END IF ELSE * -* Space needed for ZBDSQR is BDSPAC = MAX( 3*M, 5*M-4 ) +* Space needed for ZBDSQR is BDSPAC = 5*M * IF( N.GE.MNTHR ) THEN IF( WNTVN ) THEN
--- a/libcruft/lapack/zhetd2.f +++ b/libcruft/lapack/zhetd2.f @@ -3,7 +3,7 @@ * -- LAPACK routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University -* September 30, 1994 +* October 31, 1999 * * .. Scalar Arguments .. CHARACTER UPLO @@ -197,6 +197,8 @@ CALL ZHER2( UPLO, I, -ONE, A( 1, I+1 ), 1, TAU, 1, A, $ LDA ) * + ELSE + A( I, I ) = DBLE( A( I, I ) ) END IF A( I, I+1 ) = E( I ) D( I+1 ) = A( I+1, I+1 ) @@ -240,6 +242,8 @@ CALL ZHER2( UPLO, N-I, -ONE, A( I+1, I ), 1, TAU( I ), 1, $ A( I+1, I+1 ), LDA ) * + ELSE + A( I+1, I+1 ) = DBLE( A( I+1, I+1 ) ) END IF A( I+1, I ) = E( I ) D( I ) = A( I, I )