Annotation of rpl/lapack/lapack/zlahef_aa.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZLAHEF_AA
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZLAHEF_AA + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahef_aa.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahef_aa.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahef_aa.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
! 22: * H, LDH, WORK, INFO )
! 23: *
! 24: * .. Scalar Arguments ..
! 25: * CHARACTER UPLO
! 26: * INTEGER J1, M, NB, LDA, LDH, INFO
! 27: * ..
! 28: * .. Array Arguments ..
! 29: * INTEGER IPIV( * )
! 30: * COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * )
! 31: * ..
! 32: *
! 33: *
! 34: *> \par Purpose:
! 35: * =============
! 36: *>
! 37: *> \verbatim
! 38: *>
! 39: *> DLAHEF_AA factorizes a panel of a complex hermitian matrix A using
! 40: *> the Aasen's algorithm. The panel consists of a set of NB rows of A
! 41: *> when UPLO is U, or a set of NB columns when UPLO is L.
! 42: *>
! 43: *> In order to factorize the panel, the Aasen's algorithm requires the
! 44: *> last row, or column, of the previous panel. The first row, or column,
! 45: *> of A is set to be the first row, or column, of an identity matrix,
! 46: *> which is used to factorize the first panel.
! 47: *>
! 48: *> The resulting J-th row of U, or J-th column of L, is stored in the
! 49: *> (J-1)-th row, or column, of A (without the unit diagonals), while
! 50: *> the diagonal and subdiagonal of A are overwritten by those of T.
! 51: *>
! 52: *> \endverbatim
! 53: *
! 54: * Arguments:
! 55: * ==========
! 56: *
! 57: *> \param[in] UPLO
! 58: *> \verbatim
! 59: *> UPLO is CHARACTER*1
! 60: *> = 'U': Upper triangle of A is stored;
! 61: *> = 'L': Lower triangle of A is stored.
! 62: *> \endverbatim
! 63: *>
! 64: *> \param[in] J1
! 65: *> \verbatim
! 66: *> J1 is INTEGER
! 67: *> The location of the first row, or column, of the panel
! 68: *> within the submatrix of A, passed to this routine, e.g.,
! 69: *> when called by ZHETRF_AA, for the first panel, J1 is 1,
! 70: *> while for the remaining panels, J1 is 2.
! 71: *> \endverbatim
! 72: *>
! 73: *> \param[in] M
! 74: *> \verbatim
! 75: *> M is INTEGER
! 76: *> The dimension of the submatrix. M >= 0.
! 77: *> \endverbatim
! 78: *>
! 79: *> \param[in] NB
! 80: *> \verbatim
! 81: *> NB is INTEGER
! 82: *> The dimension of the panel to be facotorized.
! 83: *> \endverbatim
! 84: *>
! 85: *> \param[in,out] A
! 86: *> \verbatim
! 87: *> A is COMPLEX*16 array, dimension (LDA,M) for
! 88: *> the first panel, while dimension (LDA,M+1) for the
! 89: *> remaining panels.
! 90: *>
! 91: *> On entry, A contains the last row, or column, of
! 92: *> the previous panel, and the trailing submatrix of A
! 93: *> to be factorized, except for the first panel, only
! 94: *> the panel is passed.
! 95: *>
! 96: *> On exit, the leading panel is factorized.
! 97: *> \endverbatim
! 98: *>
! 99: *> \param[in] LDA
! 100: *> \verbatim
! 101: *> LDA is INTEGER
! 102: *> The leading dimension of the array A. LDA >= max(1,N).
! 103: *> \endverbatim
! 104: *>
! 105: *> \param[out] IPIV
! 106: *> \verbatim
! 107: *> IPIV is INTEGER array, dimension (N)
! 108: *> Details of the row and column interchanges,
! 109: *> the row and column k were interchanged with the row and
! 110: *> column IPIV(k).
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[in,out] H
! 114: *> \verbatim
! 115: *> H is COMPLEX*16 workspace, dimension (LDH,NB).
! 116: *>
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[in] LDH
! 120: *> \verbatim
! 121: *> LDH is INTEGER
! 122: *> The leading dimension of the workspace H. LDH >= max(1,M).
! 123: *> \endverbatim
! 124: *>
! 125: *> \param[out] WORK
! 126: *> \verbatim
! 127: *> WORK is COMPLEX*16 workspace, dimension (M).
! 128: *> \endverbatim
! 129: *>
! 130: *> \param[out] INFO
! 131: *> \verbatim
! 132: *> INFO is INTEGER
! 133: *> = 0: successful exit
! 134: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 135: *> > 0: if INFO = i, D(i,i) is exactly zero. The factorization
! 136: *> has been completed, but the block diagonal matrix D is
! 137: *> exactly singular, and division by zero will occur if it
! 138: *> is used to solve a system of equations.
! 139: *> \endverbatim
! 140: *
! 141: * Authors:
! 142: * ========
! 143: *
! 144: *> \author Univ. of Tennessee
! 145: *> \author Univ. of California Berkeley
! 146: *> \author Univ. of Colorado Denver
! 147: *> \author NAG Ltd.
! 148: *
! 149: *> \date December 2016
! 150: *
! 151: *> \ingroup complex16HEcomputational
! 152: *
! 153: * =====================================================================
! 154: SUBROUTINE ZLAHEF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
! 155: $ H, LDH, WORK, INFO )
! 156: *
! 157: * -- LAPACK computational routine (version 3.7.0) --
! 158: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 159: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 160: * December 2016
! 161: *
! 162: IMPLICIT NONE
! 163: *
! 164: * .. Scalar Arguments ..
! 165: CHARACTER UPLO
! 166: INTEGER M, NB, J1, LDA, LDH, INFO
! 167: * ..
! 168: * .. Array Arguments ..
! 169: INTEGER IPIV( * )
! 170: COMPLEX*16 A( LDA, * ), H( LDH, * ), WORK( * )
! 171: * ..
! 172: *
! 173: * =====================================================================
! 174: * .. Parameters ..
! 175: COMPLEX*16 ZERO, ONE
! 176: PARAMETER ( ZERO = (0.0D+0, 0.0D+0), ONE = (1.0D+0, 0.0D+0) )
! 177: *
! 178: * .. Local Scalars ..
! 179: INTEGER J, K, K1, I1, I2
! 180: COMPLEX*16 PIV, ALPHA
! 181: * ..
! 182: * .. External Functions ..
! 183: LOGICAL LSAME
! 184: INTEGER IZAMAX, ILAENV
! 185: EXTERNAL LSAME, ILAENV, IZAMAX
! 186: * ..
! 187: * .. External Subroutines ..
! 188: EXTERNAL XERBLA
! 189: * ..
! 190: * .. Intrinsic Functions ..
! 191: INTRINSIC DBLE, DCONJG, MAX
! 192: * ..
! 193: * .. Executable Statements ..
! 194: *
! 195: INFO = 0
! 196: J = 1
! 197: *
! 198: * K1 is the first column of the panel to be factorized
! 199: * i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks
! 200: *
! 201: K1 = (2-J1)+1
! 202: *
! 203: IF( LSAME( UPLO, 'U' ) ) THEN
! 204: *
! 205: * .....................................................
! 206: * Factorize A as U**T*D*U using the upper triangle of A
! 207: * .....................................................
! 208: *
! 209: 10 CONTINUE
! 210: IF ( J.GT.MIN(M, NB) )
! 211: $ GO TO 20
! 212: *
! 213: * K is the column to be factorized
! 214: * when being called from ZHETRF_AA,
! 215: * > for the first block column, J1 is 1, hence J1+J-1 is J,
! 216: * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
! 217: *
! 218: K = J1+J-1
! 219: *
! 220: * H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
! 221: * where H(J:N, J) has been initialized to be A(J, J:N)
! 222: *
! 223: IF( K.GT.2 ) THEN
! 224: *
! 225: * K is the column to be factorized
! 226: * > for the first block column, K is J, skipping the first two
! 227: * columns
! 228: * > for the rest of the columns, K is J+1, skipping only the
! 229: * first column
! 230: *
! 231: CALL ZLACGV( J-K1, A( 1, J ), 1 )
! 232: CALL ZGEMV( 'No transpose', M-J+1, J-K1,
! 233: $ -ONE, H( J, K1 ), LDH,
! 234: $ A( 1, J ), 1,
! 235: $ ONE, H( J, J ), 1 )
! 236: CALL ZLACGV( J-K1, A( 1, J ), 1 )
! 237: END IF
! 238: *
! 239: * Copy H(i:n, i) into WORK
! 240: *
! 241: CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
! 242: *
! 243: IF( J.GT.K1 ) THEN
! 244: *
! 245: * Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
! 246: * where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
! 247: *
! 248: ALPHA = -DCONJG( A( K-1, J ) )
! 249: CALL ZAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
! 250: END IF
! 251: *
! 252: * Set A(J, J) = T(J, J)
! 253: *
! 254: A( K, J ) = DBLE( WORK( 1 ) )
! 255: *
! 256: IF( J.LT.M ) THEN
! 257: *
! 258: * Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
! 259: * where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
! 260: *
! 261: IF( K.GT.1 ) THEN
! 262: ALPHA = -A( K, J )
! 263: CALL ZAXPY( M-J, ALPHA, A( K-1, J+1 ), LDA,
! 264: $ WORK( 2 ), 1 )
! 265: ENDIF
! 266: *
! 267: * Find max(|WORK(2:n)|)
! 268: *
! 269: I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1
! 270: PIV = WORK( I2 )
! 271: *
! 272: * Apply hermitian pivot
! 273: *
! 274: IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
! 275: *
! 276: * Swap WORK(I1) and WORK(I2)
! 277: *
! 278: I1 = 2
! 279: WORK( I2 ) = WORK( I1 )
! 280: WORK( I1 ) = PIV
! 281: *
! 282: * Swap A(I1, I1+1:N) with A(I1+1:N, I2)
! 283: *
! 284: I1 = I1+J-1
! 285: I2 = I2+J-1
! 286: CALL ZSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
! 287: $ A( J1+I1, I2 ), 1 )
! 288: CALL ZLACGV( I2-I1, A( J1+I1-1, I1+1 ), LDA )
! 289: CALL ZLACGV( I2-I1-1, A( J1+I1, I2 ), 1 )
! 290: *
! 291: * Swap A(I1, I2+1:N) with A(I2, I2+1:N)
! 292: *
! 293: CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
! 294: $ A( J1+I2-1, I2+1 ), LDA )
! 295: *
! 296: * Swap A(I1, I1) with A(I2,I2)
! 297: *
! 298: PIV = A( I1+J1-1, I1 )
! 299: A( J1+I1-1, I1 ) = A( J1+I2-1, I2 )
! 300: A( J1+I2-1, I2 ) = PIV
! 301: *
! 302: * Swap H(I1, 1:J1) with H(I2, 1:J1)
! 303: *
! 304: CALL ZSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
! 305: IPIV( I1 ) = I2
! 306: *
! 307: IF( I1.GT.(K1-1) ) THEN
! 308: *
! 309: * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
! 310: * skipping the first column
! 311: *
! 312: CALL ZSWAP( I1-K1+1, A( 1, I1 ), 1,
! 313: $ A( 1, I2 ), 1 )
! 314: END IF
! 315: ELSE
! 316: IPIV( J+1 ) = J+1
! 317: ENDIF
! 318: *
! 319: * Set A(J, J+1) = T(J, J+1)
! 320: *
! 321: A( K, J+1 ) = WORK( 2 )
! 322: IF( (A( K, J ).EQ.ZERO ) .AND.
! 323: $ ( (J.EQ.M) .OR. (A( K, J+1 ).EQ.ZERO))) THEN
! 324: IF(INFO .EQ. 0) THEN
! 325: INFO = J
! 326: END IF
! 327: END IF
! 328: *
! 329: IF( J.LT.NB ) THEN
! 330: *
! 331: * Copy A(J+1:N, J+1) into H(J:N, J),
! 332: *
! 333: CALL ZCOPY( M-J, A( K+1, J+1 ), LDA,
! 334: $ H( J+1, J+1 ), 1 )
! 335: END IF
! 336: *
! 337: * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
! 338: * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
! 339: *
! 340: IF( A( K, J+1 ).NE.ZERO ) THEN
! 341: ALPHA = ONE / A( K, J+1 )
! 342: CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
! 343: CALL ZSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
! 344: ELSE
! 345: CALL ZLASET( 'Full', 1, M-J-1, ZERO, ZERO,
! 346: $ A( K, J+2 ), LDA)
! 347: END IF
! 348: ELSE
! 349: IF( (A( K, J ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN
! 350: INFO = J
! 351: END IF
! 352: END IF
! 353: J = J + 1
! 354: GO TO 10
! 355: 20 CONTINUE
! 356: *
! 357: ELSE
! 358: *
! 359: * .....................................................
! 360: * Factorize A as L*D*L**T using the lower triangle of A
! 361: * .....................................................
! 362: *
! 363: 30 CONTINUE
! 364: IF( J.GT.MIN( M, NB ) )
! 365: $ GO TO 40
! 366: *
! 367: * K is the column to be factorized
! 368: * when being called from ZHETRF_AA,
! 369: * > for the first block column, J1 is 1, hence J1+J-1 is J,
! 370: * > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
! 371: *
! 372: K = J1+J-1
! 373: *
! 374: * H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
! 375: * where H(J:N, J) has been initialized to be A(J:N, J)
! 376: *
! 377: IF( K.GT.2 ) THEN
! 378: *
! 379: * K is the column to be factorized
! 380: * > for the first block column, K is J, skipping the first two
! 381: * columns
! 382: * > for the rest of the columns, K is J+1, skipping only the
! 383: * first column
! 384: *
! 385: CALL ZLACGV( J-K1, A( J, 1 ), LDA )
! 386: CALL ZGEMV( 'No transpose', M-J+1, J-K1,
! 387: $ -ONE, H( J, K1 ), LDH,
! 388: $ A( J, 1 ), LDA,
! 389: $ ONE, H( J, J ), 1 )
! 390: CALL ZLACGV( J-K1, A( J, 1 ), LDA )
! 391: END IF
! 392: *
! 393: * Copy H(J:N, J) into WORK
! 394: *
! 395: CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
! 396: *
! 397: IF( J.GT.K1 ) THEN
! 398: *
! 399: * Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
! 400: * where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
! 401: *
! 402: ALPHA = -DCONJG( A( J, K-1 ) )
! 403: CALL ZAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 )
! 404: END IF
! 405: *
! 406: * Set A(J, J) = T(J, J)
! 407: *
! 408: A( J, K ) = DBLE( WORK( 1 ) )
! 409: *
! 410: IF( J.LT.M ) THEN
! 411: *
! 412: * Compute WORK(2:N) = T(J, J) L((J+1):N, J)
! 413: * where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
! 414: *
! 415: IF( K.GT.1 ) THEN
! 416: ALPHA = -A( J, K )
! 417: CALL ZAXPY( M-J, ALPHA, A( J+1, K-1 ), 1,
! 418: $ WORK( 2 ), 1 )
! 419: ENDIF
! 420: *
! 421: * Find max(|WORK(2:n)|)
! 422: *
! 423: I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1
! 424: PIV = WORK( I2 )
! 425: *
! 426: * Apply hermitian pivot
! 427: *
! 428: IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
! 429: *
! 430: * Swap WORK(I1) and WORK(I2)
! 431: *
! 432: I1 = 2
! 433: WORK( I2 ) = WORK( I1 )
! 434: WORK( I1 ) = PIV
! 435: *
! 436: * Swap A(I1+1:N, I1) with A(I2, I1+1:N)
! 437: *
! 438: I1 = I1+J-1
! 439: I2 = I2+J-1
! 440: CALL ZSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
! 441: $ A( I2, J1+I1 ), LDA )
! 442: CALL ZLACGV( I2-I1, A( I1+1, J1+I1-1 ), 1 )
! 443: CALL ZLACGV( I2-I1-1, A( I2, J1+I1 ), LDA )
! 444: *
! 445: * Swap A(I2+1:N, I1) with A(I2+1:N, I2)
! 446: *
! 447: CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
! 448: $ A( I2+1, J1+I2-1 ), 1 )
! 449: *
! 450: * Swap A(I1, I1) with A(I2, I2)
! 451: *
! 452: PIV = A( I1, J1+I1-1 )
! 453: A( I1, J1+I1-1 ) = A( I2, J1+I2-1 )
! 454: A( I2, J1+I2-1 ) = PIV
! 455: *
! 456: * Swap H(I1, I1:J1) with H(I2, I2:J1)
! 457: *
! 458: CALL ZSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
! 459: IPIV( I1 ) = I2
! 460: *
! 461: IF( I1.GT.(K1-1) ) THEN
! 462: *
! 463: * Swap L(1:I1-1, I1) with L(1:I1-1, I2),
! 464: * skipping the first column
! 465: *
! 466: CALL ZSWAP( I1-K1+1, A( I1, 1 ), LDA,
! 467: $ A( I2, 1 ), LDA )
! 468: END IF
! 469: ELSE
! 470: IPIV( J+1 ) = J+1
! 471: ENDIF
! 472: *
! 473: * Set A(J+1, J) = T(J+1, J)
! 474: *
! 475: A( J+1, K ) = WORK( 2 )
! 476: IF( (A( J, K ).EQ.ZERO) .AND.
! 477: $ ( (J.EQ.M) .OR. (A( J+1, K ).EQ.ZERO)) ) THEN
! 478: IF (INFO .EQ. 0)
! 479: $ INFO = J
! 480: END IF
! 481: *
! 482: IF( J.LT.NB ) THEN
! 483: *
! 484: * Copy A(J+1:N, J+1) into H(J+1:N, J),
! 485: *
! 486: CALL ZCOPY( M-J, A( J+1, K+1 ), 1,
! 487: $ H( J+1, J+1 ), 1 )
! 488: END IF
! 489: *
! 490: * Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
! 491: * where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
! 492: *
! 493: IF( A( J+1, K ).NE.ZERO ) THEN
! 494: ALPHA = ONE / A( J+1, K )
! 495: CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
! 496: CALL ZSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
! 497: ELSE
! 498: CALL ZLASET( 'Full', M-J-1, 1, ZERO, ZERO,
! 499: $ A( J+2, K ), LDA )
! 500: END IF
! 501: ELSE
! 502: IF( (A( J, K ).EQ.ZERO) .AND. (J.EQ.M)
! 503: $ .AND. (INFO.EQ.0) ) INFO = J
! 504: END IF
! 505: J = J + 1
! 506: GO TO 30
! 507: 40 CONTINUE
! 508: END IF
! 509: RETURN
! 510: *
! 511: * End of ZLAHEF_AA
! 512: *
! 513: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlahef_aa.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack