slaed4.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 869 行 · 第 1/3 页

HTML
869 行
字号
123
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>slaed4.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="SLAED4.1"></a><a href="slaed4.f.html#SLAED4.1">SLAED4</a>( N, I, D, Z, DELTA, RHO, DLAM, 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">     .. Scalar Arguments ..
</span>      INTEGER            I, INFO, N
      REAL               DLAM, RHO
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      REAL               D( * ), DELTA( * ), 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 I-th updated eigenvalue of a symmetric
</span><span class="comment">*</span><span class="comment">  rank-one modification to a diagonal matrix whose elements are
</span><span class="comment">*</span><span class="comment">  given in the array d, and that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">             D(i) &lt; D(j)  for  i &lt; j
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  and that RHO &gt; 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 )  +  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 &lt;= I &lt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D      (input) REAL 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, D(I) &lt; D(J)  for I &lt; J.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Z      (input) REAL 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) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">         If N .GT. 2, DELTA contains (D(j) - lambda_I) in its  j-th
</span><span class="comment">*</span><span class="comment">         component.  If N = 1, then DELTA(1) = 1. If N = 2, see <a name="SLAED5.52"></a><a href="slaed5.f.html#SLAED5.1">SLAED5</a>
</span><span class="comment">*</span><span class="comment">         for detail. The vector DELTA contains the information necessary
</span><span class="comment">*</span><span class="comment">         to construct the eigenvectors by <a name="SLAED3.54"></a><a href="slaed3.f.html#SLAED3.1">SLAED3</a> and <a name="SLAED9.54"></a><a href="slaed9.f.html#SLAED9.1">SLAED9</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RHO    (input) REAL
</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">  DLAM   (output) REAL
</span><span class="comment">*</span><span class="comment">         The computed lambda_I, the I-th updated eigenvalue.
</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">         &gt; 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 = 30 )
      REAL               ZERO, ONE, TWO, THREE, FOUR, EIGHT, TEN
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0,
     $                   THREE = 3.0E0, FOUR = 4.0E0, EIGHT = 8.0E0,
     $                   TEN = 10.0E0 )
<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
      REAL               A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
     $                   EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
     $                   RHOINV, TAU, TEMP, TEMP1, W
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Arrays ..
</span>      REAL               ZZ( 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      REAL               <a name="SLAMCH.109"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
      EXTERNAL           <a name="SLAMCH.110"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="SLAED5.113"></a><a href="slaed5.f.html#SLAED5.1">SLAED5</a>, <a name="SLAED6.113"></a><a href="slaed6.f.html#SLAED6.1">SLAED6</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>         DLAM = D( 1 ) + RHO*Z( 1 )*Z( 1 )
         DELTA( 1 ) = ONE
         RETURN
      END IF
      IF( N.EQ.2 ) THEN
         CALL <a name="SLAED5.135"></a><a href="slaed5.f.html#SLAED5.1">SLAED5</a>( I, D, Z, DELTA, RHO, DLAM )
         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="SLAMCH.141"></a><a href="slamch.f.html#SLAMCH.1">SLAMCH</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>         MIDPT = 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>         DO 10 J = 1, N
            DELTA( J ) = ( D( J )-D( I ) ) - MIDPT
   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 )
   20    CONTINUE
<span class="comment">*</span><span class="comment">
</span>         C = RHOINV + PSI
         W = C + Z( II )*Z( II ) / DELTA( II ) +
     $       Z( N )*Z( N ) / DELTA( N )
<span class="comment">*</span><span class="comment">
</span>         IF( W.LE.ZERO ) THEN
            TEMP = Z( N-1 )*Z( N-1 ) / ( D( N )-D( N-1 )+RHO ) +
     $             Z( N )*Z( N ) / RHO
            IF( C.LE.TEMP ) THEN
               TAU = RHO
            ELSE
               DEL = D( N ) - D( N-1 )
               A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
               B = Z( N )*Z( N )*DEL
               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)+RHO/2 &lt;= LAMBDA(N) &lt; D(N)+TAU &lt;= D(N)+RHO
</span><span class="comment">*</span><span class="comment">
</span>            DLTLB = MIDPT
            DLTUB = RHO
         ELSE
            DEL = D( N ) - D( N-1 )
            A = -C*DEL + Z( N-1 )*Z( N-1 ) + Z( N )*Z( N )
            B = Z( N )*Z( N )*DEL
            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
<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) &lt; D(N)+TAU &lt; LAMBDA(N) &lt; D(N)+RHO/2
</span><span class="comment">*</span><span class="comment">
</span>            DLTLB = ZERO
            DLTUB = MIDPT
         END IF
<span class="comment">*</span><span class="comment">
</span>         DO 30 J = 1, N
            DELTA( J ) = ( D( J )-D( I ) ) - TAU
   30    CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Evaluate PSI and the derivative DPSI
</span><span class="comment">*</span><span class="comment">
</span>         DPSI = ZERO
         PSI = ZERO
         ERRETM = ZERO
         DO 40 J = 1, II
            TEMP = Z( J ) / DELTA( J )
            PSI = PSI + Z( J )*TEMP
            DPSI = DPSI + TEMP*TEMP
            ERRETM = ERRETM + PSI
   40    CONTINUE
         ERRETM = ABS( ERRETM )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Evaluate PHI and the derivative DPHI
</span><span class="comment">*</span><span class="comment">
</span>         TEMP = Z( N ) / DELTA( N )
         PHI = Z( N )*TEMP
         DPHI = TEMP*TEMP
         ERRETM = EIGHT*( -PHI-PSI ) + ERRETM - PHI + RHOINV +
     $            ABS( TAU )*( DPSI+DPHI )
<span class="comment">*</span><span class="comment">
</span>         W = RHOINV + PHI + PSI
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Test for convergence
</span><span class="comment">*</span><span class="comment">
</span>         IF( ABS( W ).LE.EPS*ERRETM ) THEN
            DLAM = D( I ) + TAU
            GO TO 250
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( W.LE.ZERO ) THEN
            DLTLB = MAX( DLTLB, TAU )
         ELSE
            DLTUB = MIN( DLTUB, TAU )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Calculate the new step
</span><span class="comment">*</span><span class="comment">
</span>         NITER = NITER + 1
         C = W - DELTA( N-1 )*DPSI - DELTA( N )*DPHI
         A = ( DELTA( N-1 )+DELTA( N ) )*W -
     $       DELTA( N-1 )*DELTA( N )*( DPSI+DPHI )
         B = DELTA( N-1 )*DELTA( N )*W
         IF( C.LT.ZERO )
     $      C = ABS( C )
         IF( C.EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">          ETA = B/A
</span><span class="comment">*</span><span class="comment">           ETA = RHO - TAU
</span>            ETA = DLTUB - TAU
         ELSE IF( A.GE.ZERO ) THEN
            ETA = ( A+SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
         ELSE
            ETA = TWO*B / ( A-SQRT( ABS( A*A-FOUR*B*C ) ) )
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Note, eta should be positive if w is negative, and
</span><span class="comment">*</span><span class="comment">        eta should be negative otherwise. However,
123

slaed4.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 slaed4.f.html 源码文件，采用 HTML 编程语言编写，共 869 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?