dlasd4.f.html
来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 915 行 · 第 1/4 页
HTML
915 行
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<title>dlasd4.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="DLASD4.1"></a><a href="dlasd4.f.html#DLASD4.1">DLASD4</a>( N, I, D, Z, DELTA, RHO, SIGMA, WORK, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> -- LAPACK auxiliary 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"> .. Scalar Arguments ..
</span> INTEGER I, INFO, N
DOUBLE PRECISION RHO, SIGMA
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Array Arguments ..
</span> DOUBLE PRECISION D( * ), DELTA( * ), WORK( * ), 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"> This subroutine computes the square root of the I-th updated
</span><span class="comment">*</span><span class="comment"> eigenvalue of a positive symmetric rank-one modification to
</span><span class="comment">*</span><span class="comment"> a positive diagonal matrix whose entries are given as the squares
</span><span class="comment">*</span><span class="comment"> of the corresponding entries in the array d, and that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> 0 <= D(i) < D(j) for i < j
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> and that RHO > 0. This is arranged by the calling routine, and is
</span><span class="comment">*</span><span class="comment"> no loss in generality. The rank-one modified system is thus
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> diag( D ) * diag( D ) + RHO * Z * Z_transpose.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> where we assume the Euclidean norm of Z is 1.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The method consists of approximating the rational functions in the
</span><span class="comment">*</span><span class="comment"> secular equation by simpler interpolating rational functions.
</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 length of all arrays.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> I (input) INTEGER
</span><span class="comment">*</span><span class="comment"> The index of the eigenvalue to be computed. 1 <= I <= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> D (input) DOUBLE PRECISION array, dimension ( N )
</span><span class="comment">*</span><span class="comment"> The original eigenvalues. It is assumed that they are in
</span><span class="comment">*</span><span class="comment"> order, 0 <= D(I) < D(J) for I < J.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Z (input) DOUBLE PRECISION array, dimension ( N )
</span><span class="comment">*</span><span class="comment"> The components of the updating vector.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> DELTA (output) DOUBLE PRECISION array, dimension ( N )
</span><span class="comment">*</span><span class="comment"> If N .ne. 1, DELTA contains (D(j) - sigma_I) in its j-th
</span><span class="comment">*</span><span class="comment"> component. If N = 1, then DELTA(1) = 1. The vector DELTA
</span><span class="comment">*</span><span class="comment"> contains the information necessary to construct the
</span><span class="comment">*</span><span class="comment"> (singular) eigenvectors.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> RHO (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> The scalar in the symmetric updating formula.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> SIGMA (output) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment"> The computed sigma_I, the I-th updated eigenvalue.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> WORK (workspace) DOUBLE PRECISION array, dimension ( N )
</span><span class="comment">*</span><span class="comment"> If N .ne. 1, WORK contains (D(j) + sigma_I) in its j-th
</span><span class="comment">*</span><span class="comment"> component. If N = 1, then WORK( 1 ) = 1.
</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 INFO = 1, the updating process failed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Internal Parameters
</span><span class="comment">*</span><span class="comment"> ===================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Logical variable ORGATI (origin-at-i?) is used for distinguishing
</span><span class="comment">*</span><span class="comment"> whether D(i) or D(i+1) is treated as the origin.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> ORGATI = .true. origin at i
</span><span class="comment">*</span><span class="comment"> ORGATI = .false. origin at i+1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Logical variable SWTCH3 (switch-for-3-poles?) is for noting
</span><span class="comment">*</span><span class="comment"> if we are working with THREE poles!
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> MAXIT is the maximum number of iterations allowed for each
</span><span class="comment">*</span><span class="comment"> eigenvalue.
</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">
</span><span class="comment">*</span><span class="comment"> Based on contributions by
</span><span class="comment">*</span><span class="comment"> Ren-Cang Li, Computer Science Division, University of California
</span><span class="comment">*</span><span class="comment"> at Berkeley, USA
</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> INTEGER MAXIT
PARAMETER ( MAXIT = 20 )
DOUBLE PRECISION ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
$ THREE = 3.0D+0, FOUR = 4.0D+0, EIGHT = 8.0D+0,
$ TEN = 10.0D+0 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Scalars ..
</span> LOGICAL ORGATI, SWTCH, SWTCH3
INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
DOUBLE PRECISION A, B, C, DELSQ, DELSQ2, DPHI, DPSI, DTIIM,
$ DTIIP, DTIPSQ, DTISQ, DTNSQ, DTNSQ1, DW, EPS,
$ ERRETM, ETA, PHI, PREW, PSI, RHOINV, SG2LB,
$ SG2UB, TAU, TEMP, TEMP1, TEMP2, W
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Local Arrays ..
</span> DOUBLE PRECISION DD( 3 ), ZZ( 3 )
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Subroutines ..
</span> EXTERNAL <a name="DLAED6.115"></a><a href="dlaed6.f.html#DLAED6.1">DLAED6</a>, <a name="DLASD5.115"></a><a href="dlasd5.f.html#DLASD5.1">DLASD5</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. External Functions ..
</span> DOUBLE PRECISION <a name="DLAMCH.118"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
EXTERNAL <a name="DLAMCH.119"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
<span class="comment">*</span><span class="comment"> ..
</span><span class="comment">*</span><span class="comment"> .. Intrinsic Functions ..
</span> INTRINSIC ABS, 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"> Since this routine is called in an inner loop, we do no argument
</span><span class="comment">*</span><span class="comment"> checking.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Quick return for N=1 and 2.
</span><span class="comment">*</span><span class="comment">
</span> INFO = 0
IF( N.EQ.1 ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Presumably, I=1 upon entry
</span><span class="comment">*</span><span class="comment">
</span> SIGMA = SQRT( D( 1 )*D( 1 )+RHO*Z( 1 )*Z( 1 ) )
DELTA( 1 ) = ONE
WORK( 1 ) = ONE
RETURN
END IF
IF( N.EQ.2 ) THEN
CALL <a name="DLASD5.142"></a><a href="dlasd5.f.html#DLASD5.1">DLASD5</a>( I, D, Z, DELTA, RHO, SIGMA, WORK )
RETURN
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Compute machine epsilon
</span><span class="comment">*</span><span class="comment">
</span> EPS = <a name="DLAMCH.148"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Epsilon'</span> )
RHOINV = ONE / RHO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The case I = N
</span><span class="comment">*</span><span class="comment">
</span> IF( I.EQ.N ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Initialize some basic variables
</span><span class="comment">*</span><span class="comment">
</span> II = N - 1
NITER = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> Calculate initial guess
</span><span class="comment">*</span><span class="comment">
</span> TEMP = RHO / TWO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> If ||Z||_2 is not one, then TEMP should be set to
</span><span class="comment">*</span><span class="comment"> RHO * ||Z||_2^2 / TWO
</span><span class="comment">*</span><span class="comment">
</span> TEMP1 = TEMP / ( D( N )+SQRT( D( N )*D( N )+TEMP ) )
DO 10 J = 1, N
WORK( J ) = D( J ) + D( N ) + TEMP1
DELTA( J ) = ( D( J )-D( N ) ) - TEMP1
10 CONTINUE
<span class="comment">*</span><span class="comment">
</span> PSI = ZERO
DO 20 J = 1, N - 2
PSI = PSI + Z( J )*Z( J ) / ( DELTA( J )*WORK( J ) )
20 CONTINUE
<span class="comment">*</span><span class="comment">
</span> C = RHOINV + PSI
W = C + Z( II )*Z( II ) / ( DELTA( II )*WORK( II ) ) +
$ Z( N )*Z( N ) / ( DELTA( N )*WORK( N ) )
<span class="comment">*</span><span class="comment">
</span> IF( W.LE.ZERO ) THEN
TEMP1 = SQRT( D( N )*D( N )+RHO )
TEMP = Z( N-1 )*Z( N-1 ) / ( ( D( N-1 )+TEMP1 )*
$ ( D( N )-D( N-1 )+RHO / ( D( N )+TEMP1 ) ) ) +
$ Z( N )*Z( N ) / RHO
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> The following TAU is to approximate
</span><span class="comment">*</span><span class="comment"> SIGMA_n^2 - D( N )*D( N )
</span><span class="comment">*</span><span class="comment">
</span> IF( C.LE.TEMP ) THEN
TAU = RHO
ELSE
DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
B = Z( N )*Z( N )*DELSQ
IF( A.LT.ZERO ) THEN
TAU = TWO*B / ( SQRT( A*A+FOUR*B*C )-A )
ELSE
TAU = ( A+SQRT( A*A+FOUR*B*C ) ) / ( TWO*C )
END IF
END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> It can be proved that
</span><span class="comment">*</span><span class="comment"> D(N)^2+RHO/2 <= SIGMA_n^2 < D(N)^2+TAU <= D(N)^2+RHO
</span><span class="comment">*</span><span class="comment">
</span> ELSE
DELSQ = ( D( N )-D( N-1 ) )*( D( N )+D( N-1 ) )
A = -C*DELSQ + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
B = Z( N )*Z( N )*DELSQ
⌨️ 快捷键说明
复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?