1*> \brief \b ZHETRF_RK computes the factorization of a complex Hermitian indefinite matrix using the bounded Bunch-Kaufman (rook) diagonal pivoting method (BLAS3 blocked algorithm). 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download ZHETRF_RK + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zhetrf_rk.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zhetrf_rk.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zhetrf_rk.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE ZHETRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK, 22* INFO ) 23* 24* .. Scalar Arguments .. 25* CHARACTER UPLO 26* INTEGER INFO, LDA, LWORK, N 27* .. 28* .. Array Arguments .. 29* INTEGER IPIV( * ) 30* COMPLEX*16 A( LDA, * ), E ( * ), WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> ZHETRF_RK computes the factorization of a complex Hermitian matrix A 39*> using the bounded Bunch-Kaufman (rook) diagonal pivoting method: 40*> 41*> A = P*U*D*(U**H)*(P**T) or A = P*L*D*(L**H)*(P**T), 42*> 43*> where U (or L) is unit upper (or lower) triangular matrix, 44*> U**H (or L**H) is the conjugate of U (or L), P is a permutation 45*> matrix, P**T is the transpose of P, and D is Hermitian and block 46*> diagonal with 1-by-1 and 2-by-2 diagonal blocks. 47*> 48*> This is the blocked version of the algorithm, calling Level 3 BLAS. 49*> For more information see Further Details section. 50*> \endverbatim 51* 52* Arguments: 53* ========== 54* 55*> \param[in] UPLO 56*> \verbatim 57*> UPLO is CHARACTER*1 58*> Specifies whether the upper or lower triangular part of the 59*> Hermitian matrix A is stored: 60*> = 'U': Upper triangular 61*> = 'L': Lower triangular 62*> \endverbatim 63*> 64*> \param[in] N 65*> \verbatim 66*> N is INTEGER 67*> The order of the matrix A. N >= 0. 68*> \endverbatim 69*> 70*> \param[in,out] A 71*> \verbatim 72*> A is COMPLEX*16 array, dimension (LDA,N) 73*> On entry, the Hermitian matrix A. 74*> If UPLO = 'U': the leading N-by-N upper triangular part 75*> of A contains the upper triangular part of the matrix A, 76*> and the strictly lower triangular part of A is not 77*> referenced. 78*> 79*> If UPLO = 'L': the leading N-by-N lower triangular part 80*> of A contains the lower triangular part of the matrix A, 81*> and the strictly upper triangular part of A is not 82*> referenced. 83*> 84*> On exit, contains: 85*> a) ONLY diagonal elements of the Hermitian block diagonal 86*> matrix D on the diagonal of A, i.e. D(k,k) = A(k,k); 87*> (superdiagonal (or subdiagonal) elements of D 88*> are stored on exit in array E), and 89*> b) If UPLO = 'U': factor U in the superdiagonal part of A. 90*> If UPLO = 'L': factor L in the subdiagonal part of A. 91*> \endverbatim 92*> 93*> \param[in] LDA 94*> \verbatim 95*> LDA is INTEGER 96*> The leading dimension of the array A. LDA >= max(1,N). 97*> \endverbatim 98*> 99*> \param[out] E 100*> \verbatim 101*> E is COMPLEX*16 array, dimension (N) 102*> On exit, contains the superdiagonal (or subdiagonal) 103*> elements of the Hermitian block diagonal matrix D 104*> with 1-by-1 or 2-by-2 diagonal blocks, where 105*> If UPLO = 'U': E(i) = D(i-1,i), i=2:N, E(1) is set to 0; 106*> If UPLO = 'L': E(i) = D(i+1,i), i=1:N-1, E(N) is set to 0. 107*> 108*> NOTE: For 1-by-1 diagonal block D(k), where 109*> 1 <= k <= N, the element E(k) is set to 0 in both 110*> UPLO = 'U' or UPLO = 'L' cases. 111*> \endverbatim 112*> 113*> \param[out] IPIV 114*> \verbatim 115*> IPIV is INTEGER array, dimension (N) 116*> IPIV describes the permutation matrix P in the factorization 117*> of matrix A as follows. The absolute value of IPIV(k) 118*> represents the index of row and column that were 119*> interchanged with the k-th row and column. The value of UPLO 120*> describes the order in which the interchanges were applied. 121*> Also, the sign of IPIV represents the block structure of 122*> the Hermitian block diagonal matrix D with 1-by-1 or 2-by-2 123*> diagonal blocks which correspond to 1 or 2 interchanges 124*> at each factorization step. For more info see Further 125*> Details section. 126*> 127*> If UPLO = 'U', 128*> ( in factorization order, k decreases from N to 1 ): 129*> a) A single positive entry IPIV(k) > 0 means: 130*> D(k,k) is a 1-by-1 diagonal block. 131*> If IPIV(k) != k, rows and columns k and IPIV(k) were 132*> interchanged in the matrix A(1:N,1:N); 133*> If IPIV(k) = k, no interchange occurred. 134*> 135*> b) A pair of consecutive negative entries 136*> IPIV(k) < 0 and IPIV(k-1) < 0 means: 137*> D(k-1:k,k-1:k) is a 2-by-2 diagonal block. 138*> (NOTE: negative entries in IPIV appear ONLY in pairs). 139*> 1) If -IPIV(k) != k, rows and columns 140*> k and -IPIV(k) were interchanged 141*> in the matrix A(1:N,1:N). 142*> If -IPIV(k) = k, no interchange occurred. 143*> 2) If -IPIV(k-1) != k-1, rows and columns 144*> k-1 and -IPIV(k-1) were interchanged 145*> in the matrix A(1:N,1:N). 146*> If -IPIV(k-1) = k-1, no interchange occurred. 147*> 148*> c) In both cases a) and b), always ABS( IPIV(k) ) <= k. 149*> 150*> d) NOTE: Any entry IPIV(k) is always NONZERO on output. 151*> 152*> If UPLO = 'L', 153*> ( in factorization order, k increases from 1 to N ): 154*> a) A single positive entry IPIV(k) > 0 means: 155*> D(k,k) is a 1-by-1 diagonal block. 156*> If IPIV(k) != k, rows and columns k and IPIV(k) were 157*> interchanged in the matrix A(1:N,1:N). 158*> If IPIV(k) = k, no interchange occurred. 159*> 160*> b) A pair of consecutive negative entries 161*> IPIV(k) < 0 and IPIV(k+1) < 0 means: 162*> D(k:k+1,k:k+1) is a 2-by-2 diagonal block. 163*> (NOTE: negative entries in IPIV appear ONLY in pairs). 164*> 1) If -IPIV(k) != k, rows and columns 165*> k and -IPIV(k) were interchanged 166*> in the matrix A(1:N,1:N). 167*> If -IPIV(k) = k, no interchange occurred. 168*> 2) If -IPIV(k+1) != k+1, rows and columns 169*> k-1 and -IPIV(k-1) were interchanged 170*> in the matrix A(1:N,1:N). 171*> If -IPIV(k+1) = k+1, no interchange occurred. 172*> 173*> c) In both cases a) and b), always ABS( IPIV(k) ) >= k. 174*> 175*> d) NOTE: Any entry IPIV(k) is always NONZERO on output. 176*> \endverbatim 177*> 178*> \param[out] WORK 179*> \verbatim 180*> WORK is COMPLEX*16 array, dimension ( MAX(1,LWORK) ). 181*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 182*> \endverbatim 183*> 184*> \param[in] LWORK 185*> \verbatim 186*> LWORK is INTEGER 187*> The length of WORK. LWORK >=1. For best performance 188*> LWORK >= N*NB, where NB is the block size returned 189*> by ILAENV. 190*> 191*> If LWORK = -1, then a workspace query is assumed; 192*> the routine only calculates the optimal size of the WORK 193*> array, returns this value as the first entry of the WORK 194*> array, and no error message related to LWORK is issued 195*> by XERBLA. 196*> \endverbatim 197*> 198*> \param[out] INFO 199*> \verbatim 200*> INFO is INTEGER 201*> = 0: successful exit 202*> 203*> < 0: If INFO = -k, the k-th argument had an illegal value 204*> 205*> > 0: If INFO = k, the matrix A is singular, because: 206*> If UPLO = 'U': column k in the upper 207*> triangular part of A contains all zeros. 208*> If UPLO = 'L': column k in the lower 209*> triangular part of A contains all zeros. 210*> 211*> Therefore D(k,k) is exactly zero, and superdiagonal 212*> elements of column k of U (or subdiagonal elements of 213*> column k of L ) are all zeros. The factorization has 214*> been completed, but the block diagonal matrix D is 215*> exactly singular, and division by zero will occur if 216*> it is used to solve a system of equations. 217*> 218*> NOTE: INFO only stores the first occurrence of 219*> a singularity, any subsequent occurrence of singularity 220*> is not stored in INFO even though the factorization 221*> always completes. 222*> \endverbatim 223* 224* Authors: 225* ======== 226* 227*> \author Univ. of Tennessee 228*> \author Univ. of California Berkeley 229*> \author Univ. of Colorado Denver 230*> \author NAG Ltd. 231* 232*> \ingroup complex16HEcomputational 233* 234*> \par Further Details: 235* ===================== 236*> 237*> \verbatim 238*> TODO: put correct description 239*> \endverbatim 240* 241*> \par Contributors: 242* ================== 243*> 244*> \verbatim 245*> 246*> December 2016, Igor Kozachenko, 247*> Computer Science Division, 248*> University of California, Berkeley 249*> 250*> September 2007, Sven Hammarling, Nicholas J. Higham, Craig Lucas, 251*> School of Mathematics, 252*> University of Manchester 253*> 254*> \endverbatim 255* 256* ===================================================================== 257 SUBROUTINE ZHETRF_RK( UPLO, N, A, LDA, E, IPIV, WORK, LWORK, 258 $ INFO ) 259* 260* -- LAPACK computational routine -- 261* -- LAPACK is a software package provided by Univ. of Tennessee, -- 262* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 263* 264* .. Scalar Arguments .. 265 CHARACTER UPLO 266 INTEGER INFO, LDA, LWORK, N 267* .. 268* .. Array Arguments .. 269 INTEGER IPIV( * ) 270 COMPLEX*16 A( LDA, * ), E( * ), WORK( * ) 271* .. 272* 273* ===================================================================== 274* 275* .. Local Scalars .. 276 LOGICAL LQUERY, UPPER 277 INTEGER I, IINFO, IP, IWS, K, KB, LDWORK, LWKOPT, 278 $ NB, NBMIN 279* .. 280* .. External Functions .. 281 LOGICAL LSAME 282 INTEGER ILAENV 283 EXTERNAL LSAME, ILAENV 284* .. 285* .. External Subroutines .. 286 EXTERNAL ZLAHEF_RK, ZHETF2_RK, ZSWAP, XERBLA 287* .. 288* .. Intrinsic Functions .. 289 INTRINSIC ABS, MAX 290* .. 291* .. Executable Statements .. 292* 293* Test the input parameters. 294* 295 INFO = 0 296 UPPER = LSAME( UPLO, 'U' ) 297 LQUERY = ( LWORK.EQ.-1 ) 298 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN 299 INFO = -1 300 ELSE IF( N.LT.0 ) THEN 301 INFO = -2 302 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 303 INFO = -4 304 ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN 305 INFO = -8 306 END IF 307* 308 IF( INFO.EQ.0 ) THEN 309* 310* Determine the block size 311* 312 NB = ILAENV( 1, 'ZHETRF_RK', UPLO, N, -1, -1, -1 ) 313 LWKOPT = N*NB 314 WORK( 1 ) = LWKOPT 315 END IF 316* 317 IF( INFO.NE.0 ) THEN 318 CALL XERBLA( 'ZHETRF_RK', -INFO ) 319 RETURN 320 ELSE IF( LQUERY ) THEN 321 RETURN 322 END IF 323* 324 NBMIN = 2 325 LDWORK = N 326 IF( NB.GT.1 .AND. NB.LT.N ) THEN 327 IWS = LDWORK*NB 328 IF( LWORK.LT.IWS ) THEN 329 NB = MAX( LWORK / LDWORK, 1 ) 330 NBMIN = MAX( 2, ILAENV( 2, 'ZHETRF_RK', 331 $ UPLO, N, -1, -1, -1 ) ) 332 END IF 333 ELSE 334 IWS = 1 335 END IF 336 IF( NB.LT.NBMIN ) 337 $ NB = N 338* 339 IF( UPPER ) THEN 340* 341* Factorize A as U*D*U**T using the upper triangle of A 342* 343* K is the main loop index, decreasing from N to 1 in steps of 344* KB, where KB is the number of columns factorized by ZLAHEF_RK; 345* KB is either NB or NB-1, or K for the last block 346* 347 K = N 348 10 CONTINUE 349* 350* If K < 1, exit from loop 351* 352 IF( K.LT.1 ) 353 $ GO TO 15 354* 355 IF( K.GT.NB ) THEN 356* 357* Factorize columns k-kb+1:k of A and use blocked code to 358* update columns 1:k-kb 359* 360 CALL ZLAHEF_RK( UPLO, K, NB, KB, A, LDA, E, 361 $ IPIV, WORK, LDWORK, IINFO ) 362 ELSE 363* 364* Use unblocked code to factorize columns 1:k of A 365* 366 CALL ZHETF2_RK( UPLO, K, A, LDA, E, IPIV, IINFO ) 367 KB = K 368 END IF 369* 370* Set INFO on the first occurrence of a zero pivot 371* 372 IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 373 $ INFO = IINFO 374* 375* No need to adjust IPIV 376* 377* 378* Apply permutations to the leading panel 1:k-1 379* 380* Read IPIV from the last block factored, i.e. 381* indices k-kb+1:k and apply row permutations to the 382* last k+1 colunms k+1:N after that block 383* (We can do the simple loop over IPIV with decrement -1, 384* since the ABS value of IPIV( I ) represents the row index 385* of the interchange with row i in both 1x1 and 2x2 pivot cases) 386* 387 IF( K.LT.N ) THEN 388 DO I = K, ( K - KB + 1 ), -1 389 IP = ABS( IPIV( I ) ) 390 IF( IP.NE.I ) THEN 391 CALL ZSWAP( N-K, A( I, K+1 ), LDA, 392 $ A( IP, K+1 ), LDA ) 393 END IF 394 END DO 395 END IF 396* 397* Decrease K and return to the start of the main loop 398* 399 K = K - KB 400 GO TO 10 401* 402* This label is the exit from main loop over K decreasing 403* from N to 1 in steps of KB 404* 405 15 CONTINUE 406* 407 ELSE 408* 409* Factorize A as L*D*L**T using the lower triangle of A 410* 411* K is the main loop index, increasing from 1 to N in steps of 412* KB, where KB is the number of columns factorized by ZLAHEF_RK; 413* KB is either NB or NB-1, or N-K+1 for the last block 414* 415 K = 1 416 20 CONTINUE 417* 418* If K > N, exit from loop 419* 420 IF( K.GT.N ) 421 $ GO TO 35 422* 423 IF( K.LE.N-NB ) THEN 424* 425* Factorize columns k:k+kb-1 of A and use blocked code to 426* update columns k+kb:n 427* 428 CALL ZLAHEF_RK( UPLO, N-K+1, NB, KB, A( K, K ), LDA, E( K ), 429 $ IPIV( K ), WORK, LDWORK, IINFO ) 430 431 432 ELSE 433* 434* Use unblocked code to factorize columns k:n of A 435* 436 CALL ZHETF2_RK( UPLO, N-K+1, A( K, K ), LDA, E( K ), 437 $ IPIV( K ), IINFO ) 438 KB = N - K + 1 439* 440 END IF 441* 442* Set INFO on the first occurrence of a zero pivot 443* 444 IF( INFO.EQ.0 .AND. IINFO.GT.0 ) 445 $ INFO = IINFO + K - 1 446* 447* Adjust IPIV 448* 449 DO I = K, K + KB - 1 450 IF( IPIV( I ).GT.0 ) THEN 451 IPIV( I ) = IPIV( I ) + K - 1 452 ELSE 453 IPIV( I ) = IPIV( I ) - K + 1 454 END IF 455 END DO 456* 457* Apply permutations to the leading panel 1:k-1 458* 459* Read IPIV from the last block factored, i.e. 460* indices k:k+kb-1 and apply row permutations to the 461* first k-1 colunms 1:k-1 before that block 462* (We can do the simple loop over IPIV with increment 1, 463* since the ABS value of IPIV( I ) represents the row index 464* of the interchange with row i in both 1x1 and 2x2 pivot cases) 465* 466 IF( K.GT.1 ) THEN 467 DO I = K, ( K + KB - 1 ), 1 468 IP = ABS( IPIV( I ) ) 469 IF( IP.NE.I ) THEN 470 CALL ZSWAP( K-1, A( I, 1 ), LDA, 471 $ A( IP, 1 ), LDA ) 472 END IF 473 END DO 474 END IF 475* 476* Increase K and return to the start of the main loop 477* 478 K = K + KB 479 GO TO 20 480* 481* This label is the exit from main loop over K increasing 482* from 1 to N in steps of KB 483* 484 35 CONTINUE 485* 486* End Lower 487* 488 END IF 489* 490 WORK( 1 ) = LWKOPT 491 RETURN 492* 493* End of ZHETRF_RK 494* 495 END 496