dlaed6.f.html

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

HTML
352 行
字号
12
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dlaed6.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="DLAED6.1"></a><a href="dlaed6.f.html#DLAED6.1">DLAED6</a>( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK routine (version 3.1.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">     February 2007
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      LOGICAL            ORGATI
      INTEGER            INFO, KNITER
      DOUBLE PRECISION   FINIT, RHO, TAU
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      DOUBLE PRECISION   D( 3 ), Z( 3 )
<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="DLAED6.19"></a><a href="dlaed6.f.html#DLAED6.1">DLAED6</a> computes the positive or negative root (closest to the origin)
</span><span class="comment">*</span><span class="comment">  of
</span><span class="comment">*</span><span class="comment">                   z(1)        z(2)        z(3)
</span><span class="comment">*</span><span class="comment">  f(x) =   rho + --------- + ---------- + ---------
</span><span class="comment">*</span><span class="comment">                  d(1)-x      d(2)-x      d(3)-x
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  It is assumed that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        if ORGATI = .true. the root is between d(2) and d(3);
</span><span class="comment">*</span><span class="comment">        otherwise it is between d(1) and d(2)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  This routine will be called by <a name="DLAED4.30"></a><a href="dlaed4.f.html#DLAED4.1">DLAED4</a> when necessary. In most cases,
</span><span class="comment">*</span><span class="comment">  the root sought is the smallest in magnitude, though it might not be
</span><span class="comment">*</span><span class="comment">  in some extremely rare situations.
</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">  KNITER       (input) INTEGER
</span><span class="comment">*</span><span class="comment">               Refer to <a name="DLAED4.38"></a><a href="dlaed4.f.html#DLAED4.1">DLAED4</a> for its significance.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ORGATI       (input) LOGICAL
</span><span class="comment">*</span><span class="comment">               If ORGATI is true, the needed root is between d(2) and
</span><span class="comment">*</span><span class="comment">               d(3); otherwise it is between d(1) and d(2).  See
</span><span class="comment">*</span><span class="comment">               <a name="DLAED4.43"></a><a href="dlaed4.f.html#DLAED4.1">DLAED4</a> for further details.
</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">               Refer to the equation f(x) above.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D            (input) DOUBLE PRECISION array, dimension (3)
</span><span class="comment">*</span><span class="comment">               D satisfies d(1) &lt; d(2) &lt; d(3).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Z            (input) DOUBLE PRECISION array, dimension (3)
</span><span class="comment">*</span><span class="comment">               Each of the elements in z must be positive.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  FINIT        (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">               The value of f at 0. It is more accurate than the one
</span><span class="comment">*</span><span class="comment">               evaluated inside this routine (if someone wants to do
</span><span class="comment">*</span><span class="comment">               so).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  TAU          (output) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">               The root of the equation f(x).
</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, failure to converge
</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">  30/06/99: 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">  10/02/03: This version has a few statements commented out for thread
</span><span class="comment">*</span><span class="comment">  safety (machine parameters are computed on each entry). SJH.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  05/10/06: Modified from a new version of Ren-Cang Li, use
</span><span class="comment">*</span><span class="comment">     Gragg-Thornton-Warner cubic convergent scheme for better stability.
</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 = 40 )
      DOUBLE PRECISION   ZERO, ONE, TWO, THREE, FOUR, EIGHT
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
     $                   THREE = 3.0D0, FOUR = 4.0D0, EIGHT = 8.0D0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMCH.89"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DLAMCH.90"></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">     .. Local Arrays ..
</span>      DOUBLE PRECISION   DSCALE( 3 ), ZSCALE( 3 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            SCALE
      INTEGER            I, ITER, NITER
      DOUBLE PRECISION   A, B, BASE, C, DDF, DF, EPS, ERRETM, ETA, F,
     $                   FC, SCLFAC, SCLINV, SMALL1, SMALL2, SMINV1,
     $                   SMINV2, TEMP, TEMP1, TEMP2, TEMP3, TEMP4, 
     $                   LBD, UBD
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, INT, LOG, 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>      INFO = 0
<span class="comment">*</span><span class="comment">
</span>      IF( ORGATI ) THEN
         LBD = D(2)
         UBD = D(3)
      ELSE
         LBD = D(1)
         UBD = D(2)
      END IF
      IF( FINIT .LT. ZERO )THEN
         LBD = ZERO
      ELSE
         UBD = ZERO 
      END IF
<span class="comment">*</span><span class="comment">
</span>      NITER = 1
      TAU = ZERO
      IF( KNITER.EQ.2 ) THEN
         IF( ORGATI ) THEN
            TEMP = ( D( 3 )-D( 2 ) ) / TWO
            C = RHO + Z( 1 ) / ( ( D( 1 )-D( 2 ) )-TEMP )
            A = C*( D( 2 )+D( 3 ) ) + Z( 2 ) + Z( 3 )
            B = C*D( 2 )*D( 3 ) + Z( 2 )*D( 3 ) + Z( 3 )*D( 2 )
         ELSE
            TEMP = ( D( 1 )-D( 2 ) ) / TWO
            C = RHO + Z( 3 ) / ( ( D( 3 )-D( 2 ) )-TEMP )
            A = C*( D( 1 )+D( 2 ) ) + Z( 1 ) + Z( 2 )
            B = C*D( 1 )*D( 2 ) + Z( 1 )*D( 2 ) + Z( 2 )*D( 1 )
         END IF
         TEMP = MAX( ABS( A ), ABS( B ), ABS( C ) )
         A = A / TEMP
         B = B / TEMP
         C = C / TEMP
         IF( C.EQ.ZERO ) THEN
            TAU = B / A
         ELSE IF( A.LE.ZERO ) THEN
            TAU = ( A-SQRT( ABS( A*A-FOUR*B*C ) ) ) / ( TWO*C )
         ELSE
            TAU = TWO*B / ( A+SQRT( ABS( A*A-FOUR*B*C ) ) )
         END IF
         IF( TAU .LT. LBD .OR. TAU .GT. UBD )
     $      TAU = ( LBD+UBD )/TWO
         IF( D(1).EQ.TAU .OR. D(2).EQ.TAU .OR. D(3).EQ.TAU ) THEN
            TAU = ZERO
         ELSE
            TEMP = FINIT + TAU*Z(1)/( D(1)*( D( 1 )-TAU ) ) +
     $                     TAU*Z(2)/( D(2)*( D( 2 )-TAU ) ) +
     $                     TAU*Z(3)/( D(3)*( D( 3 )-TAU ) )
            IF( TEMP .LE. ZERO )THEN
               LBD = TAU
12

dlaed6.f.html - 源码说明

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

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

⌨️ 快捷键说明

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