| @@ -60,12 +60,6 @@ | |||
| *> singular values which are less than RCOND times the largest singular | |||
| *> value. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -53,12 +53,6 @@ | |||
| *> | |||
| *> Note that the routine returns VT = V**H, not V. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -41,12 +41,6 @@ | |||
| *> a complex Hermitian band matrix A. If eigenvectors are desired, it | |||
| *> uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -47,12 +47,6 @@ | |||
| *> the reduction to tridiagonal. If eigenvectors are desired, it | |||
| *> uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -46,12 +46,6 @@ | |||
| *> and banded, and B is also positive definite. If eigenvectors are | |||
| *> desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -41,12 +41,6 @@ | |||
| *> complex Hermitian matrix A. If eigenvectors are desired, it uses a | |||
| *> divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -46,12 +46,6 @@ | |||
| *> the reduction to tridiagonal. If eigenvectors are desired, it uses a | |||
| *> divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -43,12 +43,6 @@ | |||
| *> B are assumed to be Hermitian and B is also positive definite. | |||
| *> If eigenvectors are desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -41,12 +41,6 @@ | |||
| *> a complex Hermitian matrix A in packed storage. If eigenvectors are | |||
| *> desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -44,12 +44,6 @@ | |||
| *> positive definite. | |||
| *> If eigenvectors are desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -18,7 +18,7 @@ | |||
| * Definition: | |||
| * =========== | |||
| * | |||
| * SUBROUTINE CLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA, | |||
| * SUBROUTINE CLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMBDA, | |||
| * Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR, | |||
| * GIVCOL, GIVNUM, INFO ) | |||
| * | |||
| @@ -29,7 +29,7 @@ | |||
| * .. Array Arguments .. | |||
| * INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), | |||
| * $ INDXQ( * ), PERM( * ) | |||
| * REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ), | |||
| * REAL D( * ), DLAMBDA( * ), GIVNUM( 2, * ), W( * ), | |||
| * $ Z( * ) | |||
| * COMPLEX Q( LDQ, * ), Q2( LDQ2, * ) | |||
| * .. | |||
| @@ -122,9 +122,9 @@ | |||
| *> destroyed during the updating process. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[out] DLAMDA | |||
| *> \param[out] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is REAL array, dimension (N) | |||
| *> DLAMBDA is REAL array, dimension (N) | |||
| *> Contains a copy of the first K eigenvalues which will be used | |||
| *> by SLAED3 to form the secular equation. | |||
| *> \endverbatim | |||
| @@ -222,7 +222,7 @@ | |||
| *> \ingroup complexOTHERcomputational | |||
| * | |||
| * ===================================================================== | |||
| SUBROUTINE CLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA, | |||
| SUBROUTINE CLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMBDA, | |||
| $ Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR, | |||
| $ GIVCOL, GIVNUM, INFO ) | |||
| * | |||
| @@ -237,7 +237,7 @@ | |||
| * .. Array Arguments .. | |||
| INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), | |||
| $ INDXQ( * ), PERM( * ) | |||
| REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ), | |||
| REAL D( * ), DLAMBDA( * ), GIVNUM( 2, * ), W( * ), | |||
| $ Z( * ) | |||
| COMPLEX Q( LDQ, * ), Q2( LDQ2, * ) | |||
| * .. | |||
| @@ -322,14 +322,14 @@ | |||
| INDXQ( I ) = INDXQ( I ) + CUTPNT | |||
| 20 CONTINUE | |||
| DO 30 I = 1, N | |||
| DLAMDA( I ) = D( INDXQ( I ) ) | |||
| DLAMBDA( I ) = D( INDXQ( I ) ) | |||
| W( I ) = Z( INDXQ( I ) ) | |||
| 30 CONTINUE | |||
| I = 1 | |||
| J = CUTPNT + 1 | |||
| CALL SLAMRG( N1, N2, DLAMDA, 1, 1, INDX ) | |||
| CALL SLAMRG( N1, N2, DLAMBDA, 1, 1, INDX ) | |||
| DO 40 I = 1, N | |||
| D( I ) = DLAMDA( INDX( I ) ) | |||
| D( I ) = DLAMBDA( INDX( I ) ) | |||
| Z( I ) = W( INDX( I ) ) | |||
| 40 CONTINUE | |||
| * | |||
| @@ -438,7 +438,7 @@ | |||
| ELSE | |||
| K = K + 1 | |||
| W( K ) = Z( JLAM ) | |||
| DLAMDA( K ) = D( JLAM ) | |||
| DLAMBDA( K ) = D( JLAM ) | |||
| INDXP( K ) = JLAM | |||
| JLAM = J | |||
| END IF | |||
| @@ -450,19 +450,19 @@ | |||
| * | |||
| K = K + 1 | |||
| W( K ) = Z( JLAM ) | |||
| DLAMDA( K ) = D( JLAM ) | |||
| DLAMBDA( K ) = D( JLAM ) | |||
| INDXP( K ) = JLAM | |||
| * | |||
| 100 CONTINUE | |||
| * | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMDA | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMBDA | |||
| * and Q2 respectively. The eigenvalues/vectors which were not | |||
| * deflated go into the first K slots of DLAMDA and Q2 respectively, | |||
| * deflated go into the first K slots of DLAMBDA and Q2 respectively, | |||
| * while those which were deflated go into the last N - K slots. | |||
| * | |||
| DO 110 J = 1, N | |||
| JP = INDXP( J ) | |||
| DLAMDA( J ) = D( JP ) | |||
| DLAMBDA( J ) = D( JP ) | |||
| PERM( J ) = INDXQ( INDX( JP ) ) | |||
| CALL CCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) | |||
| 110 CONTINUE | |||
| @@ -471,7 +471,7 @@ | |||
| * into the last N - K slots of D and Q respectively. | |||
| * | |||
| IF( K.LT.N ) THEN | |||
| CALL SCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL SCOPY( N-K, DLAMBDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL CLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, Q( 1, K+1 ), | |||
| $ LDQ ) | |||
| END IF | |||
| @@ -392,6 +392,11 @@ | |||
| $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN | |||
| RWORK( I ) = ZERO | |||
| ELSE | |||
| * | |||
| * Use calls to the subroutine SLAMC3 to enforce the | |||
| * parentheses (x+y)+z. The goal is to prevent | |||
| * optimizing compilers from doing x+(y+z). | |||
| * | |||
| RWORK( I ) = POLES( I, 2 )*Z( I ) / | |||
| $ ( SLAMC3( POLES( I, 2 ), DSIGJ )- | |||
| $ DIFLJ ) / ( POLES( I, 2 )+DJ ) | |||
| @@ -470,6 +475,11 @@ | |||
| IF( Z( J ).EQ.ZERO ) THEN | |||
| RWORK( I ) = ZERO | |||
| ELSE | |||
| * | |||
| * Use calls to the subroutine SLAMC3 to enforce the | |||
| * parentheses (x+y)+z. The goal is to prevent optimizing | |||
| * compilers from doing x+(y+z). | |||
| * | |||
| RWORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I+1, | |||
| $ 2 ) )-DIFR( I, 1 ) ) / | |||
| $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) | |||
| @@ -48,12 +48,6 @@ | |||
| *> problem; in this case a minimum norm solution is returned. | |||
| *> The actual singular values are returned in D in ascending order. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -43,12 +43,6 @@ | |||
| *> be found if CHETRD or CHPTRD or CHBTRD has been used to reduce this | |||
| *> matrix to tridiagonal form. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. See SLAED3 for details. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -45,13 +45,6 @@ | |||
| *> respectively. DBDSDC can be used to compute all singular values, | |||
| *> and optionally, singular vectors or singular vectors in compact form. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. See DLASD3 for details. | |||
| *> | |||
| *> The code currently calls DLASDQ if singular values only are desired. | |||
| *> However, it can be slightly modified to compute singular values | |||
| *> using the divide and conquer method. | |||
| @@ -59,12 +59,6 @@ | |||
| *> singular values which are less than RCOND times the largest singular | |||
| *> value. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -55,12 +55,6 @@ | |||
| *> | |||
| *> Note that the routine returns VT = V**T, not V. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -18,7 +18,7 @@ | |||
| * Definition: | |||
| * =========== | |||
| * | |||
| * SUBROUTINE DLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, | |||
| * SUBROUTINE DLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMBDA, W, | |||
| * Q2, INDX, INDXC, INDXP, COLTYP, INFO ) | |||
| * | |||
| * .. Scalar Arguments .. | |||
| @@ -28,7 +28,7 @@ | |||
| * .. Array Arguments .. | |||
| * INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ), | |||
| * $ INDXQ( * ) | |||
| * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), | |||
| * DOUBLE PRECISION D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ), | |||
| * $ W( * ), Z( * ) | |||
| * .. | |||
| * | |||
| @@ -123,9 +123,9 @@ | |||
| *> process. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[out] DLAMDA | |||
| *> \param[out] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is DOUBLE PRECISION array, dimension (N) | |||
| *> DLAMBDA is DOUBLE PRECISION array, dimension (N) | |||
| *> A copy of the first K eigenvalues which will be used by | |||
| *> DLAED3 to form the secular equation. | |||
| *> \endverbatim | |||
| @@ -148,7 +148,7 @@ | |||
| *> \param[out] INDX | |||
| *> \verbatim | |||
| *> INDX is INTEGER array, dimension (N) | |||
| *> The permutation used to sort the contents of DLAMDA into | |||
| *> The permutation used to sort the contents of DLAMBDA into | |||
| *> ascending order. | |||
| *> \endverbatim | |||
| *> | |||
| @@ -207,7 +207,7 @@ | |||
| *> Modified by Francoise Tisseur, University of Tennessee | |||
| *> | |||
| * ===================================================================== | |||
| SUBROUTINE DLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, | |||
| SUBROUTINE DLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMBDA, W, | |||
| $ Q2, INDX, INDXC, INDXP, COLTYP, INFO ) | |||
| * | |||
| * -- LAPACK computational routine -- | |||
| @@ -221,7 +221,7 @@ | |||
| * .. Array Arguments .. | |||
| INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ), | |||
| $ INDXQ( * ) | |||
| DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), | |||
| DOUBLE PRECISION D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ), | |||
| $ W( * ), Z( * ) | |||
| * .. | |||
| * | |||
| @@ -300,9 +300,9 @@ | |||
| * re-integrate the deflated parts from the last pass | |||
| * | |||
| DO 20 I = 1, N | |||
| DLAMDA( I ) = D( INDXQ( I ) ) | |||
| DLAMBDA( I ) = D( INDXQ( I ) ) | |||
| 20 CONTINUE | |||
| CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDXC ) | |||
| CALL DLAMRG( N1, N2, DLAMBDA, 1, 1, INDXC ) | |||
| DO 30 I = 1, N | |||
| INDX( I ) = INDXQ( INDXC( I ) ) | |||
| 30 CONTINUE | |||
| @@ -324,11 +324,11 @@ | |||
| DO 40 J = 1, N | |||
| I = INDX( J ) | |||
| CALL DCOPY( N, Q( 1, I ), 1, Q2( IQ2 ), 1 ) | |||
| DLAMDA( J ) = D( I ) | |||
| DLAMBDA( J ) = D( I ) | |||
| IQ2 = IQ2 + N | |||
| 40 CONTINUE | |||
| CALL DLACPY( 'A', N, N, Q2, N, Q, LDQ ) | |||
| CALL DCOPY( N, DLAMDA, 1, D, 1 ) | |||
| CALL DCOPY( N, DLAMBDA, 1, D, 1 ) | |||
| GO TO 190 | |||
| END IF | |||
| * | |||
| @@ -421,7 +421,7 @@ | |||
| PJ = NJ | |||
| ELSE | |||
| K = K + 1 | |||
| DLAMDA( K ) = D( PJ ) | |||
| DLAMBDA( K ) = D( PJ ) | |||
| W( K ) = Z( PJ ) | |||
| INDXP( K ) = PJ | |||
| PJ = NJ | |||
| @@ -433,7 +433,7 @@ | |||
| * Record the last eigenvalue. | |||
| * | |||
| K = K + 1 | |||
| DLAMDA( K ) = D( PJ ) | |||
| DLAMBDA( K ) = D( PJ ) | |||
| W( K ) = Z( PJ ) | |||
| INDXP( K ) = PJ | |||
| * | |||
| @@ -470,9 +470,9 @@ | |||
| PSM( CT ) = PSM( CT ) + 1 | |||
| 130 CONTINUE | |||
| * | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMDA | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMBDA | |||
| * and Q2 respectively. The eigenvalues/vectors which were not | |||
| * deflated go into the first K slots of DLAMDA and Q2 respectively, | |||
| * deflated go into the first K slots of DLAMBDA and Q2 respectively, | |||
| * while those which were deflated go into the last N - K slots. | |||
| * | |||
| I = 1 | |||
| @@ -18,7 +18,7 @@ | |||
| * Definition: | |||
| * =========== | |||
| * | |||
| * SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, | |||
| * SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMBDA, Q2, INDX, | |||
| * CTOT, W, S, INFO ) | |||
| * | |||
| * .. Scalar Arguments .. | |||
| @@ -27,7 +27,7 @@ | |||
| * .. | |||
| * .. Array Arguments .. | |||
| * INTEGER CTOT( * ), INDX( * ) | |||
| * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), | |||
| * DOUBLE PRECISION D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ), | |||
| * $ S( * ), W( * ) | |||
| * .. | |||
| * | |||
| @@ -44,12 +44,6 @@ | |||
| *> being combined by the matrix of eigenvectors of the K-by-K system | |||
| *> which is solved here. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -104,14 +98,12 @@ | |||
| *> RHO >= 0 required. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in,out] DLAMDA | |||
| *> \param[in] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is DOUBLE PRECISION array, dimension (K) | |||
| *> DLAMBDA is DOUBLE PRECISION array, dimension (K) | |||
| *> The first K elements of this array contain the old roots | |||
| *> of the deflated updating problem. These are the poles | |||
| *> of the secular equation. May be changed on output by | |||
| *> having lowest order bit set to zero on Cray X-MP, Cray Y-MP, | |||
| *> Cray-2, or Cray C-90, as described above. | |||
| *> of the secular equation. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in] Q2 | |||
| @@ -180,7 +172,7 @@ | |||
| *> Modified by Francoise Tisseur, University of Tennessee | |||
| *> | |||
| * ===================================================================== | |||
| SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, | |||
| SUBROUTINE DLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMBDA, Q2, INDX, | |||
| $ CTOT, W, S, INFO ) | |||
| * | |||
| * -- LAPACK computational routine -- | |||
| @@ -193,7 +185,7 @@ | |||
| * .. | |||
| * .. Array Arguments .. | |||
| INTEGER CTOT( * ), INDX( * ) | |||
| DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), | |||
| DOUBLE PRECISION D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ), | |||
| $ S( * ), W( * ) | |||
| * .. | |||
| * | |||
| @@ -208,8 +200,8 @@ | |||
| DOUBLE PRECISION TEMP | |||
| * .. | |||
| * .. External Functions .. | |||
| DOUBLE PRECISION DLAMC3, DNRM2 | |||
| EXTERNAL DLAMC3, DNRM2 | |||
| DOUBLE PRECISION DNRM2 | |||
| EXTERNAL DNRM2 | |||
| * .. | |||
| * .. External Subroutines .. | |||
| EXTERNAL DCOPY, DGEMM, DLACPY, DLAED4, DLASET, XERBLA | |||
| @@ -240,29 +232,9 @@ | |||
| IF( K.EQ.0 ) | |||
| $ RETURN | |||
| * | |||
| * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can | |||
| * be computed with high relative accuracy (barring over/underflow). | |||
| * This is a problem on machines without a guard digit in | |||
| * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). | |||
| * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), | |||
| * which on any of these machines zeros out the bottommost | |||
| * bit of DLAMDA(I) if it is 1; this makes the subsequent | |||
| * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation | |||
| * occurs. On binary machines with a guard digit (almost all | |||
| * machines) it does not change DLAMDA(I) at all. On hexadecimal | |||
| * and decimal machines with a guard digit, it slightly | |||
| * changes the bottommost bits of DLAMDA(I). It does not account | |||
| * for hexadecimal or decimal machines without guard digits | |||
| * (we know of none). We use a subroutine call to compute | |||
| * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating | |||
| * this code. | |||
| * | |||
| DO 10 I = 1, K | |||
| DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) | |||
| 10 CONTINUE | |||
| * | |||
| DO 20 J = 1, K | |||
| CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) | |||
| CALL DLAED4( K, J, DLAMBDA, W, Q( 1, J ), RHO, D( J ), INFO ) | |||
| * | |||
| * If the zero finder fails, the computation is terminated. | |||
| * | |||
| @@ -293,10 +265,10 @@ | |||
| CALL DCOPY( K, Q, LDQ+1, W, 1 ) | |||
| DO 60 J = 1, K | |||
| DO 40 I = 1, J - 1 | |||
| W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) | |||
| W( I ) = W( I )*( Q( I, J )/( DLAMBDA( I )-DLAMBDA( J ) ) ) | |||
| 40 CONTINUE | |||
| DO 50 I = J + 1, K | |||
| W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) | |||
| W( I ) = W( I )*( Q( I, J )/( DLAMBDA( I )-DLAMBDA( J ) ) ) | |||
| 50 CONTINUE | |||
| 60 CONTINUE | |||
| DO 70 I = 1, K | |||
| @@ -19,7 +19,7 @@ | |||
| * =========== | |||
| * | |||
| * SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, | |||
| * CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, | |||
| * CUTPNT, Z, DLAMBDA, Q2, LDQ2, W, PERM, GIVPTR, | |||
| * GIVCOL, GIVNUM, INDXP, INDX, INFO ) | |||
| * | |||
| * .. Scalar Arguments .. | |||
| @@ -30,7 +30,7 @@ | |||
| * .. Array Arguments .. | |||
| * INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), | |||
| * $ INDXQ( * ), PERM( * ) | |||
| * DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), | |||
| * DOUBLE PRECISION D( * ), DLAMBDA( * ), GIVNUM( 2, * ), | |||
| * $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) | |||
| * .. | |||
| * | |||
| @@ -141,9 +141,9 @@ | |||
| *> process. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[out] DLAMDA | |||
| *> \param[out] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is DOUBLE PRECISION array, dimension (N) | |||
| *> DLAMBDA is DOUBLE PRECISION array, dimension (N) | |||
| *> A copy of the first K eigenvalues which will be used by | |||
| *> DLAED3 to form the secular equation. | |||
| *> \endverbatim | |||
| @@ -238,7 +238,7 @@ | |||
| * | |||
| * ===================================================================== | |||
| SUBROUTINE DLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, | |||
| $ CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, | |||
| $ CUTPNT, Z, DLAMBDA, Q2, LDQ2, W, PERM, GIVPTR, | |||
| $ GIVCOL, GIVNUM, INDXP, INDX, INFO ) | |||
| * | |||
| * -- LAPACK computational routine -- | |||
| @@ -253,7 +253,7 @@ | |||
| * .. Array Arguments .. | |||
| INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), | |||
| $ INDXQ( * ), PERM( * ) | |||
| DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), | |||
| DOUBLE PRECISION D( * ), DLAMBDA( * ), GIVNUM( 2, * ), | |||
| $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) | |||
| * .. | |||
| * | |||
| @@ -339,14 +339,14 @@ | |||
| INDXQ( I ) = INDXQ( I ) + CUTPNT | |||
| 20 CONTINUE | |||
| DO 30 I = 1, N | |||
| DLAMDA( I ) = D( INDXQ( I ) ) | |||
| DLAMBDA( I ) = D( INDXQ( I ) ) | |||
| W( I ) = Z( INDXQ( I ) ) | |||
| 30 CONTINUE | |||
| I = 1 | |||
| J = CUTPNT + 1 | |||
| CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX ) | |||
| CALL DLAMRG( N1, N2, DLAMBDA, 1, 1, INDX ) | |||
| DO 40 I = 1, N | |||
| D( I ) = DLAMDA( INDX( I ) ) | |||
| D( I ) = DLAMBDA( INDX( I ) ) | |||
| Z( I ) = W( INDX( I ) ) | |||
| 40 CONTINUE | |||
| * | |||
| @@ -464,7 +464,7 @@ | |||
| ELSE | |||
| K = K + 1 | |||
| W( K ) = Z( JLAM ) | |||
| DLAMDA( K ) = D( JLAM ) | |||
| DLAMBDA( K ) = D( JLAM ) | |||
| INDXP( K ) = JLAM | |||
| JLAM = J | |||
| END IF | |||
| @@ -476,26 +476,26 @@ | |||
| * | |||
| K = K + 1 | |||
| W( K ) = Z( JLAM ) | |||
| DLAMDA( K ) = D( JLAM ) | |||
| DLAMBDA( K ) = D( JLAM ) | |||
| INDXP( K ) = JLAM | |||
| * | |||
| 110 CONTINUE | |||
| * | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMDA | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMBDA | |||
| * and Q2 respectively. The eigenvalues/vectors which were not | |||
| * deflated go into the first K slots of DLAMDA and Q2 respectively, | |||
| * deflated go into the first K slots of DLAMBDA and Q2 respectively, | |||
| * while those which were deflated go into the last N - K slots. | |||
| * | |||
| IF( ICOMPQ.EQ.0 ) THEN | |||
| DO 120 J = 1, N | |||
| JP = INDXP( J ) | |||
| DLAMDA( J ) = D( JP ) | |||
| DLAMBDA( J ) = D( JP ) | |||
| PERM( J ) = INDXQ( INDX( JP ) ) | |||
| 120 CONTINUE | |||
| ELSE | |||
| DO 130 J = 1, N | |||
| JP = INDXP( J ) | |||
| DLAMDA( J ) = D( JP ) | |||
| DLAMBDA( J ) = D( JP ) | |||
| PERM( J ) = INDXQ( INDX( JP ) ) | |||
| CALL DCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) | |||
| 130 CONTINUE | |||
| @@ -506,9 +506,9 @@ | |||
| * | |||
| IF( K.LT.N ) THEN | |||
| IF( ICOMPQ.EQ.0 ) THEN | |||
| CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL DCOPY( N-K, DLAMBDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| ELSE | |||
| CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL DCOPY( N-K, DLAMBDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL DLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, | |||
| $ Q( 1, K+1 ), LDQ ) | |||
| END IF | |||
| @@ -18,15 +18,15 @@ | |||
| * Definition: | |||
| * =========== | |||
| * | |||
| * SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, | |||
| * S, LDS, INFO ) | |||
| * SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMBDA, | |||
| * W, S, LDS, INFO ) | |||
| * | |||
| * .. Scalar Arguments .. | |||
| * INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N | |||
| * DOUBLE PRECISION RHO | |||
| * .. | |||
| * .. Array Arguments .. | |||
| * DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), | |||
| * DOUBLE PRECISION D( * ), DLAMBDA( * ), Q( LDQ, * ), S( LDS, * ), | |||
| * $ W( * ) | |||
| * .. | |||
| * | |||
| @@ -96,9 +96,9 @@ | |||
| *> RHO >= 0 required. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in] DLAMDA | |||
| *> \param[in] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is DOUBLE PRECISION array, dimension (K) | |||
| *> DLAMBDA is DOUBLE PRECISION array, dimension (K) | |||
| *> The first K elements of this array contain the old roots | |||
| *> of the deflated updating problem. These are the poles | |||
| *> of the secular equation. | |||
| @@ -151,8 +151,8 @@ | |||
| *> at Berkeley, USA | |||
| * | |||
| * ===================================================================== | |||
| SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, | |||
| $ S, LDS, INFO ) | |||
| SUBROUTINE DLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMBDA, | |||
| $ W, S, LDS, INFO ) | |||
| * | |||
| * -- LAPACK computational routine -- | |||
| * -- LAPACK is a software package provided by Univ. of Tennessee, -- | |||
| @@ -163,7 +163,7 @@ | |||
| DOUBLE PRECISION RHO | |||
| * .. | |||
| * .. Array Arguments .. | |||
| DOUBLE PRECISION D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), | |||
| DOUBLE PRECISION D( * ), DLAMBDA( * ), Q( LDQ, * ), S( LDS, * ), | |||
| $ W( * ) | |||
| * .. | |||
| * | |||
| @@ -174,8 +174,8 @@ | |||
| DOUBLE PRECISION TEMP | |||
| * .. | |||
| * .. External Functions .. | |||
| DOUBLE PRECISION DLAMC3, DNRM2 | |||
| EXTERNAL DLAMC3, DNRM2 | |||
| DOUBLE PRECISION DNRM2 | |||
| EXTERNAL DNRM2 | |||
| * .. | |||
| * .. External Subroutines .. | |||
| EXTERNAL DCOPY, DLAED4, XERBLA | |||
| @@ -212,30 +212,9 @@ | |||
| * | |||
| IF( K.EQ.0 ) | |||
| $ RETURN | |||
| * | |||
| * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can | |||
| * be computed with high relative accuracy (barring over/underflow). | |||
| * This is a problem on machines without a guard digit in | |||
| * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). | |||
| * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), | |||
| * which on any of these machines zeros out the bottommost | |||
| * bit of DLAMDA(I) if it is 1; this makes the subsequent | |||
| * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation | |||
| * occurs. On binary machines with a guard digit (almost all | |||
| * machines) it does not change DLAMDA(I) at all. On hexadecimal | |||
| * and decimal machines with a guard digit, it slightly | |||
| * changes the bottommost bits of DLAMDA(I). It does not account | |||
| * for hexadecimal or decimal machines without guard digits | |||
| * (we know of none). We use a subroutine call to compute | |||
| * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating | |||
| * this code. | |||
| * | |||
| DO 10 I = 1, N | |||
| DLAMDA( I ) = DLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) | |||
| 10 CONTINUE | |||
| * | |||
| DO 20 J = KSTART, KSTOP | |||
| CALL DLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) | |||
| CALL DLAED4( K, J, DLAMBDA, W, Q( 1, J ), RHO, D( J ), INFO ) | |||
| * | |||
| * If the zero finder fails, the computation is terminated. | |||
| * | |||
| @@ -261,10 +240,10 @@ | |||
| CALL DCOPY( K, Q, LDQ+1, W, 1 ) | |||
| DO 70 J = 1, K | |||
| DO 50 I = 1, J - 1 | |||
| W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) | |||
| W( I ) = W( I )*( Q( I, J )/( DLAMBDA( I )-DLAMBDA( J ) ) ) | |||
| 50 CONTINUE | |||
| DO 60 I = J + 1, K | |||
| W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) | |||
| W( I ) = W( I )*( Q( I, J )/( DLAMBDA( I )-DLAMBDA( J ) ) ) | |||
| 60 CONTINUE | |||
| 70 CONTINUE | |||
| DO 80 I = 1, K | |||
| @@ -389,6 +389,11 @@ | |||
| $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN | |||
| WORK( I ) = ZERO | |||
| ELSE | |||
| * | |||
| * Use calls to the subroutine DLAMC3 to enforce the | |||
| * parentheses (x+y)+z. The goal is to prevent | |||
| * optimizing compilers from doing x+(y+z). | |||
| * | |||
| WORK( I ) = POLES( I, 2 )*Z( I ) / | |||
| $ ( DLAMC3( POLES( I, 2 ), DSIGJ )- | |||
| $ DIFLJ ) / ( POLES( I, 2 )+DJ ) | |||
| @@ -440,6 +445,11 @@ | |||
| IF( Z( J ).EQ.ZERO ) THEN | |||
| WORK( I ) = ZERO | |||
| ELSE | |||
| * | |||
| * Use calls to the subroutine DLAMC3 to enforce the | |||
| * parentheses (x+y)+z. The goal is to prevent | |||
| * optimizing compilers from doing x+(y+z). | |||
| * | |||
| WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1, | |||
| $ 2 ) )-DIFR( I, 1 ) ) / | |||
| $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) | |||
| @@ -47,12 +47,6 @@ | |||
| *> problem; in this case a minimum norm solution is returned. | |||
| *> The actual singular values are returned in D in ascending order. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -93,9 +93,7 @@ | |||
| *> infinite. | |||
| *> | |||
| *> Overflow will not occur unless the largest singular value itself | |||
| *> overflows, or is within a few ulps of overflow. (On machines with | |||
| *> partial overflow, like the Cray, overflow may occur if the largest | |||
| *> singular value is within a factor of 2 of overflow.) | |||
| *> overflows, or is within a few ulps of overflow. | |||
| *> | |||
| *> Underflow is harmless if underflow is gradual. Otherwise, results | |||
| *> may correspond to a matrix modified by perturbations of size near | |||
| @@ -44,13 +44,6 @@ | |||
| *> appropriate calls to DLASD4 and then updates the singular | |||
| *> vectors by matrix multiplication. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> | |||
| *> DLASD3 is called from DLASD1. | |||
| *> \endverbatim | |||
| * | |||
| @@ -103,7 +96,7 @@ | |||
| *> The leading dimension of the array Q. LDQ >= K. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in,out] DSIGMA | |||
| *> \param[in] DSIGMA | |||
| *> \verbatim | |||
| *> DSIGMA is DOUBLE PRECISION array, dimension(K) | |||
| *> The first K elements of this array contain the old roots | |||
| @@ -249,8 +242,8 @@ | |||
| DOUBLE PRECISION RHO, TEMP | |||
| * .. | |||
| * .. External Functions .. | |||
| DOUBLE PRECISION DLAMC3, DNRM2 | |||
| EXTERNAL DLAMC3, DNRM2 | |||
| DOUBLE PRECISION DNRM2 | |||
| EXTERNAL DNRM2 | |||
| * .. | |||
| * .. External Subroutines .. | |||
| EXTERNAL DCOPY, DGEMM, DLACPY, DLASCL, DLASD4, XERBLA | |||
| @@ -310,27 +303,6 @@ | |||
| RETURN | |||
| END IF | |||
| * | |||
| * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can | |||
| * be computed with high relative accuracy (barring over/underflow). | |||
| * This is a problem on machines without a guard digit in | |||
| * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). | |||
| * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), | |||
| * which on any of these machines zeros out the bottommost | |||
| * bit of DSIGMA(I) if it is 1; this makes the subsequent | |||
| * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation | |||
| * occurs. On binary machines with a guard digit (almost all | |||
| * machines) it does not change DSIGMA(I) at all. On hexadecimal | |||
| * and decimal machines with a guard digit, it slightly | |||
| * changes the bottommost bits of DSIGMA(I). It does not account | |||
| * for hexadecimal or decimal machines without guard digits | |||
| * (we know of none). We use a subroutine call to compute | |||
| * 2*DSIGMA(I) to prevent optimizing compilers from eliminating | |||
| * this code. | |||
| * | |||
| DO 20 I = 1, K | |||
| DSIGMA( I ) = DLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) | |||
| 20 CONTINUE | |||
| * | |||
| * Keep a copy of Z. | |||
| * | |||
| CALL DCOPY( K, Z, 1, Q, 1 ) | |||
| @@ -121,14 +121,12 @@ | |||
| *> The leading dimension of DIFR, must be at least K. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in,out] DSIGMA | |||
| *> \param[in] DSIGMA | |||
| *> \verbatim | |||
| *> DSIGMA is DOUBLE PRECISION array, dimension ( K ) | |||
| *> On entry, the first K elements of this array contain the old | |||
| *> roots of the deflated updating problem. These are the poles | |||
| *> of the secular equation. | |||
| *> On exit, the elements of DSIGMA may be very slightly altered | |||
| *> in value. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[out] WORK | |||
| @@ -227,27 +225,6 @@ | |||
| RETURN | |||
| END IF | |||
| * | |||
| * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can | |||
| * be computed with high relative accuracy (barring over/underflow). | |||
| * This is a problem on machines without a guard digit in | |||
| * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). | |||
| * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), | |||
| * which on any of these machines zeros out the bottommost | |||
| * bit of DSIGMA(I) if it is 1; this makes the subsequent | |||
| * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation | |||
| * occurs. On binary machines with a guard digit (almost all | |||
| * machines) it does not change DSIGMA(I) at all. On hexadecimal | |||
| * and decimal machines with a guard digit, it slightly | |||
| * changes the bottommost bits of DSIGMA(I). It does not account | |||
| * for hexadecimal or decimal machines without guard digits | |||
| * (we know of none). We use a subroutine call to compute | |||
| * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating | |||
| * this code. | |||
| * | |||
| DO 10 I = 1, K | |||
| DSIGMA( I ) = DLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) | |||
| 10 CONTINUE | |||
| * | |||
| * Book keeping. | |||
| * | |||
| IWK1 = 1 | |||
| @@ -312,6 +289,11 @@ | |||
| DSIGJP = -DSIGMA( J+1 ) | |||
| END IF | |||
| WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ ) | |||
| * | |||
| * Use calls to the subroutine DLAMC3 to enforce the parentheses | |||
| * (x+y)+z. The goal is to prevent optimizing compilers | |||
| * from doing x+(y+z). | |||
| * | |||
| DO 60 I = 1, J - 1 | |||
| WORK( I ) = Z( I ) / ( DLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ ) | |||
| $ / ( DSIGMA( I )+DJ ) | |||
| @@ -124,9 +124,7 @@ | |||
| *> infinite. | |||
| *> | |||
| *> Overflow will not occur unless the largest singular value itself | |||
| *> overflows or is within a few ulps of overflow. (On machines with | |||
| *> partial overflow, like the Cray, overflow may occur if the largest | |||
| *> singular value is within a factor of 2 of overflow.) | |||
| *> overflows or is within a few ulps of overflow. | |||
| *> | |||
| *> Underflow is harmless if underflow is gradual. Otherwise, results | |||
| *> may correspond to a matrix modified by perturbations of size near | |||
| @@ -40,12 +40,6 @@ | |||
| *> a real symmetric band matrix A. If eigenvectors are desired, it uses | |||
| *> a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -45,12 +45,6 @@ | |||
| *> the reduction to tridiagonal. If eigenvectors are desired, it uses | |||
| *> a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -43,12 +43,6 @@ | |||
| *> banded, and B is also positive definite. If eigenvectors are | |||
| *> desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -40,12 +40,6 @@ | |||
| *> of a real symmetric matrix A in packed storage. If eigenvectors are | |||
| *> desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -44,12 +44,6 @@ | |||
| *> positive definite. | |||
| *> If eigenvectors are desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -42,12 +42,6 @@ | |||
| *> found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this | |||
| *> matrix to tridiagonal form. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. See DLAED3 for details. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -40,12 +40,6 @@ | |||
| *> real symmetric tridiagonal matrix. If eigenvectors are desired, it | |||
| *> uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -40,13 +40,6 @@ | |||
| *> real symmetric matrix A. If eigenvectors are desired, it uses a | |||
| *> divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> | |||
| *> Because of large use of BLAS of level 3, DSYEVD needs N**2 more | |||
| *> workspace than DSYEVX. | |||
| *> \endverbatim | |||
| @@ -45,12 +45,6 @@ | |||
| *> the reduction to tridiagonal. If eigenvectors are desired, it uses a | |||
| *> divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -42,12 +42,6 @@ | |||
| *> B are assumed to be symmetric and B is also positive definite. | |||
| *> If eigenvectors are desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -45,13 +45,6 @@ | |||
| *> respectively. SBDSDC can be used to compute all singular values, | |||
| *> and optionally, singular vectors or singular vectors in compact form. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. See SLASD3 for details. | |||
| *> | |||
| *> The code currently calls SLASDQ if singular values only are desired. | |||
| *> However, it can be slightly modified to compute singular values | |||
| *> using the divide and conquer method. | |||
| @@ -59,12 +59,6 @@ | |||
| *> singular values which are less than RCOND times the largest singular | |||
| *> value. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -55,12 +55,6 @@ | |||
| *> | |||
| *> Note that the routine returns VT = V**T, not V. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -18,7 +18,7 @@ | |||
| * Definition: | |||
| * =========== | |||
| * | |||
| * SUBROUTINE SLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, | |||
| * SUBROUTINE SLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMBDA, W, | |||
| * Q2, INDX, INDXC, INDXP, COLTYP, INFO ) | |||
| * | |||
| * .. Scalar Arguments .. | |||
| @@ -28,7 +28,7 @@ | |||
| * .. Array Arguments .. | |||
| * INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ), | |||
| * $ INDXQ( * ) | |||
| * REAL D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), | |||
| * REAL D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ), | |||
| * $ W( * ), Z( * ) | |||
| * .. | |||
| * | |||
| @@ -123,9 +123,9 @@ | |||
| *> process. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[out] DLAMDA | |||
| *> \param[out] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is REAL array, dimension (N) | |||
| *> DLAMBDA is REAL array, dimension (N) | |||
| *> A copy of the first K eigenvalues which will be used by | |||
| *> SLAED3 to form the secular equation. | |||
| *> \endverbatim | |||
| @@ -148,7 +148,7 @@ | |||
| *> \param[out] INDX | |||
| *> \verbatim | |||
| *> INDX is INTEGER array, dimension (N) | |||
| *> The permutation used to sort the contents of DLAMDA into | |||
| *> The permutation used to sort the contents of DLAMBDA into | |||
| *> ascending order. | |||
| *> \endverbatim | |||
| *> | |||
| @@ -207,7 +207,7 @@ | |||
| *> Modified by Francoise Tisseur, University of Tennessee | |||
| *> | |||
| * ===================================================================== | |||
| SUBROUTINE SLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMDA, W, | |||
| SUBROUTINE SLAED2( K, N, N1, D, Q, LDQ, INDXQ, RHO, Z, DLAMBDA, W, | |||
| $ Q2, INDX, INDXC, INDXP, COLTYP, INFO ) | |||
| * | |||
| * -- LAPACK computational routine -- | |||
| @@ -221,7 +221,7 @@ | |||
| * .. Array Arguments .. | |||
| INTEGER COLTYP( * ), INDX( * ), INDXC( * ), INDXP( * ), | |||
| $ INDXQ( * ) | |||
| REAL D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), | |||
| REAL D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ), | |||
| $ W( * ), Z( * ) | |||
| * .. | |||
| * | |||
| @@ -300,9 +300,9 @@ | |||
| * re-integrate the deflated parts from the last pass | |||
| * | |||
| DO 20 I = 1, N | |||
| DLAMDA( I ) = D( INDXQ( I ) ) | |||
| DLAMBDA( I ) = D( INDXQ( I ) ) | |||
| 20 CONTINUE | |||
| CALL SLAMRG( N1, N2, DLAMDA, 1, 1, INDXC ) | |||
| CALL SLAMRG( N1, N2, DLAMBDA, 1, 1, INDXC ) | |||
| DO 30 I = 1, N | |||
| INDX( I ) = INDXQ( INDXC( I ) ) | |||
| 30 CONTINUE | |||
| @@ -324,11 +324,11 @@ | |||
| DO 40 J = 1, N | |||
| I = INDX( J ) | |||
| CALL SCOPY( N, Q( 1, I ), 1, Q2( IQ2 ), 1 ) | |||
| DLAMDA( J ) = D( I ) | |||
| DLAMBDA( J ) = D( I ) | |||
| IQ2 = IQ2 + N | |||
| 40 CONTINUE | |||
| CALL SLACPY( 'A', N, N, Q2, N, Q, LDQ ) | |||
| CALL SCOPY( N, DLAMDA, 1, D, 1 ) | |||
| CALL SCOPY( N, DLAMBDA, 1, D, 1 ) | |||
| GO TO 190 | |||
| END IF | |||
| * | |||
| @@ -421,7 +421,7 @@ | |||
| PJ = NJ | |||
| ELSE | |||
| K = K + 1 | |||
| DLAMDA( K ) = D( PJ ) | |||
| DLAMBDA( K ) = D( PJ ) | |||
| W( K ) = Z( PJ ) | |||
| INDXP( K ) = PJ | |||
| PJ = NJ | |||
| @@ -433,7 +433,7 @@ | |||
| * Record the last eigenvalue. | |||
| * | |||
| K = K + 1 | |||
| DLAMDA( K ) = D( PJ ) | |||
| DLAMBDA( K ) = D( PJ ) | |||
| W( K ) = Z( PJ ) | |||
| INDXP( K ) = PJ | |||
| * | |||
| @@ -470,9 +470,9 @@ | |||
| PSM( CT ) = PSM( CT ) + 1 | |||
| 130 CONTINUE | |||
| * | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMDA | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMBDA | |||
| * and Q2 respectively. The eigenvalues/vectors which were not | |||
| * deflated go into the first K slots of DLAMDA and Q2 respectively, | |||
| * deflated go into the first K slots of DLAMBDA and Q2 respectively, | |||
| * while those which were deflated go into the last N - K slots. | |||
| * | |||
| I = 1 | |||
| @@ -18,7 +18,7 @@ | |||
| * Definition: | |||
| * =========== | |||
| * | |||
| * SUBROUTINE SLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, | |||
| * SUBROUTINE SLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMBDA, Q2, INDX, | |||
| * CTOT, W, S, INFO ) | |||
| * | |||
| * .. Scalar Arguments .. | |||
| @@ -27,7 +27,7 @@ | |||
| * .. | |||
| * .. Array Arguments .. | |||
| * INTEGER CTOT( * ), INDX( * ) | |||
| * REAL D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), | |||
| * REAL D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ), | |||
| * $ S( * ), W( * ) | |||
| * .. | |||
| * | |||
| @@ -44,12 +44,6 @@ | |||
| *> being combined by the matrix of eigenvectors of the K-by-K system | |||
| *> which is solved here. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -104,14 +98,12 @@ | |||
| *> RHO >= 0 required. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in,out] DLAMDA | |||
| *> \param[in] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is REAL array, dimension (K) | |||
| *> DLAMBDA is REAL array, dimension (K) | |||
| *> The first K elements of this array contain the old roots | |||
| *> of the deflated updating problem. These are the poles | |||
| *> of the secular equation. May be changed on output by | |||
| *> having lowest order bit set to zero on Cray X-MP, Cray Y-MP, | |||
| *> Cray-2, or Cray C-90, as described above. | |||
| *> of the secular equation. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in] Q2 | |||
| @@ -180,7 +172,7 @@ | |||
| *> Modified by Francoise Tisseur, University of Tennessee | |||
| *> | |||
| * ===================================================================== | |||
| SUBROUTINE SLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMDA, Q2, INDX, | |||
| SUBROUTINE SLAED3( K, N, N1, D, Q, LDQ, RHO, DLAMBDA, Q2, INDX, | |||
| $ CTOT, W, S, INFO ) | |||
| * | |||
| * -- LAPACK computational routine -- | |||
| @@ -193,7 +185,7 @@ | |||
| * .. | |||
| * .. Array Arguments .. | |||
| INTEGER CTOT( * ), INDX( * ) | |||
| REAL D( * ), DLAMDA( * ), Q( LDQ, * ), Q2( * ), | |||
| REAL D( * ), DLAMBDA( * ), Q( LDQ, * ), Q2( * ), | |||
| $ S( * ), W( * ) | |||
| * .. | |||
| * | |||
| @@ -208,8 +200,8 @@ | |||
| REAL TEMP | |||
| * .. | |||
| * .. External Functions .. | |||
| REAL SLAMC3, SNRM2 | |||
| EXTERNAL SLAMC3, SNRM2 | |||
| REAL SNRM2 | |||
| EXTERNAL SNRM2 | |||
| * .. | |||
| * .. External Subroutines .. | |||
| EXTERNAL SCOPY, SGEMM, SLACPY, SLAED4, SLASET, XERBLA | |||
| @@ -239,30 +231,9 @@ | |||
| * | |||
| IF( K.EQ.0 ) | |||
| $ RETURN | |||
| * | |||
| * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can | |||
| * be computed with high relative accuracy (barring over/underflow). | |||
| * This is a problem on machines without a guard digit in | |||
| * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). | |||
| * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), | |||
| * which on any of these machines zeros out the bottommost | |||
| * bit of DLAMDA(I) if it is 1; this makes the subsequent | |||
| * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation | |||
| * occurs. On binary machines with a guard digit (almost all | |||
| * machines) it does not change DLAMDA(I) at all. On hexadecimal | |||
| * and decimal machines with a guard digit, it slightly | |||
| * changes the bottommost bits of DLAMDA(I). It does not account | |||
| * for hexadecimal or decimal machines without guard digits | |||
| * (we know of none). We use a subroutine call to compute | |||
| * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating | |||
| * this code. | |||
| * | |||
| DO 10 I = 1, K | |||
| DLAMDA( I ) = SLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) | |||
| 10 CONTINUE | |||
| * | |||
| DO 20 J = 1, K | |||
| CALL SLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) | |||
| CALL SLAED4( K, J, DLAMBDA, W, Q( 1, J ), RHO, D( J ), INFO ) | |||
| * | |||
| * If the zero finder fails, the computation is terminated. | |||
| * | |||
| @@ -293,10 +264,10 @@ | |||
| CALL SCOPY( K, Q, LDQ+1, W, 1 ) | |||
| DO 60 J = 1, K | |||
| DO 40 I = 1, J - 1 | |||
| W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) | |||
| W( I ) = W( I )*( Q( I, J )/( DLAMBDA( I )-DLAMBDA( J ) ) ) | |||
| 40 CONTINUE | |||
| DO 50 I = J + 1, K | |||
| W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) | |||
| W( I ) = W( I )*( Q( I, J )/( DLAMBDA( I )-DLAMBDA( J ) ) ) | |||
| 50 CONTINUE | |||
| 60 CONTINUE | |||
| DO 70 I = 1, K | |||
| @@ -19,7 +19,7 @@ | |||
| * =========== | |||
| * | |||
| * SUBROUTINE SLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, | |||
| * CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, | |||
| * CUTPNT, Z, DLAMBDA, Q2, LDQ2, W, PERM, GIVPTR, | |||
| * GIVCOL, GIVNUM, INDXP, INDX, INFO ) | |||
| * | |||
| * .. Scalar Arguments .. | |||
| @@ -30,7 +30,7 @@ | |||
| * .. Array Arguments .. | |||
| * INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), | |||
| * $ INDXQ( * ), PERM( * ) | |||
| * REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ), | |||
| * REAL D( * ), DLAMBDA( * ), GIVNUM( 2, * ), | |||
| * $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) | |||
| * .. | |||
| * | |||
| @@ -141,9 +141,9 @@ | |||
| *> process. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[out] DLAMDA | |||
| *> \param[out] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is REAL array, dimension (N) | |||
| *> DLAMBDA is REAL array, dimension (N) | |||
| *> A copy of the first K eigenvalues which will be used by | |||
| *> SLAED3 to form the secular equation. | |||
| *> \endverbatim | |||
| @@ -238,7 +238,7 @@ | |||
| * | |||
| * ===================================================================== | |||
| SUBROUTINE SLAED8( ICOMPQ, K, N, QSIZ, D, Q, LDQ, INDXQ, RHO, | |||
| $ CUTPNT, Z, DLAMDA, Q2, LDQ2, W, PERM, GIVPTR, | |||
| $ CUTPNT, Z, DLAMBDA, Q2, LDQ2, W, PERM, GIVPTR, | |||
| $ GIVCOL, GIVNUM, INDXP, INDX, INFO ) | |||
| * | |||
| * -- LAPACK computational routine -- | |||
| @@ -253,7 +253,7 @@ | |||
| * .. Array Arguments .. | |||
| INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), | |||
| $ INDXQ( * ), PERM( * ) | |||
| REAL D( * ), DLAMDA( * ), GIVNUM( 2, * ), | |||
| REAL D( * ), DLAMBDA( * ), GIVNUM( 2, * ), | |||
| $ Q( LDQ, * ), Q2( LDQ2, * ), W( * ), Z( * ) | |||
| * .. | |||
| * | |||
| @@ -339,14 +339,14 @@ | |||
| INDXQ( I ) = INDXQ( I ) + CUTPNT | |||
| 20 CONTINUE | |||
| DO 30 I = 1, N | |||
| DLAMDA( I ) = D( INDXQ( I ) ) | |||
| DLAMBDA( I ) = D( INDXQ( I ) ) | |||
| W( I ) = Z( INDXQ( I ) ) | |||
| 30 CONTINUE | |||
| I = 1 | |||
| J = CUTPNT + 1 | |||
| CALL SLAMRG( N1, N2, DLAMDA, 1, 1, INDX ) | |||
| CALL SLAMRG( N1, N2, DLAMBDA, 1, 1, INDX ) | |||
| DO 40 I = 1, N | |||
| D( I ) = DLAMDA( INDX( I ) ) | |||
| D( I ) = DLAMBDA( INDX( I ) ) | |||
| Z( I ) = W( INDX( I ) ) | |||
| 40 CONTINUE | |||
| * | |||
| @@ -464,7 +464,7 @@ | |||
| ELSE | |||
| K = K + 1 | |||
| W( K ) = Z( JLAM ) | |||
| DLAMDA( K ) = D( JLAM ) | |||
| DLAMBDA( K ) = D( JLAM ) | |||
| INDXP( K ) = JLAM | |||
| JLAM = J | |||
| END IF | |||
| @@ -476,26 +476,26 @@ | |||
| * | |||
| K = K + 1 | |||
| W( K ) = Z( JLAM ) | |||
| DLAMDA( K ) = D( JLAM ) | |||
| DLAMBDA( K ) = D( JLAM ) | |||
| INDXP( K ) = JLAM | |||
| * | |||
| 110 CONTINUE | |||
| * | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMDA | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMBDA | |||
| * and Q2 respectively. The eigenvalues/vectors which were not | |||
| * deflated go into the first K slots of DLAMDA and Q2 respectively, | |||
| * deflated go into the first K slots of DLAMBDA and Q2 respectively, | |||
| * while those which were deflated go into the last N - K slots. | |||
| * | |||
| IF( ICOMPQ.EQ.0 ) THEN | |||
| DO 120 J = 1, N | |||
| JP = INDXP( J ) | |||
| DLAMDA( J ) = D( JP ) | |||
| DLAMBDA( J ) = D( JP ) | |||
| PERM( J ) = INDXQ( INDX( JP ) ) | |||
| 120 CONTINUE | |||
| ELSE | |||
| DO 130 J = 1, N | |||
| JP = INDXP( J ) | |||
| DLAMDA( J ) = D( JP ) | |||
| DLAMBDA( J ) = D( JP ) | |||
| PERM( J ) = INDXQ( INDX( JP ) ) | |||
| CALL SCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) | |||
| 130 CONTINUE | |||
| @@ -506,9 +506,9 @@ | |||
| * | |||
| IF( K.LT.N ) THEN | |||
| IF( ICOMPQ.EQ.0 ) THEN | |||
| CALL SCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL SCOPY( N-K, DLAMBDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| ELSE | |||
| CALL SCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL SCOPY( N-K, DLAMBDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL SLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, | |||
| $ Q( 1, K+1 ), LDQ ) | |||
| END IF | |||
| @@ -18,15 +18,15 @@ | |||
| * Definition: | |||
| * =========== | |||
| * | |||
| * SUBROUTINE SLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, | |||
| * S, LDS, INFO ) | |||
| * SUBROUTINE SLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMBDA, | |||
| * W, S, LDS, INFO ) | |||
| * | |||
| * .. Scalar Arguments .. | |||
| * INTEGER INFO, K, KSTART, KSTOP, LDQ, LDS, N | |||
| * REAL RHO | |||
| * .. | |||
| * .. Array Arguments .. | |||
| * REAL D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), | |||
| * REAL D( * ), DLAMBDA( * ), Q( LDQ, * ), S( LDS, * ), | |||
| * $ W( * ) | |||
| * .. | |||
| * | |||
| @@ -96,9 +96,9 @@ | |||
| *> RHO >= 0 required. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in] DLAMDA | |||
| *> \param[in] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is REAL array, dimension (K) | |||
| *> DLAMBDA is REAL array, dimension (K) | |||
| *> The first K elements of this array contain the old roots | |||
| *> of the deflated updating problem. These are the poles | |||
| *> of the secular equation. | |||
| @@ -151,8 +151,8 @@ | |||
| *> at Berkeley, USA | |||
| * | |||
| * ===================================================================== | |||
| SUBROUTINE SLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMDA, W, | |||
| $ S, LDS, INFO ) | |||
| SUBROUTINE SLAED9( K, KSTART, KSTOP, N, D, Q, LDQ, RHO, DLAMBDA, | |||
| $ W, S, LDS, INFO ) | |||
| * | |||
| * -- LAPACK computational routine -- | |||
| * -- LAPACK is a software package provided by Univ. of Tennessee, -- | |||
| @@ -163,7 +163,7 @@ | |||
| REAL RHO | |||
| * .. | |||
| * .. Array Arguments .. | |||
| REAL D( * ), DLAMDA( * ), Q( LDQ, * ), S( LDS, * ), | |||
| REAL D( * ), DLAMBDA( * ), Q( LDQ, * ), S( LDS, * ), | |||
| $ W( * ) | |||
| * .. | |||
| * | |||
| @@ -174,8 +174,8 @@ | |||
| REAL TEMP | |||
| * .. | |||
| * .. External Functions .. | |||
| REAL SLAMC3, SNRM2 | |||
| EXTERNAL SLAMC3, SNRM2 | |||
| REAL SNRM2 | |||
| EXTERNAL SNRM2 | |||
| * .. | |||
| * .. External Subroutines .. | |||
| EXTERNAL SCOPY, SLAED4, XERBLA | |||
| @@ -212,30 +212,9 @@ | |||
| * | |||
| IF( K.EQ.0 ) | |||
| $ RETURN | |||
| * | |||
| * Modify values DLAMDA(i) to make sure all DLAMDA(i)-DLAMDA(j) can | |||
| * be computed with high relative accuracy (barring over/underflow). | |||
| * This is a problem on machines without a guard digit in | |||
| * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). | |||
| * The following code replaces DLAMDA(I) by 2*DLAMDA(I)-DLAMDA(I), | |||
| * which on any of these machines zeros out the bottommost | |||
| * bit of DLAMDA(I) if it is 1; this makes the subsequent | |||
| * subtractions DLAMDA(I)-DLAMDA(J) unproblematic when cancellation | |||
| * occurs. On binary machines with a guard digit (almost all | |||
| * machines) it does not change DLAMDA(I) at all. On hexadecimal | |||
| * and decimal machines with a guard digit, it slightly | |||
| * changes the bottommost bits of DLAMDA(I). It does not account | |||
| * for hexadecimal or decimal machines without guard digits | |||
| * (we know of none). We use a subroutine call to compute | |||
| * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating | |||
| * this code. | |||
| * | |||
| DO 10 I = 1, N | |||
| DLAMDA( I ) = SLAMC3( DLAMDA( I ), DLAMDA( I ) ) - DLAMDA( I ) | |||
| 10 CONTINUE | |||
| * | |||
| DO 20 J = KSTART, KSTOP | |||
| CALL SLAED4( K, J, DLAMDA, W, Q( 1, J ), RHO, D( J ), INFO ) | |||
| CALL SLAED4( K, J, DLAMBDA, W, Q( 1, J ), RHO, D( J ), INFO ) | |||
| * | |||
| * If the zero finder fails, the computation is terminated. | |||
| * | |||
| @@ -261,10 +240,10 @@ | |||
| CALL SCOPY( K, Q, LDQ+1, W, 1 ) | |||
| DO 70 J = 1, K | |||
| DO 50 I = 1, J - 1 | |||
| W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) | |||
| W( I ) = W( I )*( Q( I, J )/( DLAMBDA( I )-DLAMBDA( J ) ) ) | |||
| 50 CONTINUE | |||
| DO 60 I = J + 1, K | |||
| W( I ) = W( I )*( Q( I, J ) / ( DLAMDA( I )-DLAMDA( J ) ) ) | |||
| W( I ) = W( I )*( Q( I, J )/( DLAMBDA( I )-DLAMBDA( J ) ) ) | |||
| 60 CONTINUE | |||
| 70 CONTINUE | |||
| DO 80 I = 1, K | |||
| @@ -389,6 +389,11 @@ | |||
| $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN | |||
| WORK( I ) = ZERO | |||
| ELSE | |||
| * | |||
| * Use calls to the subroutine SLAMC3 to enforce the | |||
| * parentheses (x+y)+z. The goal is to prevent | |||
| * optimizing compilers from doing x+(y+z). | |||
| * | |||
| WORK( I ) = POLES( I, 2 )*Z( I ) / | |||
| $ ( SLAMC3( POLES( I, 2 ), DSIGJ )- | |||
| $ DIFLJ ) / ( POLES( I, 2 )+DJ ) | |||
| @@ -440,6 +445,11 @@ | |||
| IF( Z( J ).EQ.ZERO ) THEN | |||
| WORK( I ) = ZERO | |||
| ELSE | |||
| * | |||
| * Use calls to the subroutine SLAMC3 to enforce the | |||
| * parentheses (x+y)+z. The goal is to prevent | |||
| * optimizing compilers from doing x+(y+z). | |||
| * | |||
| WORK( I ) = Z( J ) / ( SLAMC3( DSIGJ, -POLES( I+1, | |||
| $ 2 ) )-DIFR( I, 1 ) ) / | |||
| $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) | |||
| @@ -47,12 +47,6 @@ | |||
| *> problem; in this case a minimum norm solution is returned. | |||
| *> The actual singular values are returned in D in ascending order. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -93,9 +93,7 @@ | |||
| *> infinite. | |||
| *> | |||
| *> Overflow will not occur unless the largest singular value itself | |||
| *> overflows, or is within a few ulps of overflow. (On machines with | |||
| *> partial overflow, like the Cray, overflow may occur if the largest | |||
| *> singular value is within a factor of 2 of overflow.) | |||
| *> overflows, or is within a few ulps of overflow. | |||
| *> | |||
| *> Underflow is harmless if underflow is gradual. Otherwise, results | |||
| *> may correspond to a matrix modified by perturbations of size near | |||
| @@ -44,13 +44,6 @@ | |||
| *> appropriate calls to SLASD4 and then updates the singular | |||
| *> vectors by matrix multiplication. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> | |||
| *> SLASD3 is called from SLASD1. | |||
| *> \endverbatim | |||
| * | |||
| @@ -103,7 +96,7 @@ | |||
| *> The leading dimension of the array Q. LDQ >= K. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in,out] DSIGMA | |||
| *> \param[in] DSIGMA | |||
| *> \verbatim | |||
| *> DSIGMA is REAL array, dimension(K) | |||
| *> The first K elements of this array contain the old roots | |||
| @@ -249,8 +242,8 @@ | |||
| REAL RHO, TEMP | |||
| * .. | |||
| * .. External Functions .. | |||
| REAL SLAMC3, SNRM2 | |||
| EXTERNAL SLAMC3, SNRM2 | |||
| REAL SNRM2 | |||
| EXTERNAL SNRM2 | |||
| * .. | |||
| * .. External Subroutines .. | |||
| EXTERNAL SCOPY, SGEMM, SLACPY, SLASCL, SLASD4, XERBLA | |||
| @@ -310,27 +303,6 @@ | |||
| RETURN | |||
| END IF | |||
| * | |||
| * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can | |||
| * be computed with high relative accuracy (barring over/underflow). | |||
| * This is a problem on machines without a guard digit in | |||
| * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). | |||
| * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), | |||
| * which on any of these machines zeros out the bottommost | |||
| * bit of DSIGMA(I) if it is 1; this makes the subsequent | |||
| * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation | |||
| * occurs. On binary machines with a guard digit (almost all | |||
| * machines) it does not change DSIGMA(I) at all. On hexadecimal | |||
| * and decimal machines with a guard digit, it slightly | |||
| * changes the bottommost bits of DSIGMA(I). It does not account | |||
| * for hexadecimal or decimal machines without guard digits | |||
| * (we know of none). We use a subroutine call to compute | |||
| * 2*DSIGMA(I) to prevent optimizing compilers from eliminating | |||
| * this code. | |||
| * | |||
| DO 20 I = 1, K | |||
| DSIGMA( I ) = SLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) | |||
| 20 CONTINUE | |||
| * | |||
| * Keep a copy of Z. | |||
| * | |||
| CALL SCOPY( K, Z, 1, Q, 1 ) | |||
| @@ -121,14 +121,12 @@ | |||
| *> The leading dimension of DIFR, must be at least K. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[in,out] DSIGMA | |||
| *> \param[in] DSIGMA | |||
| *> \verbatim | |||
| *> DSIGMA is REAL array, dimension ( K ) | |||
| *> On entry, the first K elements of this array contain the old | |||
| *> roots of the deflated updating problem. These are the poles | |||
| *> of the secular equation. | |||
| *> On exit, the elements of DSIGMA may be very slightly altered | |||
| *> in value. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[out] WORK | |||
| @@ -227,27 +225,6 @@ | |||
| RETURN | |||
| END IF | |||
| * | |||
| * Modify values DSIGMA(i) to make sure all DSIGMA(i)-DSIGMA(j) can | |||
| * be computed with high relative accuracy (barring over/underflow). | |||
| * This is a problem on machines without a guard digit in | |||
| * add/subtract (Cray XMP, Cray YMP, Cray C 90 and Cray 2). | |||
| * The following code replaces DSIGMA(I) by 2*DSIGMA(I)-DSIGMA(I), | |||
| * which on any of these machines zeros out the bottommost | |||
| * bit of DSIGMA(I) if it is 1; this makes the subsequent | |||
| * subtractions DSIGMA(I)-DSIGMA(J) unproblematic when cancellation | |||
| * occurs. On binary machines with a guard digit (almost all | |||
| * machines) it does not change DSIGMA(I) at all. On hexadecimal | |||
| * and decimal machines with a guard digit, it slightly | |||
| * changes the bottommost bits of DSIGMA(I). It does not account | |||
| * for hexadecimal or decimal machines without guard digits | |||
| * (we know of none). We use a subroutine call to compute | |||
| * 2*DLAMBDA(I) to prevent optimizing compilers from eliminating | |||
| * this code. | |||
| * | |||
| DO 10 I = 1, K | |||
| DSIGMA( I ) = SLAMC3( DSIGMA( I ), DSIGMA( I ) ) - DSIGMA( I ) | |||
| 10 CONTINUE | |||
| * | |||
| * Book keeping. | |||
| * | |||
| IWK1 = 1 | |||
| @@ -312,6 +289,11 @@ | |||
| DSIGJP = -DSIGMA( J+1 ) | |||
| END IF | |||
| WORK( J ) = -Z( J ) / DIFLJ / ( DSIGMA( J )+DJ ) | |||
| * | |||
| * Use calls to the subroutine SLAMC3 to enforce the parentheses | |||
| * (x+y)+z. The goal is to prevent optimizing compilers | |||
| * from doing x+(y+z). | |||
| * | |||
| DO 60 I = 1, J - 1 | |||
| WORK( I ) = Z( I ) / ( SLAMC3( DSIGMA( I ), DSIGJ )-DIFLJ ) | |||
| $ / ( DSIGMA( I )+DJ ) | |||
| @@ -124,9 +124,7 @@ | |||
| *> infinite. | |||
| *> | |||
| *> Overflow will not occur unless the largest singular value itself | |||
| *> overflows or is within a few ulps of overflow. (On machines with | |||
| *> partial overflow, like the Cray, overflow may occur if the largest | |||
| *> singular value is within a factor of 2 of overflow.) | |||
| *> overflows or is within a few ulps of overflow. | |||
| *> | |||
| *> Underflow is harmless if underflow is gradual. Otherwise, results | |||
| *> may correspond to a matrix modified by perturbations of size near | |||
| @@ -40,12 +40,6 @@ | |||
| *> a real symmetric band matrix A. If eigenvectors are desired, it uses | |||
| *> a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -45,12 +45,6 @@ | |||
| *> the reduction to tridiagonal. If eigenvectors are desired, it uses | |||
| *> a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -43,12 +43,6 @@ | |||
| *> banded, and B is also positive definite. If eigenvectors are | |||
| *> desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -40,12 +40,6 @@ | |||
| *> of a real symmetric matrix A in packed storage. If eigenvectors are | |||
| *> desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -44,12 +44,6 @@ | |||
| *> positive definite. | |||
| *> If eigenvectors are desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -42,12 +42,6 @@ | |||
| *> found if SSYTRD or SSPTRD or SSBTRD has been used to reduce this | |||
| *> matrix to tridiagonal form. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. See SLAED3 for details. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -40,12 +40,6 @@ | |||
| *> real symmetric tridiagonal matrix. If eigenvectors are desired, it | |||
| *> uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -40,13 +40,6 @@ | |||
| *> real symmetric matrix A. If eigenvectors are desired, it uses a | |||
| *> divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> | |||
| *> Because of large use of BLAS of level 3, SSYEVD needs N**2 more | |||
| *> workspace than SSYEVX. | |||
| *> \endverbatim | |||
| @@ -45,12 +45,6 @@ | |||
| *> the reduction to tridiagonal. If eigenvectors are desired, it uses a | |||
| *> divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -42,12 +42,6 @@ | |||
| *> B are assumed to be symmetric and B is also positive definite. | |||
| *> If eigenvectors are desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -60,12 +60,6 @@ | |||
| *> singular values which are less than RCOND times the largest singular | |||
| *> value. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -53,12 +53,6 @@ | |||
| *> | |||
| *> Note that the routine returns VT = V**H, not V. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -41,12 +41,6 @@ | |||
| *> a complex Hermitian band matrix A. If eigenvectors are desired, it | |||
| *> uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -47,12 +47,6 @@ | |||
| *> the reduction to tridiagonal. If eigenvectors are desired, it | |||
| *> uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -46,12 +46,6 @@ | |||
| *> and banded, and B is also positive definite. If eigenvectors are | |||
| *> desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -41,12 +41,6 @@ | |||
| *> complex Hermitian matrix A. If eigenvectors are desired, it uses a | |||
| *> divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -46,12 +46,6 @@ | |||
| *> the reduction to tridiagonal. If eigenvectors are desired, it uses a | |||
| *> divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -43,12 +43,6 @@ | |||
| *> B are assumed to be Hermitian and B is also positive definite. | |||
| *> If eigenvectors are desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -41,12 +41,6 @@ | |||
| *> a complex Hermitian matrix A in packed storage. If eigenvectors are | |||
| *> desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -44,12 +44,6 @@ | |||
| *> positive definite. | |||
| *> If eigenvectors are desired, it uses a divide and conquer algorithm. | |||
| *> | |||
| *> The divide and conquer algorithm makes very mild assumptions about | |||
| *> floating point arithmetic. It will work on machines with a guard | |||
| *> digit in add/subtract, or on those binary machines without guard | |||
| *> digits which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or | |||
| *> Cray-2. It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -18,7 +18,7 @@ | |||
| * Definition: | |||
| * =========== | |||
| * | |||
| * SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA, | |||
| * SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMBDA, | |||
| * Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR, | |||
| * GIVCOL, GIVNUM, INFO ) | |||
| * | |||
| @@ -29,7 +29,7 @@ | |||
| * .. Array Arguments .. | |||
| * INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), | |||
| * $ INDXQ( * ), PERM( * ) | |||
| * DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ), | |||
| * DOUBLE PRECISION D( * ), DLAMBDA( * ), GIVNUM( 2, * ), W( * ), | |||
| * $ Z( * ) | |||
| * COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * ) | |||
| * .. | |||
| @@ -122,9 +122,9 @@ | |||
| *> destroyed during the updating process. | |||
| *> \endverbatim | |||
| *> | |||
| *> \param[out] DLAMDA | |||
| *> \param[out] DLAMBDA | |||
| *> \verbatim | |||
| *> DLAMDA is DOUBLE PRECISION array, dimension (N) | |||
| *> DLAMBDA is DOUBLE PRECISION array, dimension (N) | |||
| *> Contains a copy of the first K eigenvalues which will be used | |||
| *> by DLAED3 to form the secular equation. | |||
| *> \endverbatim | |||
| @@ -222,7 +222,7 @@ | |||
| *> \ingroup complex16OTHERcomputational | |||
| * | |||
| * ===================================================================== | |||
| SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMDA, | |||
| SUBROUTINE ZLAED8( K, N, QSIZ, Q, LDQ, D, RHO, CUTPNT, Z, DLAMBDA, | |||
| $ Q2, LDQ2, W, INDXP, INDX, INDXQ, PERM, GIVPTR, | |||
| $ GIVCOL, GIVNUM, INFO ) | |||
| * | |||
| @@ -237,7 +237,7 @@ | |||
| * .. Array Arguments .. | |||
| INTEGER GIVCOL( 2, * ), INDX( * ), INDXP( * ), | |||
| $ INDXQ( * ), PERM( * ) | |||
| DOUBLE PRECISION D( * ), DLAMDA( * ), GIVNUM( 2, * ), W( * ), | |||
| DOUBLE PRECISION D( * ), DLAMBDA( * ), GIVNUM( 2, * ), W( * ), | |||
| $ Z( * ) | |||
| COMPLEX*16 Q( LDQ, * ), Q2( LDQ2, * ) | |||
| * .. | |||
| @@ -322,14 +322,14 @@ | |||
| INDXQ( I ) = INDXQ( I ) + CUTPNT | |||
| 20 CONTINUE | |||
| DO 30 I = 1, N | |||
| DLAMDA( I ) = D( INDXQ( I ) ) | |||
| DLAMBDA( I ) = D( INDXQ( I ) ) | |||
| W( I ) = Z( INDXQ( I ) ) | |||
| 30 CONTINUE | |||
| I = 1 | |||
| J = CUTPNT + 1 | |||
| CALL DLAMRG( N1, N2, DLAMDA, 1, 1, INDX ) | |||
| CALL DLAMRG( N1, N2, DLAMBDA, 1, 1, INDX ) | |||
| DO 40 I = 1, N | |||
| D( I ) = DLAMDA( INDX( I ) ) | |||
| D( I ) = DLAMBDA( INDX( I ) ) | |||
| Z( I ) = W( INDX( I ) ) | |||
| 40 CONTINUE | |||
| * | |||
| @@ -438,7 +438,7 @@ | |||
| ELSE | |||
| K = K + 1 | |||
| W( K ) = Z( JLAM ) | |||
| DLAMDA( K ) = D( JLAM ) | |||
| DLAMBDA( K ) = D( JLAM ) | |||
| INDXP( K ) = JLAM | |||
| JLAM = J | |||
| END IF | |||
| @@ -450,19 +450,19 @@ | |||
| * | |||
| K = K + 1 | |||
| W( K ) = Z( JLAM ) | |||
| DLAMDA( K ) = D( JLAM ) | |||
| DLAMBDA( K ) = D( JLAM ) | |||
| INDXP( K ) = JLAM | |||
| * | |||
| 100 CONTINUE | |||
| * | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMDA | |||
| * Sort the eigenvalues and corresponding eigenvectors into DLAMBDA | |||
| * and Q2 respectively. The eigenvalues/vectors which were not | |||
| * deflated go into the first K slots of DLAMDA and Q2 respectively, | |||
| * deflated go into the first K slots of DLAMBDA and Q2 respectively, | |||
| * while those which were deflated go into the last N - K slots. | |||
| * | |||
| DO 110 J = 1, N | |||
| JP = INDXP( J ) | |||
| DLAMDA( J ) = D( JP ) | |||
| DLAMBDA( J ) = D( JP ) | |||
| PERM( J ) = INDXQ( INDX( JP ) ) | |||
| CALL ZCOPY( QSIZ, Q( 1, PERM( J ) ), 1, Q2( 1, J ), 1 ) | |||
| 110 CONTINUE | |||
| @@ -471,7 +471,7 @@ | |||
| * into the last N - K slots of D and Q respectively. | |||
| * | |||
| IF( K.LT.N ) THEN | |||
| CALL DCOPY( N-K, DLAMDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL DCOPY( N-K, DLAMBDA( K+1 ), 1, D( K+1 ), 1 ) | |||
| CALL ZLACPY( 'A', QSIZ, N-K, Q2( 1, K+1 ), LDQ2, Q( 1, K+1 ), | |||
| $ LDQ ) | |||
| END IF | |||
| @@ -392,6 +392,11 @@ | |||
| $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN | |||
| RWORK( I ) = ZERO | |||
| ELSE | |||
| * | |||
| * Use calls to the subroutine DLAMC3 to enforce the | |||
| * parentheses (x+y)+z. The goal is to prevent | |||
| * optimizing compilers from doing x+(y+z). | |||
| * | |||
| RWORK( I ) = POLES( I, 2 )*Z( I ) / | |||
| $ ( DLAMC3( POLES( I, 2 ), DSIGJ )- | |||
| $ DIFLJ ) / ( POLES( I, 2 )+DJ ) | |||
| @@ -470,6 +475,11 @@ | |||
| IF( Z( J ).EQ.ZERO ) THEN | |||
| RWORK( I ) = ZERO | |||
| ELSE | |||
| * | |||
| * Use calls to the subroutine DLAMC3 to enforce the | |||
| * parentheses (x+y)+z. The goal is to prevent | |||
| * optimizing compilers from doing x+(y+z). | |||
| * | |||
| RWORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1, | |||
| $ 2 ) )-DIFR( I, 1 ) ) / | |||
| $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) | |||
| @@ -48,12 +48,6 @@ | |||
| *> problem; in this case a minimum norm solution is returned. | |||
| *> The actual singular values are returned in D in ascending order. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||
| @@ -43,12 +43,6 @@ | |||
| *> be found if ZHETRD or ZHPTRD or ZHBTRD has been used to reduce this | |||
| *> matrix to tridiagonal form. | |||
| *> | |||
| *> This code makes very mild assumptions about floating point | |||
| *> arithmetic. It will work on machines with a guard digit in | |||
| *> add/subtract, or on those binary machines without guard digits | |||
| *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2. | |||
| *> It could conceivably fail on hexadecimal or decimal machines | |||
| *> without guard digits, but we know of none. See DLAED3 for details. | |||
| *> \endverbatim | |||
| * | |||
| * Arguments: | |||