SRC/clasyf_aa.f

   1 *> \brief \b CLASYF_AA
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download CLASYF_AA + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clasyf_aa.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clasyf_aa.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clasyf_aa.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE CLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
  22 *                             H, LDH, WORK )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          UPLO
  26 *       INTEGER            J1, M, NB, LDA, LDH
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IPIV( * )
  30 *       COMPLEX            A( LDA, * ), H( LDH, * ), WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DLATRF_AA factorizes a panel of a complex symmetric 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 CSYTRF_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 REAL 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,M).
 103 *> \endverbatim
 104 *>
 105 *> \param[out] IPIV
 106 *> \verbatim
 107 *>          IPIV is INTEGER array, dimension (M)
 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 REAL 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 REAL workspace, dimension (M).
 128 *> \endverbatim
 129 *>
 130 *
 131 *  Authors:
 132 *  ========
 133 *
 134 *> \author Univ. of Tennessee
 135 *> \author Univ. of California Berkeley
 136 *> \author Univ. of Colorado Denver
 137 *> \author NAG Ltd.
 138 *
 139 *> \date December 2016
 140 *
 141 *> \ingroup complexSYcomputational
 142 *
 143 *  =====================================================================
 144       SUBROUTINE CLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV,
 145      $                      H, LDH, WORK )
 146 *
 147 *  -- LAPACK computational routine (version 3.7.0) --
 148 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 149 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 150 *     December 2016
 151 *
 152       IMPLICIT NONE
 153 *
 154 *     .. Scalar Arguments ..
 155       CHARACTER          UPLO
 156       INTEGER            M, NB, J1, LDA, LDH
 157 *     ..
 158 *     .. Array Arguments ..
 159       INTEGER            IPIV( * )
 160       COMPLEX            A( LDA, * ), H( LDH, * ), WORK( * )
 161 *     ..
 162 *
 163 *  =====================================================================
 164 *     .. Parameters ..
 165       COMPLEX            ZERO, ONE
 166       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
 167 *
 168 *     .. Local Scalars ..
 169       INTEGER            J, K, K1, I1, I2
 170       COMPLEX            PIV, ALPHA
 171 *     ..
 172 *     .. External Functions ..
 173       LOGICAL            LSAME
 174       INTEGER            ICAMAX, ILAENV
 175       EXTERNAL           LSAME, ILAENV, ICAMAX
 176 *     ..
 177 *     .. External Subroutines ..
 178       EXTERNAL           XERBLA
 179 *     ..
 180 *     .. Intrinsic Functions ..
 181       INTRINSIC          MAX
 182 *     ..
 183 *     .. Executable Statements ..
 184 *
 185       J = 1
 186 *
 187 *     K1 is the first column of the panel to be factorized
 188 *     i.e.,  K1 is 2 for the first block column, and 1 for the rest of the blocks
 189 *
 190       K1 = (2-J1)+1
 191 *
 192       IF( LSAME( UPLO, 'U' ) ) THEN
 193 *
 194 *        .....................................................
 195 *        Factorize A as U**T*D*U using the upper triangle of A
 196 *        .....................................................
 197 *
 198  10      CONTINUE
 199          IF ( J.GT.MIN(M, NB) )
 200      $      GO TO 20
 201 *
 202 *        K is the column to be factorized
 203 *         when being called from CSYTRF_AA,
 204 *         > for the first block column, J1 is 1, hence J1+J-1 is J,
 205 *         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
 206 *
 207          K = J1+J-1
 208 *
 209 *        H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J),
 210 *         where H(J:M, J) has been initialized to be A(J, J:M)
 211 *
 212          IF( K.GT.2 ) THEN
 213 *
 214 *        K is the column to be factorized
 215 *         > for the first block column, K is J, skipping the first two
 216 *           columns
 217 *         > for the rest of the columns, K is J+1, skipping only the
 218 *           first column
 219 *
 220             CALL CGEMV( 'No transpose', M-J+1, J-K1,
 221      $                 -ONE, H( J, K1 ), LDH,
 222      $                       A( 1, J ), 1,
 223      $                  ONE, H( J, J ), 1 )
 224          END IF
 225 *
 226 *        Copy H(i:M, i) into WORK
 227 *
 228          CALL CCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 229 *
 230          IF( J.GT.K1 ) THEN
 231 *
 232 *           Compute WORK := WORK - L(J-1, J:M) * T(J-1,J),
 233 *            where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M)
 234 *
 235             ALPHA = -A( K-1, J )
 236             CALL CAXPY( M-J+1, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 )
 237          END IF
 238 *
 239 *        Set A(J, J) = T(J, J)
 240 *
 241          A( K, J ) = WORK( 1 )
 242 *
 243          IF( J.LT.M ) THEN
 244 *
 245 *           Compute WORK(2:M) = T(J, J) L(J, (J+1):M)
 246 *            where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M)
 247 *
 248             IF( K.GT.1 ) THEN
 249                ALPHA = -A( K, J )
 250                CALL CAXPY( M-J, ALPHA, A( K-1, J+1 ), LDA,
 251      $                                 WORK( 2 ), 1 )
 252             ENDIF
 253 *
 254 *           Find max(|WORK(2:M)|)
 255 *
 256             I2 = ICAMAX( M-J, WORK( 2 ), 1 ) + 1
 257             PIV = WORK( I2 )
 258 *
 259 *           Apply symmetric pivot
 260 *
 261             IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
 262 *
 263 *              Swap WORK(I1) and WORK(I2)
 264 *
 265                I1 = 2
 266                WORK( I2 ) = WORK( I1 )
 267                WORK( I1 ) = PIV
 268 *
 269 *              Swap A(I1, I1+1:M) with A(I1+1:M, I2)
 270 *
 271                I1 = I1+J-1
 272                I2 = I2+J-1
 273                CALL CSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA,
 274      $                              A( J1+I1, I2 ), 1 )
 275 *
 276 *              Swap A(I1, I2+1:M) with A(I2, I2+1:M)
 277 *
 278                CALL CSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA,
 279      $                           A( J1+I2-1, I2+1 ), LDA )
 280 *
 281 *              Swap A(I1, I1) with A(I2,I2)
 282 *
 283                PIV = A( I1+J1-1, I1 )
 284                A( J1+I1-1, I1 ) = A( J1+I2-1, I2 )
 285                A( J1+I2-1, I2 ) = PIV
 286 *
 287 *              Swap H(I1, 1:J1) with H(I2, 1:J1)
 288 *
 289                CALL CSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
 290                IPIV( I1 ) = I2
 291 *
 292                IF( I1.GT.(K1-1) ) THEN
 293 *
 294 *                 Swap L(1:I1-1, I1) with L(1:I1-1, I2),
 295 *                  skipping the first column
 296 *
 297                   CALL CSWAP( I1-K1+1, A( 1, I1 ), 1,
 298      $                                 A( 1, I2 ), 1 )
 299                END IF
 300             ELSE
 301                IPIV( J+1 ) = J+1
 302             ENDIF
 303 *
 304 *           Set A(J, J+1) = T(J, J+1)
 305 *
 306             A( K, J+1 ) = WORK( 2 )
 307 *
 308             IF( J.LT.NB ) THEN
 309 *
 310 *              Copy A(J+1:M, J+1) into H(J:M, J),
 311 *
 312                CALL CCOPY( M-J, A( K+1, J+1 ), LDA,
 313      $                          H( J+1, J+1 ), 1 )
 314             END IF
 315 *
 316 *           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
 317 *            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 318 *
 319             IF( A( K, J+1 ).NE.ZERO ) THEN
 320                ALPHA = ONE / A( K, J+1 )
 321                CALL CCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA )
 322                CALL CSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA )
 323             ELSE
 324                CALL CLASET( 'Full', 1, M-J-1, ZERO, ZERO,
 325      $                      A( K, J+2 ), LDA)
 326             END IF
 327          END IF
 328          J = J + 1
 329          GO TO 10
 330  20      CONTINUE
 331 *
 332       ELSE
 333 *
 334 *        .....................................................
 335 *        Factorize A as L*D*L**T using the lower triangle of A
 336 *        .....................................................
 337 *
 338  30      CONTINUE
 339          IF( J.GT.MIN( M, NB ) )
 340      $      GO TO 40
 341 *
 342 *        K is the column to be factorized
 343 *         when being called from CSYTRF_AA,
 344 *         > for the first block column, J1 is 1, hence J1+J-1 is J,
 345 *         > for the rest of the columns, J1 is 2, and J1+J-1 is J+1,
 346 *
 347          K = J1+J-1
 348 *
 349 *        H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T,
 350 *         where H(J:M, J) has been initialized to be A(J:M, J)
 351 *
 352          IF( K.GT.2 ) THEN
 353 *
 354 *        K is the column to be factorized
 355 *         > for the first block column, K is J, skipping the first two
 356 *           columns
 357 *         > for the rest of the columns, K is J+1, skipping only the
 358 *           first column
 359 *
 360             CALL CGEMV( 'No transpose', M-J+1, J-K1,
 361      $                 -ONE, H( J, K1 ), LDH,
 362      $                       A( J, 1 ), LDA,
 363      $                  ONE, H( J, J ), 1 )
 364          END IF
 365 *
 366 *        Copy H(J:M, J) into WORK
 367 *
 368          CALL CCOPY( M-J+1, H( J, J ), 1, WORK( 1 ), 1 )
 369 *
 370          IF( J.GT.K1 ) THEN
 371 *
 372 *           Compute WORK := WORK - L(J:M, J-1) * T(J-1,J),
 373 *            where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1)
 374 *
 375             ALPHA = -A( J, K-1 )
 376             CALL CAXPY( M-J+1, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 )
 377          END IF
 378 *
 379 *        Set A(J, J) = T(J, J)
 380 *
 381          A( J, K ) = WORK( 1 )
 382 *
 383          IF( J.LT.M ) THEN
 384 *
 385 *           Compute WORK(2:M) = T(J, J) L((J+1):M, J)
 386 *            where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J)
 387 *
 388             IF( K.GT.1 ) THEN
 389                ALPHA = -A( J, K )
 390                CALL CAXPY( M-J, ALPHA, A( J+1, K-1 ), 1,
 391      $                                 WORK( 2 ), 1 )
 392             ENDIF
 393 *
 394 *           Find max(|WORK(2:M)|)
 395 *
 396             I2 = ICAMAX( M-J, WORK( 2 ), 1 ) + 1
 397             PIV = WORK( I2 )
 398 *
 399 *           Apply symmetric pivot
 400 *
 401             IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN
 402 *
 403 *              Swap WORK(I1) and WORK(I2)
 404 *
 405                I1 = 2
 406                WORK( I2 ) = WORK( I1 )
 407                WORK( I1 ) = PIV
 408 *
 409 *              Swap A(I1+1:M, I1) with A(I2, I1+1:M)
 410 *
 411                I1 = I1+J-1
 412                I2 = I2+J-1
 413                CALL CSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1,
 414      $                              A( I2, J1+I1 ), LDA )
 415 *
 416 *              Swap A(I2+1:M, I1) with A(I2+1:M, I2)
 417 *
 418                CALL CSWAP( M-I2, A( I2+1, J1+I1-1 ), 1,
 419      $                           A( I2+1, J1+I2-1 ), 1 )
 420 *
 421 *              Swap A(I1, I1) with A(I2, I2)
 422 *
 423                PIV = A( I1, J1+I1-1 )
 424                A( I1, J1+I1-1 ) = A( I2, J1+I2-1 )
 425                A( I2, J1+I2-1 ) = PIV
 426 *
 427 *              Swap H(I1, I1:J1) with H(I2, I2:J1)
 428 *
 429                CALL CSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH )
 430                IPIV( I1 ) = I2
 431 *
 432                IF( I1.GT.(K1-1) ) THEN
 433 *
 434 *                 Swap L(1:I1-1, I1) with L(1:I1-1, I2),
 435 *                  skipping the first column
 436 *
 437                   CALL CSWAP( I1-K1+1, A( I1, 1 ), LDA,
 438      $                                 A( I2, 1 ), LDA )
 439                END IF
 440             ELSE
 441                IPIV( J+1 ) = J+1
 442             ENDIF
 443 *
 444 *           Set A(J+1, J) = T(J+1, J)
 445 *
 446             A( J+1, K ) = WORK( 2 )
 447 *
 448             IF( J.LT.NB ) THEN
 449 *
 450 *              Copy A(J+1:M, J+1) into H(J+1:M, J),
 451 *
 452                CALL CCOPY( M-J, A( J+1, K+1 ), 1,
 453      $                          H( J+1, J+1 ), 1 )
 454             END IF
 455 *
 456 *           Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1),
 457 *            where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1)
 458 *
 459             IF( A( J+1, K ).NE.ZERO ) THEN
 460                ALPHA = ONE / A( J+1, K )
 461                CALL CCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 )
 462                CALL CSCAL( M-J-1, ALPHA, A( J+2, K ), 1 )
 463             ELSE
 464                CALL CLASET( 'Full', M-J-1, 1, ZERO, ZERO,
 465      $                      A( J+2, K ), LDA )
 466             END IF
 467          END IF
 468          J = J + 1
 469          GO TO 30
 470  40      CONTINUE
 471       END IF
 472       RETURN
 473 *
 474 *     End of CLASYF_AA
 475 *
 476       END