![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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