dsptrf.f.html

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

HTML
572 行
字号
123
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>dsptrf.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="DSPTRF.1"></a><a href="dsptrf.f.html#DSPTRF.1">DSPTRF</a>( UPLO, N, AP, IPIV, 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>      CHARACTER          UPLO
      INTEGER            INFO, N
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IPIV( * )
      DOUBLE PRECISION   AP( * )
<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="DSPTRF.19"></a><a href="dsptrf.f.html#DSPTRF.1">DSPTRF</a> computes the factorization of a real symmetric matrix A stored
</span><span class="comment">*</span><span class="comment">  in packed format using the Bunch-Kaufman diagonal pivoting method:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     A = U*D*U**T  or  A = L*D*L**T
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where U (or L) is a product of permutation and unit upper (lower)
</span><span class="comment">*</span><span class="comment">  triangular matrices, and D is symmetric and block diagonal with
</span><span class="comment">*</span><span class="comment">  1-by-1 and 2-by-2 diagonal blocks.
</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">  UPLO    (input) CHARACTER*1
</span><span class="comment">*</span><span class="comment">          = 'U':  Upper triangle of A is stored;
</span><span class="comment">*</span><span class="comment">          = 'L':  Lower triangle of A is stored.
</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 A.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  AP      (input/output) DOUBLE PRECISION array, dimension (N*(N+1)/2)
</span><span class="comment">*</span><span class="comment">          On entry, the upper or lower triangle of the symmetric matrix
</span><span class="comment">*</span><span class="comment">          A, packed columnwise in a linear array.  The j-th column of A
</span><span class="comment">*</span><span class="comment">          is stored in the array AP as follows:
</span><span class="comment">*</span><span class="comment">          if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1&lt;=i&lt;=j;
</span><span class="comment">*</span><span class="comment">          if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j&lt;=i&lt;=n.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          On exit, the block diagonal matrix D and the multipliers used
</span><span class="comment">*</span><span class="comment">          to obtain the factor U or L, stored as a packed triangular
</span><span class="comment">*</span><span class="comment">          matrix overwriting A (see below for further details).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IPIV    (output) INTEGER array, dimension (N)
</span><span class="comment">*</span><span class="comment">          Details of the interchanges and the block structure of D.
</span><span class="comment">*</span><span class="comment">          If IPIV(k) &gt; 0, then rows and columns k and IPIV(k) were
</span><span class="comment">*</span><span class="comment">          interchanged and D(k,k) is a 1-by-1 diagonal block.
</span><span class="comment">*</span><span class="comment">          If UPLO = 'U' and IPIV(k) = IPIV(k-1) &lt; 0, then rows and
</span><span class="comment">*</span><span class="comment">          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
</span><span class="comment">*</span><span class="comment">          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
</span><span class="comment">*</span><span class="comment">          IPIV(k+1) &lt; 0, then rows and columns k+1 and -IPIV(k) were
</span><span class="comment">*</span><span class="comment">          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
</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">          &lt; 0: if INFO = -i, the i-th argument had an illegal value
</span><span class="comment">*</span><span class="comment">          &gt; 0: if INFO = i, D(i,i) is exactly zero.  The factorization
</span><span class="comment">*</span><span class="comment">               has been completed, but the block diagonal matrix D is
</span><span class="comment">*</span><span class="comment">               exactly singular, and division by zero will occur if it
</span><span class="comment">*</span><span class="comment">               is used to solve a system of equations.
</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">  5-96 - Based on modifications by J. Lewis, Boeing Computer Services
</span><span class="comment">*</span><span class="comment">         Company
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If UPLO = 'U', then A = U*D*U', where
</span><span class="comment">*</span><span class="comment">     U = P(n)*U(n)* ... *P(k)U(k)* ...,
</span><span class="comment">*</span><span class="comment">  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
</span><span class="comment">*</span><span class="comment">  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
</span><span class="comment">*</span><span class="comment">  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
</span><span class="comment">*</span><span class="comment">  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
</span><span class="comment">*</span><span class="comment">  that if the diagonal block D(k) is of order s (s = 1 or 2), then
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">             (   I    v    0   )   k-s
</span><span class="comment">*</span><span class="comment">     U(k) =  (   0    I    0   )   s
</span><span class="comment">*</span><span class="comment">             (   0    0    I   )   n-k
</span><span class="comment">*</span><span class="comment">                k-s   s   n-k
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
</span><span class="comment">*</span><span class="comment">  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
</span><span class="comment">*</span><span class="comment">  and A(k,k), and v overwrites A(1:k-2,k-1:k).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If UPLO = 'L', then A = L*D*L', where
</span><span class="comment">*</span><span class="comment">     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
</span><span class="comment">*</span><span class="comment">  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
</span><span class="comment">*</span><span class="comment">  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
</span><span class="comment">*</span><span class="comment">  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
</span><span class="comment">*</span><span class="comment">  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
</span><span class="comment">*</span><span class="comment">  that if the diagonal block D(k) is of order s (s = 1 or 2), then
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">             (   I    0     0   )  k-1
</span><span class="comment">*</span><span class="comment">     L(k) =  (   0    I     0   )  s
</span><span class="comment">*</span><span class="comment">             (   0    v     I   )  n-k-s+1
</span><span class="comment">*</span><span class="comment">                k-1   s  n-k-s+1
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
</span><span class="comment">*</span><span class="comment">  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
</span><span class="comment">*</span><span class="comment">  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+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   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
      DOUBLE PRECISION   EIGHT, SEVTEN
      PARAMETER          ( EIGHT = 8.0D+0, SEVTEN = 17.0D+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            UPPER
      INTEGER            I, IMAX, J, JMAX, K, KC, KK, KNC, KP, KPC,
     $                   KSTEP, KX, NPP
      DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, D11, D12, D21, D22, R1,
     $                   ROWMAX, T, WK, WKM1, WKP1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.123"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            IDAMAX
      EXTERNAL           <a name="LSAME.125"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, IDAMAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           DSCAL, DSPR, DSWAP, <a name="XERBLA.128"></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">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, 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 parameters.
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      UPPER = <a name="LSAME.138"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> )
      IF( .NOT.UPPER .AND. .NOT.<a name="LSAME.139"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'L'</span> ) ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      END IF
      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.145"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="DSPTRF.145"></a><a href="dsptrf.f.html#DSPTRF.1">DSPTRF</a>'</span>, -INFO )
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Initialize ALPHA for use in choosing pivot block size.
</span><span class="comment">*</span><span class="comment">
</span>      ALPHA = ( ONE+SQRT( SEVTEN ) ) / EIGHT
<span class="comment">*</span><span class="comment">
</span>      IF( UPPER ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Factorize A as U*D*U' using the upper triangle of A
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        K is the main loop index, decreasing from N to 1 in steps of
</span><span class="comment">*</span><span class="comment">        1 or 2
</span><span class="comment">*</span><span class="comment">
</span>         K = N
         KC = ( N-1 )*N / 2 + 1
   10    CONTINUE
         KNC = KC
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        If K &lt; 1, exit from loop
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.LT.1 )
     $      GO TO 110
         KSTEP = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Determine rows and columns to be interchanged and whether
</span><span class="comment">*</span><span class="comment">        a 1-by-1 or 2-by-2 pivot block will be used
123

dsptrf.f.html - 源码说明

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

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

⌨️ 快捷键说明

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