zlahef.f.html

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

HTML
672 行
字号
123
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>zlahef.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="ZLAHEF.1"></a><a href="zlahef.f.html#ZLAHEF.1">ZLAHEF</a>( UPLO, N, NB, KB, A, LDA, IPIV, W, LDW, 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, KB, LDA, LDW, N, NB
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      INTEGER            IPIV( * )
      COMPLEX*16         A( LDA, * ), W( LDW, * )
<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="ZLAHEF.19"></a><a href="zlahef.f.html#ZLAHEF.1">ZLAHEF</a> computes a partial factorization of a complex Hermitian
</span><span class="comment">*</span><span class="comment">  matrix A using the Bunch-Kaufman diagonal pivoting method. The
</span><span class="comment">*</span><span class="comment">  partial factorization has the form:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A  =  ( I  U12 ) ( A11  0  ) (  I    0   )  if UPLO = 'U', or:
</span><span class="comment">*</span><span class="comment">        ( 0  U22 ) (  0   D  ) ( U12' U22' )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A  =  ( L11  0 ) (  D   0  ) ( L11' L21' )  if UPLO = 'L'
</span><span class="comment">*</span><span class="comment">        ( L21  I ) (  0  A22 ) (  0    I   )
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where the order of D is at most NB. The actual order is returned in
</span><span class="comment">*</span><span class="comment">  the argument KB, and is either NB or NB-1, or N if N &lt;= NB.
</span><span class="comment">*</span><span class="comment">  Note that U' denotes the conjugate transpose of U.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="ZLAHEF.33"></a><a href="zlahef.f.html#ZLAHEF.1">ZLAHEF</a> is an auxiliary routine called by <a name="ZHETRF.33"></a><a href="zhetrf.f.html#ZHETRF.1">ZHETRF</a>. It uses blocked code
</span><span class="comment">*</span><span class="comment">  (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
</span><span class="comment">*</span><span class="comment">  A22 (if UPLO = 'L').
</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">          Specifies whether the upper or lower triangular part of the
</span><span class="comment">*</span><span class="comment">          Hermitian matrix A is stored:
</span><span class="comment">*</span><span class="comment">          = 'U':  Upper triangular
</span><span class="comment">*</span><span class="comment">          = 'L':  Lower triangular
</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">  NB      (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The maximum number of columns of the matrix A that should be
</span><span class="comment">*</span><span class="comment">          factored.  NB should be at least 2 to allow for 2-by-2 pivot
</span><span class="comment">*</span><span class="comment">          blocks.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  KB      (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns of A that were actually factored.
</span><span class="comment">*</span><span class="comment">          KB is either NB-1 or NB, or N if N &lt;= NB.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
</span><span class="comment">*</span><span class="comment">          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
</span><span class="comment">*</span><span class="comment">          n-by-n upper triangular part of A contains the upper
</span><span class="comment">*</span><span class="comment">          triangular part of the matrix A, and the strictly lower
</span><span class="comment">*</span><span class="comment">          triangular part of A is not referenced.  If UPLO = 'L', the
</span><span class="comment">*</span><span class="comment">          leading n-by-n lower triangular part of A contains the lower
</span><span class="comment">*</span><span class="comment">          triangular part of the matrix A, and the strictly upper
</span><span class="comment">*</span><span class="comment">          triangular part of A is not referenced.
</span><span class="comment">*</span><span class="comment">          On exit, A contains details of the partial factorization.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A.  LDA &gt;= max(1,N).
</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 UPLO = 'U', only the last KB elements of IPIV are set;
</span><span class="comment">*</span><span class="comment">          if UPLO = 'L', only the first KB elements are set.
</span><span class="comment">*</span><span class="comment">
</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">  W       (workspace) COMPLEX*16 array, dimension (LDW,NB)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDW     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array W.  LDW &gt;= max(1,N).
</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 = k, D(k,k) 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.
</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 )
      COMPLEX*16         CONE
      PARAMETER          ( CONE = ( 1.0D+0, 0.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>      INTEGER            IMAX, J, JB, JJ, JMAX, JP, K, KK, KKW, KP,
     $                   KSTEP, KW
      DOUBLE PRECISION   ABSAKK, ALPHA, COLMAX, R1, ROWMAX, T
      COMPLEX*16         D11, D21, D22, Z
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      LOGICAL            <a name="LSAME.112"></a><a href="lsame.f.html#LSAME.1">LSAME</a>
      INTEGER            IZAMAX
      EXTERNAL           <a name="LSAME.114"></a><a href="lsame.f.html#LSAME.1">LSAME</a>, IZAMAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           ZCOPY, ZDSCAL, ZGEMM, ZGEMV, <a name="ZLACGV.117"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>, ZSWAP
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, DBLE, DCONJG, DIMAG, MAX, MIN, SQRT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Functions ..
</span>      DOUBLE PRECISION   CABS1
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Statement Function definitions ..
</span>      CABS1( Z ) = ABS( DBLE( Z ) ) + ABS( DIMAG( Z ) )
<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><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( <a name="LSAME.136"></a><a href="lsame.f.html#LSAME.1">LSAME</a>( UPLO, <span class="string">'U'</span> ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Factorize the trailing columns of A using the upper triangle
</span><span class="comment">*</span><span class="comment">        of A and working backwards, and compute the matrix W = U12*D
</span><span class="comment">*</span><span class="comment">        for use in updating A11 (note that conjg(W) is actually stored)
</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 in steps of 1 or 2
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        KW is the column of W which corresponds to column K of A
</span><span class="comment">*</span><span class="comment">
</span>         K = N
   10    CONTINUE
         KW = NB + K - N
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Exit from loop
</span><span class="comment">*</span><span class="comment">
</span>         IF( ( K.LE.N-NB+1 .AND. NB.LT.N ) .OR. K.LT.1 )
     $      GO TO 30
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy column K of A to column KW of W and update it
</span><span class="comment">*</span><span class="comment">
</span>         CALL ZCOPY( K-1, A( 1, K ), 1, W( 1, KW ), 1 )
         W( K, KW ) = DBLE( A( K, K ) )
         IF( K.LT.N ) THEN
            CALL ZGEMV( <span class="string">'No transpose'</span>, K, N-K, -CONE, A( 1, K+1 ), LDA,
     $                  W( K, KW+1 ), LDW, CONE, W( 1, KW ), 1 )
            W( K, KW ) = DBLE( W( K, KW ) )
         END IF
<span class="comment">*</span><span class="comment">
</span>         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
</span><span class="comment">*</span><span class="comment">
</span>         ABSAKK = ABS( DBLE( W( K, KW ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        IMAX is the row-index of the largest off-diagonal element in
</span><span class="comment">*</span><span class="comment">        column K, and COLMAX is its absolute value
</span><span class="comment">*</span><span class="comment">
</span>         IF( K.GT.1 ) THEN
            IMAX = IZAMAX( K-1, W( 1, KW ), 1 )
            COLMAX = CABS1( W( IMAX, KW ) )
         ELSE
            COLMAX = ZERO
         END IF
<span class="comment">*</span><span class="comment">
</span>         IF( MAX( ABSAKK, COLMAX ).EQ.ZERO ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           Column K is zero: set INFO and continue
</span><span class="comment">*</span><span class="comment">
</span>            IF( INFO.EQ.0 )
     $         INFO = K
            KP = K
            A( K, K ) = DBLE( A( K, K ) )
         ELSE
            IF( ABSAKK.GE.ALPHA*COLMAX ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              no interchange, use 1-by-1 pivot block
</span><span class="comment">*</span><span class="comment">
</span>               KP = K
            ELSE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">              Copy column IMAX to column KW-1 of W and update it
</span><span class="comment">*</span><span class="comment">
</span>               CALL ZCOPY( IMAX-1, A( 1, IMAX ), 1, W( 1, KW-1 ), 1 )
               W( IMAX, KW-1 ) = DBLE( A( IMAX, IMAX ) )
               CALL ZCOPY( K-IMAX, A( IMAX, IMAX+1 ), LDA,
     $                     W( IMAX+1, KW-1 ), 1 )
               CALL <a name="ZLACGV.204"></a><a href="zlacgv.f.html#ZLACGV.1">ZLACGV</a>( K-IMAX, W( IMAX+1, KW-1 ), 1 )
               IF( K.LT.N ) THEN
123

zlahef.f.html - 源码说明

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

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

⌨️ 快捷键说明

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