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