SRC/zlahef_aasen.f

   1 *> \brief \b ZLAHEF_AASEN
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLAHEF_AASEN + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlahef_aasen.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlahef_aasen.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlahef_aasen.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLAHEF_AASEN( 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 *> DLATRF_AASEN factorizes a panel of a real 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 diatonals), 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_AASEN, 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 November 2016
 150 *
 151 *> \ingroup complex16SYcomputational
 152 *
 153 *  @precisions fortran z -> c
 154 *
 155 *  =====================================================================
 156       SUBROUTINE ZLAHEF_AASEN( UPLO, J1, M, NB, A, LDA, IPIV,
 157      $                         H, LDH, WORK, INFO )
 158 *
 159 *  -- LAPACK computational routine (version 3.4.0) --
 160 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 161 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 162 *     November 2016
 163 *
 164       IMPLICIT NONE
 165 *
 166 *     .. Scalar Arguments ..
 167       CHARACTER    UPLO
 168       INTEGER      M, NB, J1, LDA, LDH, INFO
 169 *     ..
 170 *     .. Array Arguments ..
 171       INTEGER      IPIV( * )
 172       COMPLEX*16   A( LDA, * ), H( LDH, * ), WORK( * )
 173 *     ..
 174 *
 175 *  =====================================================================
 176 *     .. Parameters ..
 177       COMPLEX*16   ZERO, ONE
 178       PARAMETER    ( ZERO = (0.0D+0, 0.0D+0), ONE = (1.0D+0, 0.0D+0) )
 179 *
 180 *     .. Local Scalars ..
 181       INTEGER      J, K, K1, I1, I2
 182       COMPLEX*16   PIV, ALPHA
 183 *     ..
 184 *     .. External Functions ..
 185       LOGICAL      LSAME
 186       INTEGER      IZAMAX, ILAENV
 187       EXTERNAL     LSAME, ILAENV, IZAMAX
 188 *     ..
 189 *     .. External Subroutines ..
 190       EXTERNAL     XERBLA
 191 *     ..
 192 *     .. Intrinsic Functions ..
 193       INTRINSIC    DBLE, DCONJG, MAX
 194 *     ..
 195 *     .. Executable Statements ..
 196 *
 197       INFO = 0
 198       J = 1
 199 *
 200 *     K1 is the first column of the panel to be factorized
 201 *     i.e.,  K1 is 2 for the first block column, and 1 for the rest of the blocks
 202 *
 203       K1 = (2-J1)+1
 204 *
 205       IF( LSAME( UPLO, 'U' ) ) THEN
 206 *
 207 *        .....................................................
 208 *        Factorize A as U**T*D*U using the upper triangle of A
 209 *        .....................................................
 210 *
 211  10      CONTINUE
 212          IF ( J.GT.MIN(M, NB) )
 213      $      GO TO 20
 214 *
 215 *        K is the column to be factorized
 216 *         when being called from ZHETRF_AASEN,
 217 *         > for the first block column, J1 is 1, hence J1+J-1 is J,
 218 *         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
 219 *
 220          K = J1+J-1
 221 *
 222 *        H(J:N, J) := A(J, J:N) - H(J:N, 1:(J-1)) * L(J1:(J-1), J),
 223 *         where H(J:N, J) has been initialized to be A(J, J:N)
 224 *
 225          IF( K.GT.2 ) THEN
 226 *
 227 *        K is the column to be factorized
 228 *         > for the first block column, K is J, skipping the first two
 229 *           columns
 230 *         > for the rest of the columns, K is J+1, skipping only the
 231 *           first column
 232 *
 233             CALL ZLACGV( J-K1, A( 1, J ), 1 )
 234             CALL ZGEMV( 'No transpose', M-J+1, J-K1,
 235      $                 -ONE, H( J, K1 ), LDH,
 236      $                       A( 1, J ), 1,
 237      $                  ONE, H( J, J ), 1 )
 238             CALL ZLACGV( J-K1, A( 1, J ), 1 )
 239          END IF
 240 *
 241 *        Copy H(i:n, i) into WORK
 242 *
 243          CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 244 *
 245          IF( J.GT.K1 ) THEN
 246 *
 247 *           Compute WORK := WORK - L(J-1, J:N) * T(J-1,J),
 248 *            where A(J-1, J) stores T(J-1, J) and A(J-2, J:N) stores U(J-1, J:N)
 249 *
 250             ALPHA = -DCONJG( A( K-1, J ) )
 251             CALL ZAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
 252          END IF
 253 *
 254 *        Set A(J, J) = T(J, J)
 255 *
 256          A( K, J ) = DBLE( WORK( 1 ) )
 257 *
 258          IF( J.LT.M ) THEN
 259 *
 260 *           Compute WORK(2:N) = T(J, J) L(J, (J+1):N)
 261 *            where A(J, J) stores T(J, J) and A(J-1, (J+1):N) stores U(J, (J+1):N)
 262 *
 263             IF( (J1+J-1).GT.1 ) THEN
 264                ALPHA = -A( K, J )
 265                CALL ZAXPY( M-J, ALPHA, A( K-1, J+1 ), LDA,
 266      $                                 WORK( 2 ), 1 )
 267             ENDIF
 268 *
 269 *           Find max(|WORK(2:n)|)
 270 *
 271             I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1
 272             PIV = WORK( I2 )
 273 *
 274 *           Apply hermitian pivot
 275 *
 276             IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
 277 *
 278 *              Swap WORK(I1) and WORK(I2)
 279 *
 280                I1 = 2
 281                WORK( I2 ) = WORK( I1 )
 282                WORK( I1 ) = PIV
 283 *
 284 *              Swap A(I1, I1+1:N) with A(I1+1:N, I2)
 285 *
 286                I1 = I1+J-1
 287                I2 = I2+J-1
 288                CALL ZSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
 289      $                              A( J1+I1, I2 ), 1 )
 290                CALL ZLACGV( I2-I1, A( J1+I1-1, I1+1 ), LDA )
 291                CALL ZLACGV( I2-I1-1, A( J1+I1, I2 ), 1 )
 292 *
 293 *              Swap A(I1, I2+1:N) with A(I2, I2+1:N)
 294 *
 295                CALL ZSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
 296      $                           A( J1+I2-1, I2+1 ), LDA )
 297 *
 298 *              Swap A(I1, I1) with A(I2,I2)
 299 *
 300                PIV = A( I1+J1-1, I1 )
 301                A( J1+I1-1, I1 ) = A( J1+I2-1, I2 )
 302                A( J1+I2-1, I2 ) = PIV
 303 *
 304 *              Swap H(I1, 1:J1) with H(I2, 1:J1)
 305 *
 306                CALL ZSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
 307                IPIV( I1 ) = I2
 308 *
 309                IF( I1.GT.(K1-1) ) THEN
 310 *
 311 *                 Swap L(1:I1-1, I1) with L(1:I1-1, I2),
 312 *                  skipping the first column
 313 *
 314                   CALL ZSWAP( I1-K1+1, A( 1, I1 ), 1,
 315      $                                 A( 1, I2 ), 1 )
 316                END IF
 317             ELSE
 318                IPIV( J+1 ) = J+1
 319             ENDIF
 320 *
 321 *           Set A(J, J+1) = T(J, J+1)
 322 *
 323             A( K, J+1 ) = WORK( 2 )
 324             IF( (A( K, J ).EQ.ZERO ) .AND.
 325      $        ( (J.EQ.M) .OR. (A( K, J+1 ).EQ.ZERO))) THEN
 326                 IF(INFO .EQ. 0) THEN
 327                     INFO = J
 328                 END IF
 329             END IF
 330 *
 331             IF( J.LT.NB ) THEN
 332 *
 333 *              Copy A(J+1:N, J+1) into H(J:N, J),
 334 *
 335                CALL ZCOPY( M-J, A( K+1, J+1 ), LDA,
 336      $                          H( J+1, J+1 ), 1 )
 337             END IF
 338 *
 339 *           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
 340 *            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
 341 *
 342             IF( A( K, J+1 ).NE.ZERO ) THEN
 343                ALPHA = ONE / A( K, J+1 )
 344                CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
 345                CALL ZSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
 346             ELSE
 347                CALL ZLASET( 'Full', 1, M-J-1, ZERO, ZERO,
 348      $                      A( K, J+2 ), LDA)
 349             END IF
 350          ELSE
 351             IF( (A( K, J ).EQ.ZERO) .AND. (INFO.EQ.0) ) THEN
 352                INFO = J
 353             END IF
 354          END IF
 355          J = J + 1
 356          GO TO 10
 357  20      CONTINUE
 358 *
 359       ELSE
 360 *
 361 *        .....................................................
 362 *        Factorize A as L*D*L**T using the lower triangle of A
 363 *        .....................................................
 364 *
 365  30      CONTINUE
 366          IF( J.GT.MIN( M, NB ) )
 367      $      GO TO 40
 368 *
 369 *        K is the column to be factorized
 370 *         when being called from ZHETRF_AASEN,
 371 *         > for the first block column, J1 is 1, hence J1+J-1 is J,
 372 *         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
 373 *
 374          K = J1+J-1
 375 *
 376 *        H(J:N, J) := A(J:N, J) - H(J:N, 1:(J-1)) * L(J, J1:(J-1))^T,
 377 *         where H(J:N, J) has been initialized to be A(J:N, J)
 378 *
 379          IF( K.GT.2 ) THEN
 380 *
 381 *        K is the column to be factorized
 382 *         > for the first block column, K is J, skipping the first two
 383 *           columns
 384 *         > for the rest of the columns, K is J+1, skipping only the
 385 *           first column
 386 *
 387             CALL ZLACGV( J-K1, A( J, 1 ), LDA )
 388             CALL ZGEMV( 'No transpose', M-J+1, J-K1,
 389      $                 -ONE, H( J, K1 ), LDH,
 390      $                       A( J, 1 ), LDA,
 391      $                  ONE, H( J, J ), 1 )
 392             CALL ZLACGV( J-K1, A( J, 1 ), LDA )
 393          END IF
 394 *
 395 *        Copy H(J:N, J) into WORK
 396 *
 397          CALL ZCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 398 *
 399          IF( J.GT.K1 ) THEN
 400 *
 401 *           Compute WORK := WORK - L(J:N, J-1) * T(J-1,J),
 402 *            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
 403 *
 404             ALPHA = -DCONJG( A( J, K-1 ) )
 405             CALL ZAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 )
 406          END IF
 407 *
 408 *        Set A(J, J) = T(J, J)
 409 *
 410          A( J, K ) = DBLE( WORK( 1 ) )
 411 *
 412          IF( J.LT.M ) THEN
 413 *
 414 *           Compute WORK(2:N) = T(J, J) L((J+1):N, J)
 415 *            where A(J, J) = T(J, J) and A((J+1):N, J-1) = L((J+1):N, J)
 416 *
 417             IF( (J1+J-1).GT.1 ) THEN
 418                ALPHA = -A( J, K )
 419                CALL ZAXPY( M-J, ALPHA, A( J+1, K-1 ), 1,
 420      $                                 WORK( 2 ), 1 )
 421             ENDIF
 422 *
 423 *           Find max(|WORK(2:n)|)
 424 *
 425             I2 = IZAMAX( M-J, WORK( 2 ), 1 ) + 1
 426             PIV = WORK( I2 )
 427 *
 428 *           Apply hermitian pivot
 429 *
 430             IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
 431 *
 432 *              Swap WORK(I1) and WORK(I2)
 433 *
 434                I1 = 2
 435                WORK( I2 ) = WORK( I1 )
 436                WORK( I1 ) = PIV
 437 *
 438 *              Swap A(I1+1:N, I1) with A(I2, I1+1:N)
 439 *
 440                I1 = I1+J-1
 441                I2 = I2+J-1
 442                CALL ZSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
 443      $                              A( I2, J1+I1 ), LDA )
 444                CALL ZLACGV( I2-I1, A( I1+1, J1+I1-1 ), 1 )
 445                CALL ZLACGV( I2-I1-1, A( I2, J1+I1 ), LDA )
 446 *
 447 *              Swap A(I2+1:N, I1) with A(I2+1:N, I2)
 448 *
 449                CALL ZSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
 450      $                           A( I2+1, J1+I2-1 ), 1 )
 451 *
 452 *              Swap A(I1, I1) with A(I2, I2)
 453 *
 454                PIV = A( I1, J1+I1-1 )
 455                A( I1, J1+I1-1 ) = A( I2, J1+I2-1 )
 456                A( I2, J1+I2-1 ) = PIV
 457 *
 458 *              Swap H(I1, I1:J1) with H(I2, I2:J1)
 459 *
 460                CALL ZSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
 461                IPIV( I1 ) = I2
 462 *
 463                IF( I1.GT.(K1-1) ) THEN
 464 *
 465 *                 Swap L(1:I1-1, I1) with L(1:I1-1, I2),
 466 *                  skipping the first column
 467 *
 468                   CALL ZSWAP( I1-K1+1, A( I1, 1 ), LDA,
 469      $                                 A( I2, 1 ), LDA )
 470                END IF
 471             ELSE
 472                IPIV( J+1 ) = J+1
 473             ENDIF
 474 *
 475 *           Set A(J+1, J) = T(J+1, J)
 476 *
 477             A( J+1, K ) = WORK( 2 )
 478             IF( (A( J, K ).EQ.ZERO) .AND.
 479      $        ( (J.EQ.M) .OR. (A( J+1, K ).EQ.ZERO)) ) THEN
 480                 IF (INFO .EQ. 0)
 481      $              INFO = J
 482             END IF
 483 *
 484             IF( J.LT.NB ) THEN
 485 *
 486 *              Copy A(J+1:N, J+1) into H(J+1:N, J),
 487 *
 488                CALL ZCOPY( M-J, A( J+1, K+1 ), 1,
 489      $                          H( J+1, J+1 ), 1 )
 490             END IF
 491 *
 492 *           Compute L(J+2, J+1) = WORK( 3:N ) / T(J, J+1),
 493 *            where A(J, J+1) = T(J, J+1) and A(J+2:N, J) = L(J+2:N, J+1)
 494 *
 495             IF( A( J+1, K ).NE.ZERO ) THEN
 496                ALPHA = ONE / A( J+1, K )
 497                CALL ZCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
 498                CALL ZSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
 499             ELSE
 500                CALL ZLASET( 'Full', M-J-1, 1, ZERO, ZERO,
 501      $                      A( J+2, K ), LDA )
 502             END IF
 503          ELSE
 504             IF( (A( J, K ).EQ.ZERO) .AND. (J.EQ.M)
 505      $          .AND. (INFO.EQ.0) ) INFO = J
 506          END IF
 507          J = J + 1
 508          GO TO 30
 509  40      CONTINUE
 510       END IF
 511       RETURN
 512 *
 513 *     End of ZLAHEF_AASEN
 514 *
 515       END