dlasq2.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 473 行 · 第 1/2 页
HTML
473 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dlasq2.f</title>
<meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!--
body { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.default a { background: rgb(255, 255, 255); color: rgb(0, 0, 0); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.string { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.string a { color: rgb(188, 143, 143); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
span.comment { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: none; }
span.comment a { color: rgb(178, 34, 34); background: rgb(255, 255, 255); font-style: normal; font-weight: 500; font-stretch: normal; font-family: adobe-courier; font-size: 11pt; text-decoration: underline; }
--></style>
</head>
<body>
<pre>
SUBROUTINE <a name="DLASQ2.1"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>( N, Z, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK routine (version 3.1) --
</span><span class="comment">*</span><span class="comment"> Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment"> November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Modified to call <a name="DLAZQ3.7"></a><a href="dlazq3.f.html#DLAZQ3.1">DLAZQ3</a> in place of <a name="DLASQ3.7"></a><a href="dlasq3.f.html#DLASQ3.1">DLASQ3</a>, 13 Feb 03, SJH.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Scalar Arguments ..
</span> INTEGER INFO, N
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> DOUBLE PRECISION Z( * )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Purpose
</span><span class="comment">*</span><span class="comment"> =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> <a name="DLASQ2.19"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a> computes all the eigenvalues of the symmetric positive
</span><span class="comment">*</span><span class="comment"> definite tridiagonal matrix associated with the qd array Z to high
</span><span class="comment">*</span><span class="comment"> relative accuracy are computed to high relative accuracy, in the
</span><span class="comment">*</span><span class="comment"> absence of denormalization, underflow and overflow.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> To see the relation of Z to the tridiagonal matrix, let L be a
</span><span class="comment">*</span><span class="comment"> unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and
</span><span class="comment">*</span><span class="comment"> let U be an upper bidiagonal matrix with 1's above and diagonal
</span><span class="comment">*</span><span class="comment"> Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
</span><span class="comment">*</span><span class="comment"> symmetric tridiagonal to which it is similar.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Note : <a name="DLASQ2.30"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a> defines a logical variable, IEEE, which is true
</span><span class="comment">*</span><span class="comment"> on machines which follow ieee-754 floating-point standard in their
</span><span class="comment">*</span><span class="comment"> handling of infinities and NaNs, and false otherwise. This variable
</span><span class="comment">*</span><span class="comment"> is passed to <a name="DLAZQ3.33"></a><a href="dlazq3.f.html#DLAZQ3.1">DLAZQ3</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Arguments
</span><span class="comment">*</span><span class="comment"> =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> N (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The number of rows and columns in the matrix. N >= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Z (workspace) DOUBLE PRECISION array, dimension ( 4*N )
</span><span class="comment">*</span><span class="comment"> On entry Z holds the qd array. On exit, entries 1 to N hold
</span><span class="comment">*</span><span class="comment"> the eigenvalues in decreasing order, Z( 2*N+1 ) holds the
</span><span class="comment">*</span><span class="comment"> trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If
</span><span class="comment">*</span><span class="comment"> N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )
</span><span class="comment">*</span><span class="comment"> holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of
</span><span class="comment">*</span><span class="comment"> shifts that failed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> INFO (output) INTEGER
</span><span class="comment">*</span><span class="comment"> = 0: successful exit
</span><span class="comment">*</span><span class="comment"> < 0: if the i-th argument is a scalar and had an illegal
</span><span class="comment">*</span><span class="comment"> value, then INFO = -i, if the i-th argument is an
</span><span class="comment">*</span><span class="comment"> array and the j-entry had an illegal value, then
</span><span class="comment">*</span><span class="comment"> INFO = -(i*100+j)
</span><span class="comment">*</span><span class="comment"> > 0: the algorithm failed
</span><span class="comment">*</span><span class="comment"> = 1, a split was marked by a positive value in E
</span><span class="comment">*</span><span class="comment"> = 2, current block of Z not diagonalized after 30*N
</span><span class="comment">*</span><span class="comment"> iterations (in inner while loop)
</span><span class="comment">*</span><span class="comment"> = 3, termination criterion of outer while loop not met
</span><span class="comment">*</span><span class="comment"> (program created more than N unreduced blocks)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Further Details
</span><span class="comment">*</span><span class="comment"> ===============
</span><span class="comment">*</span><span class="comment"> Local Variables: I0:N0 defines a current unreduced segment of Z.
</span><span class="comment">*</span><span class="comment"> The shifts are accumulated in SIGMA. Iteration count is in ITER.
</span><span class="comment">*</span><span class="comment"> Ping-pong is controlled by PP (alternates between 0 and 1).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> .. Parameters ..
</span> DOUBLE PRECISION CBIAS
PARAMETER ( CBIAS = 1.50D0 )
DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD
PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0,
$ TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL IEEE
INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K,
$ N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE
DOUBLE PRECISION D, DESIG, DMIN, DMIN1, DMIN2, DN, DN1, DN2, E,
$ EMAX, EMIN, EPS, OLDEMN, QMAX, QMIN, S, SAFMIN,
$ SIGMA, T, TAU, TEMP, TOL, TOL2, TRACE, ZMAX
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="DLAZQ3.86"></a><a href="dlazq3.f.html#DLAZQ3.1">DLAZQ3</a>, <a name="DLASRT.86"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>, <a name="XERBLA.86"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> INTEGER <a name="ILAENV.89"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
DOUBLE PRECISION <a name="DLAMCH.90"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
EXTERNAL <a name="DLAMCH.91"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>, <a name="ILAENV.91"></a><a href="hfy-index.html#ILAENV">ILAENV</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC ABS, DBLE, MAX, MIN, SQRT
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Test the input arguments.
</span><span class="comment">*</span><span class="comment"> (in case <a name="DLASQ2.99"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a> is not called by <a name="DLASQ1.99"></a><a href="dlasq1.f.html#DLASQ1.1">DLASQ1</a>)
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
EPS = <a name="DLAMCH.102"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
SAFMIN = <a name="DLAMCH.103"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Safe minimum'</span> )
TOL = EPS*HUNDRD
TOL2 = TOL**2
<span class="comment">*</span><span class="comment">
</span> IF( N.LT.0 ) THEN
INFO = -1
CALL <a name="XERBLA.109"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLASQ2.109"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>'</span>, 1 )
RETURN
ELSE IF( N.EQ.0 ) THEN
RETURN
ELSE IF( N.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> 1-by-1 case.
</span><span class="comment">*</span><span class="comment">
</span> IF( Z( 1 ).LT.ZERO ) THEN
INFO = -201
CALL <a name="XERBLA.119"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLASQ2.119"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>'</span>, 2 )
END IF
RETURN
ELSE IF( N.EQ.2 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> 2-by-2 case.
</span><span class="comment">*</span><span class="comment">
</span> IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN
INFO = -2
CALL <a name="XERBLA.128"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLASQ2.128"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>'</span>, 2 )
RETURN
ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN
D = Z( 3 )
Z( 3 ) = Z( 1 )
Z( 1 ) = D
END IF
Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 )
IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN
T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) )
S = Z( 3 )*( Z( 2 ) / T )
IF( S.LE.T ) THEN
S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) )
ELSE
S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) )
END IF
T = Z( 1 ) + ( S+Z( 2 ) )
Z( 3 ) = Z( 3 )*( Z( 1 ) / T )
Z( 1 ) = T
END IF
Z( 2 ) = Z( 3 )
Z( 6 ) = Z( 2 ) + Z( 1 )
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Check for negative data and compute sums of q's and e's.
</span><span class="comment">*</span><span class="comment">
</span> Z( 2*N ) = ZERO
EMIN = Z( 2 )
QMAX = ZERO
ZMAX = ZERO
D = ZERO
E = ZERO
<span class="comment">*</span><span class="comment">
</span> DO 10 K = 1, 2*( N-1 ), 2
IF( Z( K ).LT.ZERO ) THEN
INFO = -( 200+K )
CALL <a name="XERBLA.165"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLASQ2.165"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>'</span>, 2 )
RETURN
ELSE IF( Z( K+1 ).LT.ZERO ) THEN
INFO = -( 200+K+1 )
CALL <a name="XERBLA.169"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLASQ2.169"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>'</span>, 2 )
RETURN
END IF
D = D + Z( K )
E = E + Z( K+1 )
QMAX = MAX( QMAX, Z( K ) )
EMIN = MIN( EMIN, Z( K+1 ) )
ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) )
10 CONTINUE
IF( Z( 2*N-1 ).LT.ZERO ) THEN
INFO = -( 200+2*N-1 )
CALL <a name="XERBLA.180"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DLASQ2.180"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>'</span>, 2 )
RETURN
END IF
D = D + Z( 2*N-1 )
QMAX = MAX( QMAX, Z( 2*N-1 ) )
ZMAX = MAX( QMAX, ZMAX )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Check for diagonality.
</span><span class="comment">*</span><span class="comment">
</span> IF( E.EQ.ZERO ) THEN
DO 20 K = 2, N
Z( K ) = Z( 2*K-1 )
20 CONTINUE
CALL <a name="DLASRT.193"></a><a href="dlasrt.f.html#DLASRT.1">DLASRT</a>( <span class="string">'D'</span>, N, Z, IINFO )
Z( 2*N-1 ) = D
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span> TRACE = D + E
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Check for zero data.
</span><span class="comment">*</span><span class="comment">
</span> IF( TRACE.EQ.ZERO ) THEN
Z( 2*N-1 ) = ZERO
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Check whether the machine is IEEE conformable.
</span><span class="comment">*</span><span class="comment">
</span> IEEE = <a name="ILAENV.209"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 10, <span class="string">'<a name="DLASQ2.209"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>'</span>, <span class="string">'N'</span>, 1, 2, 3, 4 ).EQ.1 .AND.
$ <a name="ILAENV.210"></a><a href="hfy-index.html#ILAENV">ILAENV</a>( 11, <span class="string">'<a name="DLASQ2.210"></a><a href="dlasq2.f.html#DLASQ2.1">DLASQ2</a>'</span>, <span class="string">'N'</span>, 1, 2, 3, 4 ).EQ.1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
</span><span class="comment">*</span><span class="comment">
</span> DO 30 K = 2*N, 2, -2
Z( 2*K ) = ZERO
Z( 2*K-1 ) = Z( K )
Z( 2*K-2 ) = ZERO
Z( 2*K-3 ) = Z( K-1 )
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?