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 COMPLEX 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 COMPLEX 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 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*> \ingroup complexSYcomputational 140* 141* ===================================================================== 142 SUBROUTINE CLASYF_AA( UPLO, J1, M, NB, A, LDA, IPIV, 143 $ H, LDH, WORK ) 144* 145* -- LAPACK computational routine -- 146* -- LAPACK is a software package provided by Univ. of Tennessee, -- 147* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 148* 149 IMPLICIT NONE 150* 151* .. Scalar Arguments .. 152 CHARACTER UPLO 153 INTEGER M, NB, J1, LDA, LDH 154* .. 155* .. Array Arguments .. 156 INTEGER IPIV( * ) 157 COMPLEX A( LDA, * ), H( LDH, * ), WORK( * ) 158* .. 159* 160* ===================================================================== 161* .. Parameters .. 162 COMPLEX ZERO, ONE 163 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 164* 165* .. Local Scalars .. 166 INTEGER J, K, K1, I1, I2, MJ 167 COMPLEX PIV, ALPHA 168* .. 169* .. External Functions .. 170 LOGICAL LSAME 171 INTEGER ICAMAX, ILAENV 172 EXTERNAL LSAME, ILAENV, ICAMAX 173* .. 174* .. External Subroutines .. 175 EXTERNAL CAXPY, CGEMV, CSCAL, CCOPY, CSWAP, CLASET, 176 $ XERBLA 177* .. 178* .. Intrinsic Functions .. 179 INTRINSIC MAX 180* .. 181* .. Executable Statements .. 182* 183 J = 1 184* 185* K1 is the first column of the panel to be factorized 186* i.e., K1 is 2 for the first block column, and 1 for the rest of the blocks 187* 188 K1 = (2-J1)+1 189* 190 IF( LSAME( UPLO, 'U' ) ) THEN 191* 192* ..................................................... 193* Factorize A as U**T*D*U using the upper triangle of A 194* ..................................................... 195* 196 10 CONTINUE 197 IF ( J.GT.MIN(M, NB) ) 198 $ GO TO 20 199* 200* K is the column to be factorized 201* when being called from CSYTRF_AA, 202* > for the first block column, J1 is 1, hence J1+J-1 is J, 203* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, 204* 205 K = J1+J-1 206 IF( J.EQ.M ) THEN 207* 208* Only need to compute T(J, J) 209* 210 MJ = 1 211 ELSE 212 MJ = M-J+1 213 END IF 214* 215* H(J:M, J) := A(J, J:M) - H(J:M, 1:(J-1)) * L(J1:(J-1), J), 216* where H(J:M, J) has been initialized to be A(J, J:M) 217* 218 IF( K.GT.2 ) THEN 219* 220* K is the column to be factorized 221* > for the first block column, K is J, skipping the first two 222* columns 223* > for the rest of the columns, K is J+1, skipping only the 224* first column 225* 226 CALL CGEMV( 'No transpose', MJ, J-K1, 227 $ -ONE, H( J, K1 ), LDH, 228 $ A( 1, J ), 1, 229 $ ONE, H( J, J ), 1 ) 230 END IF 231* 232* Copy H(i:M, i) into WORK 233* 234 CALL CCOPY( MJ, H( J, J ), 1, WORK( 1 ), 1 ) 235* 236 IF( J.GT.K1 ) THEN 237* 238* Compute WORK := WORK - L(J-1, J:M) * T(J-1,J), 239* where A(J-1, J) stores T(J-1, J) and A(J-2, J:M) stores U(J-1, J:M) 240* 241 ALPHA = -A( K-1, J ) 242 CALL CAXPY( MJ, ALPHA, A( K-2, J ), LDA, WORK( 1 ), 1 ) 243 END IF 244* 245* Set A(J, J) = T(J, J) 246* 247 A( K, J ) = WORK( 1 ) 248* 249 IF( J.LT.M ) THEN 250* 251* Compute WORK(2:M) = T(J, J) L(J, (J+1):M) 252* where A(J, J) stores T(J, J) and A(J-1, (J+1):M) stores U(J, (J+1):M) 253* 254 IF( K.GT.1 ) THEN 255 ALPHA = -A( K, J ) 256 CALL CAXPY( M-J, ALPHA, A( K-1, J+1 ), LDA, 257 $ WORK( 2 ), 1 ) 258 ENDIF 259* 260* Find max(|WORK(2:M)|) 261* 262 I2 = ICAMAX( M-J, WORK( 2 ), 1 ) + 1 263 PIV = WORK( I2 ) 264* 265* Apply symmetric pivot 266* 267 IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN 268* 269* Swap WORK(I1) and WORK(I2) 270* 271 I1 = 2 272 WORK( I2 ) = WORK( I1 ) 273 WORK( I1 ) = PIV 274* 275* Swap A(I1, I1+1:M) with A(I1+1:M, I2) 276* 277 I1 = I1+J-1 278 I2 = I2+J-1 279 CALL CSWAP( I2-I1-1, A( J1+I1-1, I1+1 ), LDA, 280 $ A( J1+I1, I2 ), 1 ) 281* 282* Swap A(I1, I2+1:M) with A(I2, I2+1:M) 283* 284 IF( I2.LT.M ) 285 $ CALL CSWAP( M-I2, A( J1+I1-1, I2+1 ), LDA, 286 $ A( J1+I2-1, I2+1 ), LDA ) 287* 288* Swap A(I1, I1) with A(I2,I2) 289* 290 PIV = A( I1+J1-1, I1 ) 291 A( J1+I1-1, I1 ) = A( J1+I2-1, I2 ) 292 A( J1+I2-1, I2 ) = PIV 293* 294* Swap H(I1, 1:J1) with H(I2, 1:J1) 295* 296 CALL CSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH ) 297 IPIV( I1 ) = I2 298* 299 IF( I1.GT.(K1-1) ) THEN 300* 301* Swap L(1:I1-1, I1) with L(1:I1-1, I2), 302* skipping the first column 303* 304 CALL CSWAP( I1-K1+1, A( 1, I1 ), 1, 305 $ A( 1, I2 ), 1 ) 306 END IF 307 ELSE 308 IPIV( J+1 ) = J+1 309 ENDIF 310* 311* Set A(J, J+1) = T(J, J+1) 312* 313 A( K, J+1 ) = WORK( 2 ) 314* 315 IF( J.LT.NB ) THEN 316* 317* Copy A(J+1:M, J+1) into H(J:M, J), 318* 319 CALL CCOPY( M-J, A( K+1, J+1 ), LDA, 320 $ H( J+1, J+1 ), 1 ) 321 END IF 322* 323* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1), 324* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1) 325* 326 IF( J.LT.(M-1) ) THEN 327 IF( A( K, J+1 ).NE.ZERO ) THEN 328 ALPHA = ONE / A( K, J+1 ) 329 CALL CCOPY( M-J-1, WORK( 3 ), 1, A( K, J+2 ), LDA ) 330 CALL CSCAL( M-J-1, ALPHA, A( K, J+2 ), LDA ) 331 ELSE 332 CALL CLASET( 'Full', 1, M-J-1, ZERO, ZERO, 333 $ A( K, J+2 ), LDA) 334 END IF 335 END IF 336 END IF 337 J = J + 1 338 GO TO 10 339 20 CONTINUE 340* 341 ELSE 342* 343* ..................................................... 344* Factorize A as L*D*L**T using the lower triangle of A 345* ..................................................... 346* 347 30 CONTINUE 348 IF( J.GT.MIN( M, NB ) ) 349 $ GO TO 40 350* 351* K is the column to be factorized 352* when being called from CSYTRF_AA, 353* > for the first block column, J1 is 1, hence J1+J-1 is J, 354* > for the rest of the columns, J1 is 2, and J1+J-1 is J+1, 355* 356 K = J1+J-1 357 IF( J.EQ.M ) THEN 358* 359* Only need to compute T(J, J) 360* 361 MJ = 1 362 ELSE 363 MJ = M-J+1 364 END IF 365* 366* H(J:M, J) := A(J:M, J) - H(J:M, 1:(J-1)) * L(J, J1:(J-1))^T, 367* where H(J:M, J) has been initialized to be A(J:M, J) 368* 369 IF( K.GT.2 ) THEN 370* 371* K is the column to be factorized 372* > for the first block column, K is J, skipping the first two 373* columns 374* > for the rest of the columns, K is J+1, skipping only the 375* first column 376* 377 CALL CGEMV( 'No transpose', MJ, J-K1, 378 $ -ONE, H( J, K1 ), LDH, 379 $ A( J, 1 ), LDA, 380 $ ONE, H( J, J ), 1 ) 381 END IF 382* 383* Copy H(J:M, J) into WORK 384* 385 CALL CCOPY( MJ, H( J, J ), 1, WORK( 1 ), 1 ) 386* 387 IF( J.GT.K1 ) THEN 388* 389* Compute WORK := WORK - L(J:M, J-1) * T(J-1,J), 390* where A(J-1, J) = T(J-1, J) and A(J, J-2) = L(J, J-1) 391* 392 ALPHA = -A( J, K-1 ) 393 CALL CAXPY( MJ, ALPHA, A( J, K-2 ), 1, WORK( 1 ), 1 ) 394 END IF 395* 396* Set A(J, J) = T(J, J) 397* 398 A( J, K ) = WORK( 1 ) 399* 400 IF( J.LT.M ) THEN 401* 402* Compute WORK(2:M) = T(J, J) L((J+1):M, J) 403* where A(J, J) = T(J, J) and A((J+1):M, J-1) = L((J+1):M, J) 404* 405 IF( K.GT.1 ) THEN 406 ALPHA = -A( J, K ) 407 CALL CAXPY( M-J, ALPHA, A( J+1, K-1 ), 1, 408 $ WORK( 2 ), 1 ) 409 ENDIF 410* 411* Find max(|WORK(2:M)|) 412* 413 I2 = ICAMAX( M-J, WORK( 2 ), 1 ) + 1 414 PIV = WORK( I2 ) 415* 416* Apply symmetric pivot 417* 418 IF( (I2.NE.2) .AND. (PIV.NE.0) ) THEN 419* 420* Swap WORK(I1) and WORK(I2) 421* 422 I1 = 2 423 WORK( I2 ) = WORK( I1 ) 424 WORK( I1 ) = PIV 425* 426* Swap A(I1+1:M, I1) with A(I2, I1+1:M) 427* 428 I1 = I1+J-1 429 I2 = I2+J-1 430 CALL CSWAP( I2-I1-1, A( I1+1, J1+I1-1 ), 1, 431 $ A( I2, J1+I1 ), LDA ) 432* 433* Swap A(I2+1:M, I1) with A(I2+1:M, I2) 434* 435 IF( I2.LT.M ) 436 $ CALL CSWAP( M-I2, A( I2+1, J1+I1-1 ), 1, 437 $ A( I2+1, J1+I2-1 ), 1 ) 438* 439* Swap A(I1, I1) with A(I2, I2) 440* 441 PIV = A( I1, J1+I1-1 ) 442 A( I1, J1+I1-1 ) = A( I2, J1+I2-1 ) 443 A( I2, J1+I2-1 ) = PIV 444* 445* Swap H(I1, I1:J1) with H(I2, I2:J1) 446* 447 CALL CSWAP( I1-1, H( I1, 1 ), LDH, H( I2, 1 ), LDH ) 448 IPIV( I1 ) = I2 449* 450 IF( I1.GT.(K1-1) ) THEN 451* 452* Swap L(1:I1-1, I1) with L(1:I1-1, I2), 453* skipping the first column 454* 455 CALL CSWAP( I1-K1+1, A( I1, 1 ), LDA, 456 $ A( I2, 1 ), LDA ) 457 END IF 458 ELSE 459 IPIV( J+1 ) = J+1 460 ENDIF 461* 462* Set A(J+1, J) = T(J+1, J) 463* 464 A( J+1, K ) = WORK( 2 ) 465* 466 IF( J.LT.NB ) THEN 467* 468* Copy A(J+1:M, J+1) into H(J+1:M, J), 469* 470 CALL CCOPY( M-J, A( J+1, K+1 ), 1, 471 $ H( J+1, J+1 ), 1 ) 472 END IF 473* 474* Compute L(J+2, J+1) = WORK( 3:M ) / T(J, J+1), 475* where A(J, J+1) = T(J, J+1) and A(J+2:M, J) = L(J+2:M, J+1) 476* 477 IF( J.LT.(M-1) ) THEN 478 IF( A( J+1, K ).NE.ZERO ) THEN 479 ALPHA = ONE / A( J+1, K ) 480 CALL CCOPY( M-J-1, WORK( 3 ), 1, A( J+2, K ), 1 ) 481 CALL CSCAL( M-J-1, ALPHA, A( J+2, K ), 1 ) 482 ELSE 483 CALL CLASET( 'Full', M-J-1, 1, ZERO, ZERO, 484 $ A( J+2, K ), LDA ) 485 END IF 486 END IF 487 END IF 488 J = J + 1 489 GO TO 30 490 40 CONTINUE 491 END IF 492 RETURN 493* 494* End of CLASYF_AA 495* 496 END 497