LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zhetrf.f
Go to the documentation of this file.
00001 *> \brief \b ZHETRF
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZHETRF + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZHETRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
00022 * 
00023 *       .. Scalar Arguments ..
00024 *       CHARACTER          UPLO
00025 *       INTEGER            INFO, LDA, LWORK, N
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            IPIV( * )
00029 *       COMPLEX*16         A( LDA, * ), WORK( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> ZHETRF computes the factorization of a complex Hermitian matrix A
00039 *> using the Bunch-Kaufman diagonal pivoting method.  The form of the
00040 *> factorization is
00041 *>
00042 *>    A = U*D*U**H  or  A = L*D*L**H
00043 *>
00044 *> where U (or L) is a product of permutation and unit upper (lower)
00045 *> triangular matrices, and D is Hermitian and block diagonal with
00046 *> 1-by-1 and 2-by-2 diagonal blocks.
00047 *>
00048 *> This is the blocked version of the algorithm, calling Level 3 BLAS.
00049 *> \endverbatim
00050 *
00051 *  Arguments:
00052 *  ==========
00053 *
00054 *> \param[in] UPLO
00055 *> \verbatim
00056 *>          UPLO is CHARACTER*1
00057 *>          = 'U':  Upper triangle of A is stored;
00058 *>          = 'L':  Lower triangle of A is stored.
00059 *> \endverbatim
00060 *>
00061 *> \param[in] N
00062 *> \verbatim
00063 *>          N is INTEGER
00064 *>          The order of the matrix A.  N >= 0.
00065 *> \endverbatim
00066 *>
00067 *> \param[in,out] A
00068 *> \verbatim
00069 *>          A is COMPLEX*16 array, dimension (LDA,N)
00070 *>          On entry, the Hermitian matrix A.  If UPLO = 'U', the leading
00071 *>          N-by-N upper triangular part of A contains the upper
00072 *>          triangular part of the matrix A, and the strictly lower
00073 *>          triangular part of A is not referenced.  If UPLO = 'L', the
00074 *>          leading N-by-N lower triangular part of A contains the lower
00075 *>          triangular part of the matrix A, and the strictly upper
00076 *>          triangular part of A is not referenced.
00077 *>
00078 *>          On exit, the block diagonal matrix D and the multipliers used
00079 *>          to obtain the factor U or L (see below for further details).
00080 *> \endverbatim
00081 *>
00082 *> \param[in] LDA
00083 *> \verbatim
00084 *>          LDA is INTEGER
00085 *>          The leading dimension of the array A.  LDA >= max(1,N).
00086 *> \endverbatim
00087 *>
00088 *> \param[out] IPIV
00089 *> \verbatim
00090 *>          IPIV is INTEGER array, dimension (N)
00091 *>          Details of the interchanges and the block structure of D.
00092 *>          If IPIV(k) > 0, then rows and columns k and IPIV(k) were
00093 *>          interchanged and D(k,k) is a 1-by-1 diagonal block.
00094 *>          If UPLO = 'U' and IPIV(k) = IPIV(k-1) < 0, then rows and
00095 *>          columns k-1 and -IPIV(k) were interchanged and D(k-1:k,k-1:k)
00096 *>          is a 2-by-2 diagonal block.  If UPLO = 'L' and IPIV(k) =
00097 *>          IPIV(k+1) < 0, then rows and columns k+1 and -IPIV(k) were
00098 *>          interchanged and D(k:k+1,k:k+1) is a 2-by-2 diagonal block.
00099 *> \endverbatim
00100 *>
00101 *> \param[out] WORK
00102 *> \verbatim
00103 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
00104 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00105 *> \endverbatim
00106 *>
00107 *> \param[in] LWORK
00108 *> \verbatim
00109 *>          LWORK is INTEGER
00110 *>          The length of WORK.  LWORK >=1.  For best performance
00111 *>          LWORK >= N*NB, where NB is the block size returned by ILAENV.
00112 *> \endverbatim
00113 *>
00114 *> \param[out] INFO
00115 *> \verbatim
00116 *>          INFO is INTEGER
00117 *>          = 0:  successful exit
00118 *>          < 0:  if INFO = -i, the i-th argument had an illegal value
00119 *>          > 0:  if INFO = i, D(i,i) is exactly zero.  The factorization
00120 *>                has been completed, but the block diagonal matrix D is
00121 *>                exactly singular, and division by zero will occur if it
00122 *>                is used to solve a system of equations.
00123 *> \endverbatim
00124 *
00125 *  Authors:
00126 *  ========
00127 *
00128 *> \author Univ. of Tennessee 
00129 *> \author Univ. of California Berkeley 
00130 *> \author Univ. of Colorado Denver 
00131 *> \author NAG Ltd. 
00132 *
00133 *> \date November 2011
00134 *
00135 *> \ingroup complex16HEcomputational
00136 *
00137 *> \par Further Details:
00138 *  =====================
00139 *>
00140 *> \verbatim
00141 *>
00142 *>  If UPLO = 'U', then A = U*D*U**H, where
00143 *>     U = P(n)*U(n)* ... *P(k)U(k)* ...,
00144 *>  i.e., U is a product of terms P(k)*U(k), where k decreases from n to
00145 *>  1 in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00146 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00147 *>  defined by IPIV(k), and U(k) is a unit upper triangular matrix, such
00148 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00149 *>
00150 *>             (   I    v    0   )   k-s
00151 *>     U(k) =  (   0    I    0   )   s
00152 *>             (   0    0    I   )   n-k
00153 *>                k-s   s   n-k
00154 *>
00155 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(1:k-1,k).
00156 *>  If s = 2, the upper triangle of D(k) overwrites A(k-1,k-1), A(k-1,k),
00157 *>  and A(k,k), and v overwrites A(1:k-2,k-1:k).
00158 *>
00159 *>  If UPLO = 'L', then A = L*D*L**H, where
00160 *>     L = P(1)*L(1)* ... *P(k)*L(k)* ...,
00161 *>  i.e., L is a product of terms P(k)*L(k), where k increases from 1 to
00162 *>  n in steps of 1 or 2, and D is a block diagonal matrix with 1-by-1
00163 *>  and 2-by-2 diagonal blocks D(k).  P(k) is a permutation matrix as
00164 *>  defined by IPIV(k), and L(k) is a unit lower triangular matrix, such
00165 *>  that if the diagonal block D(k) is of order s (s = 1 or 2), then
00166 *>
00167 *>             (   I    0     0   )  k-1
00168 *>     L(k) =  (   0    I     0   )  s
00169 *>             (   0    v     I   )  n-k-s+1
00170 *>                k-1   s  n-k-s+1
00171 *>
00172 *>  If s = 1, D(k) overwrites A(k,k), and v overwrites A(k+1:n,k).
00173 *>  If s = 2, the lower triangle of D(k) overwrites A(k,k), A(k+1,k),
00174 *>  and A(k+1,k+1), and v overwrites A(k+2:n,k:k+1).
00175 *> \endverbatim
00176 *>
00177 *  =====================================================================
00178       SUBROUTINE ZHETRF( UPLO, N, A, LDA, IPIV, WORK, LWORK, INFO )
00179 *
00180 *  -- LAPACK computational routine (version 3.4.0) --
00181 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00182 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00183 *     November 2011
00184 *
00185 *     .. Scalar Arguments ..
00186       CHARACTER          UPLO
00187       INTEGER            INFO, LDA, LWORK, N
00188 *     ..
00189 *     .. Array Arguments ..
00190       INTEGER            IPIV( * )
00191       COMPLEX*16         A( LDA, * ), WORK( * )
00192 *     ..
00193 *
00194 *  =====================================================================
00195 *
00196 *     .. Local Scalars ..
00197       LOGICAL            LQUERY, UPPER
00198       INTEGER            IINFO, IWS, J, K, KB, LDWORK, LWKOPT, NB, NBMIN
00199 *     ..
00200 *     .. External Functions ..
00201       LOGICAL            LSAME
00202       INTEGER            ILAENV
00203       EXTERNAL           LSAME, ILAENV
00204 *     ..
00205 *     .. External Subroutines ..
00206       EXTERNAL           XERBLA, ZHETF2, ZLAHEF
00207 *     ..
00208 *     .. Intrinsic Functions ..
00209       INTRINSIC          MAX
00210 *     ..
00211 *     .. Executable Statements ..
00212 *
00213 *     Test the input parameters.
00214 *
00215       INFO = 0
00216       UPPER = LSAME( UPLO, 'U' )
00217       LQUERY = ( LWORK.EQ.-1 )
00218       IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
00219          INFO = -1
00220       ELSE IF( N.LT.0 ) THEN
00221          INFO = -2
00222       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00223          INFO = -4
00224       ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
00225          INFO = -7
00226       END IF
00227 *
00228       IF( INFO.EQ.0 ) THEN
00229 *
00230 *        Determine the block size
00231 *
00232          NB = ILAENV( 1, 'ZHETRF', UPLO, N, -1, -1, -1 )
00233          LWKOPT = N*NB
00234          WORK( 1 ) = LWKOPT
00235       END IF
00236 *
00237       IF( INFO.NE.0 ) THEN
00238          CALL XERBLA( 'ZHETRF', -INFO )
00239          RETURN
00240       ELSE IF( LQUERY ) THEN
00241          RETURN
00242       END IF
00243 *
00244       NBMIN = 2
00245       LDWORK = N
00246       IF( NB.GT.1 .AND. NB.LT.N ) THEN
00247          IWS = LDWORK*NB
00248          IF( LWORK.LT.IWS ) THEN
00249             NB = MAX( LWORK / LDWORK, 1 )
00250             NBMIN = MAX( 2, ILAENV( 2, 'ZHETRF', UPLO, N, -1, -1, -1 ) )
00251          END IF
00252       ELSE
00253          IWS = 1
00254       END IF
00255       IF( NB.LT.NBMIN )
00256      $   NB = N
00257 *
00258       IF( UPPER ) THEN
00259 *
00260 *        Factorize A as U*D*U**H using the upper triangle of A
00261 *
00262 *        K is the main loop index, decreasing from N to 1 in steps of
00263 *        KB, where KB is the number of columns factorized by ZLAHEF;
00264 *        KB is either NB or NB-1, or K for the last block
00265 *
00266          K = N
00267    10    CONTINUE
00268 *
00269 *        If K < 1, exit from loop
00270 *
00271          IF( K.LT.1 )
00272      $      GO TO 40
00273 *
00274          IF( K.GT.NB ) THEN
00275 *
00276 *           Factorize columns k-kb+1:k of A and use blocked code to
00277 *           update columns 1:k-kb
00278 *
00279             CALL ZLAHEF( UPLO, K, NB, KB, A, LDA, IPIV, WORK, N, IINFO )
00280          ELSE
00281 *
00282 *           Use unblocked code to factorize columns 1:k of A
00283 *
00284             CALL ZHETF2( UPLO, K, A, LDA, IPIV, IINFO )
00285             KB = K
00286          END IF
00287 *
00288 *        Set INFO on the first occurrence of a zero pivot
00289 *
00290          IF( INFO.EQ.0 .AND. IINFO.GT.0 )
00291      $      INFO = IINFO
00292 *
00293 *        Decrease K and return to the start of the main loop
00294 *
00295          K = K - KB
00296          GO TO 10
00297 *
00298       ELSE
00299 *
00300 *        Factorize A as L*D*L**H using the lower triangle of A
00301 *
00302 *        K is the main loop index, increasing from 1 to N in steps of
00303 *        KB, where KB is the number of columns factorized by ZLAHEF;
00304 *        KB is either NB or NB-1, or N-K+1 for the last block
00305 *
00306          K = 1
00307    20    CONTINUE
00308 *
00309 *        If K > N, exit from loop
00310 *
00311          IF( K.GT.N )
00312      $      GO TO 40
00313 *
00314          IF( K.LE.N-NB ) THEN
00315 *
00316 *           Factorize columns k:k+kb-1 of A and use blocked code to
00317 *           update columns k+kb:n
00318 *
00319             CALL ZLAHEF( UPLO, N-K+1, NB, KB, A( K, K ), LDA, IPIV( K ),
00320      $                   WORK, N, IINFO )
00321          ELSE
00322 *
00323 *           Use unblocked code to factorize columns k:n of A
00324 *
00325             CALL ZHETF2( UPLO, N-K+1, A( K, K ), LDA, IPIV( K ), IINFO )
00326             KB = N - K + 1
00327          END IF
00328 *
00329 *        Set INFO on the first occurrence of a zero pivot
00330 *
00331          IF( INFO.EQ.0 .AND. IINFO.GT.0 )
00332      $      INFO = IINFO + K - 1
00333 *
00334 *        Adjust IPIV
00335 *
00336          DO 30 J = K, K + KB - 1
00337             IF( IPIV( J ).GT.0 ) THEN
00338                IPIV( J ) = IPIV( J ) + K - 1
00339             ELSE
00340                IPIV( J ) = IPIV( J ) - K + 1
00341             END IF
00342    30    CONTINUE
00343 *
00344 *        Increase K and return to the start of the main loop
00345 *
00346          K = K + KB
00347          GO TO 20
00348 *
00349       END IF
00350 *
00351    40 CONTINUE
00352       WORK( 1 ) = LWKOPT
00353       RETURN
00354 *
00355 *     End of ZHETRF
00356 *
00357       END
 All Files Functions