dlarrf.f.html

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

HTML
398 行
字号
12
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dlarrf.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="DLARRF.1"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a>( N, D, L, LD, CLSTRT, CLEND,
     $                   W, WGAP, WERR,
     $                   SPDIAM, CLGAPL, CLGAPR, PIVMIN, SIGMA,
     $                   DPLUS, LPLUS, 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            CLSTRT, CLEND, INFO, N
      DOUBLE PRECISION   CLGAPL, CLGAPR, PIVMIN, SIGMA, SPDIAM
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      DOUBLE PRECISION   D( * ), DPLUS( * ), L( * ), LD( * ),
     $          LPLUS( * ), W( * ), WGAP( * ), WERR( * ), WORK( * )
<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">  Given the initial representation L D L^T and its cluster of close
</span><span class="comment">*</span><span class="comment">  eigenvalues (in a relative measure), W( CLSTRT ), W( CLSTRT+1 ), ...
</span><span class="comment">*</span><span class="comment">  W( CLEND ), <a name="DLARRF.24"></a><a href="dlarrf.f.html#DLARRF.1">DLARRF</a> finds a new relatively robust representation
</span><span class="comment">*</span><span class="comment">  L D L^T - SIGMA I = L(+) D(+) L(+)^T such that at least one of the
</span><span class="comment">*</span><span class="comment">  eigenvalues of L(+) D(+) L(+)^T is relatively isolated.
</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 order of the matrix (subblock, if the matrix splitted).
</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 N diagonal elements of the diagonal matrix D.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  L       (input) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          The (N-1) subdiagonal elements of the unit bidiagonal
</span><span class="comment">*</span><span class="comment">          matrix L.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LD      (input) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          The (N-1) elements L(i)*D(i).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  CLSTRT  (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The index of the first eigenvalue in the cluster.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  CLEND   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The index of the last eigenvalue in the cluster.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  W       (input) DOUBLE PRECISION array, dimension &gt;=  (CLEND-CLSTRT+1)
</span><span class="comment">*</span><span class="comment">          The eigenvalue APPROXIMATIONS of L D L^T in ascending order.
</span><span class="comment">*</span><span class="comment">          W( CLSTRT ) through W( CLEND ) form the cluster of relatively
</span><span class="comment">*</span><span class="comment">          close eigenalues.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WGAP    (input/output) DOUBLE PRECISION array, dimension &gt;=  (CLEND-CLSTRT+1)
</span><span class="comment">*</span><span class="comment">          The separation from the right neighbor eigenvalue in W.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WERR    (input) DOUBLE PRECISION array, dimension &gt;=  (CLEND-CLSTRT+1)
</span><span class="comment">*</span><span class="comment">          WERR contain the semiwidth of the uncertainty
</span><span class="comment">*</span><span class="comment">          interval of the corresponding eigenvalue APPROXIMATION in W
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  SPDIAM (input) estimate of the spectral diameter obtained from the
</span><span class="comment">*</span><span class="comment">          Gerschgorin intervals
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  CLGAPL, CLGAPR (input) absolute gap on each end of the cluster.
</span><span class="comment">*</span><span class="comment">          Set by the calling routine to protect against shifts too close
</span><span class="comment">*</span><span class="comment">          to eigenvalues outside the cluster.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  PIVMIN  (input) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          The minimum pivot allowed in the Sturm sequence.
</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 shift used to form L(+) D(+) L(+)^T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  DPLUS   (output) DOUBLE PRECISION array, dimension (N)
</span><span class="comment">*</span><span class="comment">          The N diagonal elements of the diagonal matrix D(+).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LPLUS   (output) DOUBLE PRECISION array, dimension (N-1)
</span><span class="comment">*</span><span class="comment">          The first (N-1) elements of LPLUS contain the subdiagonal
</span><span class="comment">*</span><span class="comment">          elements of the unit bidiagonal matrix L(+).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)
</span><span class="comment">*</span><span class="comment">          Workspace.
</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">     Beresford Parlett, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Jim Demmel, University of California, Berkeley, USA
</span><span class="comment">*</span><span class="comment">     Inderjit Dhillon, University of Texas, Austin, USA
</span><span class="comment">*</span><span class="comment">     Osni Marques, LBNL/NERSC, USA
</span><span class="comment">*</span><span class="comment">     Christof Voemel, University of California, 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>      DOUBLE PRECISION   FOUR, MAXGROWTH1, MAXGROWTH2, ONE, QUART, TWO,
     $                   ZERO
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0,
     $                     FOUR = 4.0D0, QUART = 0.25D0,
     $                     MAXGROWTH1 = 8.D0,
     $                     MAXGROWTH2 = 8.D0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL   DORRR1, FORCER, NOFAIL, SAWNAN1, SAWNAN2, TRYRRR1
      INTEGER            I, INDX, KTRY, KTRYMAX, SLEFT, SRIGHT, SHIFT
      PARAMETER          ( KTRYMAX = 1, SLEFT = 1, SRIGHT = 2 )
      DOUBLE PRECISION   AVGAP, BESTSHIFT, CLWDTH, EPS, FACT, FAIL,
     $                   FAIL2, GROWTHBOUND, LDELTA, LDMAX, LSIGMA,
     $                   MAX1, MAX2, MINGAP, OLDP, PROD, RDELTA, RDMAX,
     $                   RRR1, RRR2, RSIGMA, S, SMLGROWTH, TMP, ZNM2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL <a name="DISNAN.115"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>
      DOUBLE PRECISION   <a name="DLAMCH.116"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>
      EXTERNAL           <a name="DISNAN.117"></a><a href="disnan.f.html#DISNAN.1">DISNAN</a>, <a name="DLAMCH.117"></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">     .. External Subroutines ..
</span>      EXTERNAL           DCOPY
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS
<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
      FACT = DBLE(2**KTRYMAX)
      EPS = <a name="DLAMCH.129"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'Precision'</span> )
      SHIFT = 0
      FORCER = .FALSE.


<span class="comment">*</span><span class="comment">     Note that we cannot guarantee that for any of the shifts tried,
</span><span class="comment">*</span><span class="comment">     the factorization has a small or even moderate element growth.
</span><span class="comment">*</span><span class="comment">     There could be Ritz values at both ends of the cluster and despite
</span><span class="comment">*</span><span class="comment">     backing off, there are examples where all factorizations tried
</span><span class="comment">*</span><span class="comment">     (in IEEE mode, allowing zero pivots &amp; infinities) have INFINITE
</span><span class="comment">*</span><span class="comment">     element growth.
</span><span class="comment">*</span><span class="comment">     For this reason, we should use PIVMIN in this subroutine so that at
</span><span class="comment">*</span><span class="comment">     least the L D L^T factorization exists. It can be checked afterwards
</span><span class="comment">*</span><span class="comment">     whether the element growth caused bad residuals/orthogonality.
</span>
<span class="comment">*</span><span class="comment">     Decide whether the code should accept the best among all
</span><span class="comment">*</span><span class="comment">     representations despite large element growth or signal INFO=1
</span>      NOFAIL = .TRUE.
<span class="comment">*</span><span class="comment">
</span>
<span class="comment">*</span><span class="comment">     Compute the average gap length of the cluster
</span>      CLWDTH = ABS(W(CLEND)-W(CLSTRT)) + WERR(CLEND) + WERR(CLSTRT)
      AVGAP = CLWDTH / DBLE(CLEND-CLSTRT)
      MINGAP = MIN(CLGAPL, CLGAPR)
<span class="comment">*</span><span class="comment">     Initial values for shifts to both ends of cluster
</span>      LSIGMA = MIN(W( CLSTRT ),W( CLEND )) - WERR( CLSTRT )
      RSIGMA = MAX(W( CLSTRT ),W( CLEND )) + WERR( CLEND )

<span class="comment">*</span><span class="comment">     Use a small fudge to make sure that we really shift to the outside
</span>      LSIGMA = LSIGMA - ABS(LSIGMA)* FOUR * EPS
      RSIGMA = RSIGMA + ABS(RSIGMA)* FOUR * EPS

<span class="comment">*</span><span class="comment">     Compute upper bounds for how much to back off the initial shifts
</span>      LDMAX = QUART * MINGAP + TWO * PIVMIN
      RDMAX = QUART * MINGAP + TWO * PIVMIN

      LDELTA = MAX(AVGAP,WGAP( CLSTRT ))/FACT
      RDELTA = MAX(AVGAP,WGAP( CLEND-1 ))/FACT
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize the record of the best representation found
</span><span class="comment">*</span><span class="comment">
</span>      S = <a name="DLAMCH.170"></a><a href="dlamch.f.html#DLAMCH.1">DLAMCH</a>( <span class="string">'S'</span> )
      SMLGROWTH = ONE / S
      FAIL = DBLE(N-1)*MINGAP/(SPDIAM*EPS)
      FAIL2 = DBLE(N-1)*MINGAP/(SPDIAM*SQRT(EPS))
      BESTSHIFT = LSIGMA
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     while (KTRY &lt;= KTRYMAX)
</span>      KTRY = 0
      GROWTHBOUND = MAXGROWTH1*SPDIAM

 5    CONTINUE
12

dlarrf.f.html - 源码说明

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

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

⌨️ 快捷键说明

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