Annotation of rpl/lapack/lapack/zlaqz2.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZLAQZ2
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download ZLAQZ2 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ZLAQZ2.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ZLAQZ2.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ZLAQZ2.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE ZLAQZ2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B,
! 22: * $ LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
! 23: * $ WORK, LWORK, RWORK, REC, INFO )
! 24: * IMPLICIT NONE
! 25: *
! 26: * Arguments
! 27: * LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
! 28: * INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
! 29: * $ LDQC, LDZC, LWORK, REC
! 30: *
! 31: * COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
! 32: * $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
! 33: * INTEGER, INTENT( OUT ) :: NS, ND, INFO
! 34: * COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
! 35: * DOUBLE PRECISION :: RWORK( * )
! 36: * ..
! 37: *
! 38: *
! 39: *> \par Purpose:
! 40: * =============
! 41: *>
! 42: *> \verbatim
! 43: *>
! 44: *> ZLAQZ2 performs AED
! 45: *> \endverbatim
! 46: *
! 47: * Arguments:
! 48: * ==========
! 49: *
! 50: *> \param[in] ILSCHUR
! 51: *> \verbatim
! 52: *> ILSCHUR is LOGICAL
! 53: *> Determines whether or not to update the full Schur form
! 54: *> \endverbatim
! 55: *> \param[in] ILQ
! 56: *> \verbatim
! 57: *> ILQ is LOGICAL
! 58: *> Determines whether or not to update the matrix Q
! 59: *> \endverbatim
! 60: *>
! 61: *> \param[in] ILZ
! 62: *> \verbatim
! 63: *> ILZ is LOGICAL
! 64: *> Determines whether or not to update the matrix Z
! 65: *> \endverbatim
! 66: *>
! 67: *> \param[in] N
! 68: *> \verbatim
! 69: *> N is INTEGER
! 70: *> The order of the matrices A, B, Q, and Z. N >= 0.
! 71: *> \endverbatim
! 72: *>
! 73: *> \param[in] ILO
! 74: *> \verbatim
! 75: *> ILO is INTEGER
! 76: *> \endverbatim
! 77: *>
! 78: *> \param[in] IHI
! 79: *> \verbatim
! 80: *> IHI is INTEGER
! 81: *> ILO and IHI mark the rows and columns of (A,B) which
! 82: *> are to be normalized
! 83: *> \endverbatim
! 84: *>
! 85: *> \param[in] NW
! 86: *> \verbatim
! 87: *> NW is INTEGER
! 88: *> The desired size of the deflation window.
! 89: *> \endverbatim
! 90: *>
! 91: *> \param[in,out] A
! 92: *> \verbatim
! 93: *> A is COMPLEX*16 array, dimension (LDA, N)
! 94: *> \endverbatim
! 95: *>
! 96: *> \param[in] LDA
! 97: *> \verbatim
! 98: *> LDA is INTEGER
! 99: *> The leading dimension of the array A. LDA >= max( 1, N ).
! 100: *> \endverbatim
! 101: *>
! 102: *> \param[in,out] B
! 103: *> \verbatim
! 104: *> B is COMPLEX*16 array, dimension (LDB, N)
! 105: *> \endverbatim
! 106: *>
! 107: *> \param[in] LDB
! 108: *> \verbatim
! 109: *> LDB is INTEGER
! 110: *> The leading dimension of the array B. LDB >= max( 1, N ).
! 111: *> \endverbatim
! 112: *>
! 113: *> \param[in,out] Q
! 114: *> \verbatim
! 115: *> Q is COMPLEX*16 array, dimension (LDQ, N)
! 116: *> \endverbatim
! 117: *>
! 118: *> \param[in] LDQ
! 119: *> \verbatim
! 120: *> LDQ is INTEGER
! 121: *> \endverbatim
! 122: *>
! 123: *> \param[in,out] Z
! 124: *> \verbatim
! 125: *> Z is COMPLEX*16 array, dimension (LDZ, N)
! 126: *> \endverbatim
! 127: *>
! 128: *> \param[in] LDZ
! 129: *> \verbatim
! 130: *> LDZ is INTEGER
! 131: *> \endverbatim
! 132: *>
! 133: *> \param[out] NS
! 134: *> \verbatim
! 135: *> NS is INTEGER
! 136: *> The number of unconverged eigenvalues available to
! 137: *> use as shifts.
! 138: *> \endverbatim
! 139: *>
! 140: *> \param[out] ND
! 141: *> \verbatim
! 142: *> ND is INTEGER
! 143: *> The number of converged eigenvalues found.
! 144: *> \endverbatim
! 145: *>
! 146: *> \param[out] ALPHA
! 147: *> \verbatim
! 148: *> ALPHA is COMPLEX*16 array, dimension (N)
! 149: *> Each scalar alpha defining an eigenvalue
! 150: *> of GNEP.
! 151: *> \endverbatim
! 152: *>
! 153: *> \param[out] BETA
! 154: *> \verbatim
! 155: *> BETA is COMPLEX*16 array, dimension (N)
! 156: *> The scalars beta that define the eigenvalues of GNEP.
! 157: *> Together, the quantities alpha = ALPHA(j) and
! 158: *> beta = BETA(j) represent the j-th eigenvalue of the matrix
! 159: *> pair (A,B), in one of the forms lambda = alpha/beta or
! 160: *> mu = beta/alpha. Since either lambda or mu may overflow,
! 161: *> they should not, in general, be computed.
! 162: *> \endverbatim
! 163: *>
! 164: *> \param[in,out] QC
! 165: *> \verbatim
! 166: *> QC is COMPLEX*16 array, dimension (LDQC, NW)
! 167: *> \endverbatim
! 168: *>
! 169: *> \param[in] LDQC
! 170: *> \verbatim
! 171: *> LDQC is INTEGER
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[in,out] ZC
! 175: *> \verbatim
! 176: *> ZC is COMPLEX*16 array, dimension (LDZC, NW)
! 177: *> \endverbatim
! 178: *>
! 179: *> \param[in] LDZC
! 180: *> \verbatim
! 181: *> LDZ is INTEGER
! 182: *> \endverbatim
! 183: *>
! 184: *> \param[out] WORK
! 185: *> \verbatim
! 186: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
! 187: *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
! 188: *> \endverbatim
! 189: *>
! 190: *> \param[in] LWORK
! 191: *> \verbatim
! 192: *> LWORK is INTEGER
! 193: *> The dimension of the array WORK. LWORK >= max(1,N).
! 194: *>
! 195: *> If LWORK = -1, then a workspace query is assumed; the routine
! 196: *> only calculates the optimal size of the WORK array, returns
! 197: *> this value as the first entry of the WORK array, and no error
! 198: *> message related to LWORK is issued by XERBLA.
! 199: *> \endverbatim
! 200: *>
! 201: *> \param[out] RWORK
! 202: *> \verbatim
! 203: *> RWORK is DOUBLE PRECISION array, dimension (N)
! 204: *> \endverbatim
! 205: *>
! 206: *> \param[in] REC
! 207: *> \verbatim
! 208: *> REC is INTEGER
! 209: *> REC indicates the current recursion level. Should be set
! 210: *> to 0 on first call.
! 211: *> \endverbatim
! 212: *>
! 213: *> \param[out] INFO
! 214: *> \verbatim
! 215: *> INFO is INTEGER
! 216: *> = 0: successful exit
! 217: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 218: *> \endverbatim
! 219: *
! 220: * Authors:
! 221: * ========
! 222: *
! 223: *> \author Thijs Steel, KU Leuven
! 224: *
! 225: *> \date May 2020
! 226: *
! 227: *> \ingroup complex16GEcomputational
! 228: *>
! 229: * =====================================================================
! 230: RECURSIVE SUBROUTINE ZLAQZ2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
! 231: $ A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
! 232: $ ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
! 233: $ WORK, LWORK, RWORK, REC, INFO )
! 234: IMPLICIT NONE
! 235:
! 236: * Arguments
! 237: LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
! 238: INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
! 239: $ LDQC, LDZC, LWORK, REC
! 240:
! 241: COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
! 242: $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
! 243: INTEGER, INTENT( OUT ) :: NS, ND, INFO
! 244: COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
! 245: DOUBLE PRECISION :: RWORK( * )
! 246:
! 247: * Parameters
! 248: COMPLEX*16 CZERO, CONE
! 249: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), CONE = ( 1.0D+0,
! 250: $ 0.0D+0 ) )
! 251: DOUBLE PRECISION :: ZERO, ONE, HALF
! 252: PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 )
! 253:
! 254: * Local Scalars
! 255: INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, ZTGEXC_INFO,
! 256: $ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
! 257: DOUBLE PRECISION ::SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
! 258: COMPLEX*16 :: S, S1, TEMP
! 259:
! 260: * External Functions
! 261: EXTERNAL :: XERBLA, ZLAQZ0, ZLAQZ1, DLABAD, ZLACPY, ZLASET, ZGEMM,
! 262: $ ZTGEXC, ZLARTG, ZROT
! 263: DOUBLE PRECISION, EXTERNAL :: DLAMCH
! 264:
! 265: INFO = 0
! 266:
! 267: * Set up deflation window
! 268: JW = MIN( NW, IHI-ILO+1 )
! 269: KWTOP = IHI-JW+1
! 270: IF ( KWTOP .EQ. ILO ) THEN
! 271: S = CZERO
! 272: ELSE
! 273: S = A( KWTOP, KWTOP-1 )
! 274: END IF
! 275:
! 276: * Determine required workspace
! 277: IFST = 1
! 278: ILST = JW
! 279: CALL ZLAQZ0( 'S', 'V', 'V', JW, 1, JW, A( KWTOP, KWTOP ), LDA,
! 280: $ B( KWTOP, KWTOP ), LDB, ALPHA, BETA, QC, LDQC, ZC,
! 281: $ LDZC, WORK, -1, RWORK, REC+1, QZ_SMALL_INFO )
! 282: LWORKREQ = INT( WORK( 1 ) )+2*JW**2
! 283: LWORKREQ = MAX( LWORKREQ, N*NW, 2*NW**2+N )
! 284: IF ( LWORK .EQ.-1 ) THEN
! 285: * workspace query, quick return
! 286: WORK( 1 ) = LWORKREQ
! 287: RETURN
! 288: ELSE IF ( LWORK .LT. LWORKREQ ) THEN
! 289: INFO = -26
! 290: END IF
! 291:
! 292: IF( INFO.NE.0 ) THEN
! 293: CALL XERBLA( 'ZLAQZ2', -INFO )
! 294: RETURN
! 295: END IF
! 296:
! 297: * Get machine constants
! 298: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
! 299: SAFMAX = ONE/SAFMIN
! 300: CALL DLABAD( SAFMIN, SAFMAX )
! 301: ULP = DLAMCH( 'PRECISION' )
! 302: SMLNUM = SAFMIN*( DBLE( N )/ULP )
! 303:
! 304: IF ( IHI .EQ. KWTOP ) THEN
! 305: * 1 by 1 deflation window, just try a regular deflation
! 306: ALPHA( KWTOP ) = A( KWTOP, KWTOP )
! 307: BETA( KWTOP ) = B( KWTOP, KWTOP )
! 308: NS = 1
! 309: ND = 0
! 310: IF ( ABS( S ) .LE. MAX( SMLNUM, ULP*ABS( A( KWTOP,
! 311: $ KWTOP ) ) ) ) THEN
! 312: NS = 0
! 313: ND = 1
! 314: IF ( KWTOP .GT. ILO ) THEN
! 315: A( KWTOP, KWTOP-1 ) = CZERO
! 316: END IF
! 317: END IF
! 318: END IF
! 319:
! 320:
! 321: * Store window in case of convergence failure
! 322: CALL ZLACPY( 'ALL', JW, JW, A( KWTOP, KWTOP ), LDA, WORK, JW )
! 323: CALL ZLACPY( 'ALL', JW, JW, B( KWTOP, KWTOP ), LDB, WORK( JW**2+
! 324: $ 1 ), JW )
! 325:
! 326: * Transform window to real schur form
! 327: CALL ZLASET( 'FULL', JW, JW, CZERO, CONE, QC, LDQC )
! 328: CALL ZLASET( 'FULL', JW, JW, CZERO, CONE, ZC, LDZC )
! 329: CALL ZLAQZ0( 'S', 'V', 'V', JW, 1, JW, A( KWTOP, KWTOP ), LDA,
! 330: $ B( KWTOP, KWTOP ), LDB, ALPHA, BETA, QC, LDQC, ZC,
! 331: $ LDZC, WORK( 2*JW**2+1 ), LWORK-2*JW**2, RWORK,
! 332: $ REC+1, QZ_SMALL_INFO )
! 333:
! 334: IF( QZ_SMALL_INFO .NE. 0 ) THEN
! 335: * Convergence failure, restore the window and exit
! 336: ND = 0
! 337: NS = JW-QZ_SMALL_INFO
! 338: CALL ZLACPY( 'ALL', JW, JW, WORK, JW, A( KWTOP, KWTOP ), LDA )
! 339: CALL ZLACPY( 'ALL', JW, JW, WORK( JW**2+1 ), JW, B( KWTOP,
! 340: $ KWTOP ), LDB )
! 341: RETURN
! 342: END IF
! 343:
! 344: * Deflation detection loop
! 345: IF ( KWTOP .EQ. ILO .OR. S .EQ. CZERO ) THEN
! 346: KWBOT = KWTOP-1
! 347: ELSE
! 348: KWBOT = IHI
! 349: K = 1
! 350: K2 = 1
! 351: DO WHILE ( K .LE. JW )
! 352: * Try to deflate eigenvalue
! 353: TEMPR = ABS( A( KWBOT, KWBOT ) )
! 354: IF( TEMPR .EQ. ZERO ) THEN
! 355: TEMPR = ABS( S )
! 356: END IF
! 357: IF ( ( ABS( S*QC( 1, KWBOT-KWTOP+1 ) ) ) .LE. MAX( ULP*
! 358: $ TEMPR, SMLNUM ) ) THEN
! 359: * Deflatable
! 360: KWBOT = KWBOT-1
! 361: ELSE
! 362: * Not deflatable, move out of the way
! 363: IFST = KWBOT-KWTOP+1
! 364: ILST = K2
! 365: CALL ZTGEXC( .TRUE., .TRUE., JW, A( KWTOP, KWTOP ),
! 366: $ LDA, B( KWTOP, KWTOP ), LDB, QC, LDQC,
! 367: $ ZC, LDZC, IFST, ILST, ZTGEXC_INFO )
! 368: K2 = K2+1
! 369: END IF
! 370:
! 371: K = K+1
! 372: END DO
! 373: END IF
! 374:
! 375: * Store eigenvalues
! 376: ND = IHI-KWBOT
! 377: NS = JW-ND
! 378: K = KWTOP
! 379: DO WHILE ( K .LE. IHI )
! 380: ALPHA( K ) = A( K, K )
! 381: BETA( K ) = B( K, K )
! 382: K = K+1
! 383: END DO
! 384:
! 385: IF ( KWTOP .NE. ILO .AND. S .NE. CZERO ) THEN
! 386: * Reflect spike back, this will create optimally packed bulges
! 387: A( KWTOP:KWBOT, KWTOP-1 ) = A( KWTOP, KWTOP-1 ) *DCONJG( QC( 1,
! 388: $ 1:JW-ND ) )
! 389: DO K = KWBOT-1, KWTOP, -1
! 390: CALL ZLARTG( A( K, KWTOP-1 ), A( K+1, KWTOP-1 ), C1, S1,
! 391: $ TEMP )
! 392: A( K, KWTOP-1 ) = TEMP
! 393: A( K+1, KWTOP-1 ) = CZERO
! 394: K2 = MAX( KWTOP, K-1 )
! 395: CALL ZROT( IHI-K2+1, A( K, K2 ), LDA, A( K+1, K2 ), LDA, C1,
! 396: $ S1 )
! 397: CALL ZROT( IHI-( K-1 )+1, B( K, K-1 ), LDB, B( K+1, K-1 ),
! 398: $ LDB, C1, S1 )
! 399: CALL ZROT( JW, QC( 1, K-KWTOP+1 ), 1, QC( 1, K+1-KWTOP+1 ),
! 400: $ 1, C1, DCONJG( S1 ) )
! 401: END DO
! 402:
! 403: * Chase bulges down
! 404: ISTARTM = KWTOP
! 405: ISTOPM = IHI
! 406: K = KWBOT-1
! 407: DO WHILE ( K .GE. KWTOP )
! 408:
! 409: * Move bulge down and remove it
! 410: DO K2 = K, KWBOT-1
! 411: CALL ZLAQZ1( .TRUE., .TRUE., K2, KWTOP, KWTOP+JW-1,
! 412: $ KWBOT, A, LDA, B, LDB, JW, KWTOP, QC, LDQC,
! 413: $ JW, KWTOP, ZC, LDZC )
! 414: END DO
! 415:
! 416: K = K-1
! 417: END DO
! 418:
! 419: END IF
! 420:
! 421: * Apply Qc and Zc to rest of the matrix
! 422: IF ( ILSCHUR ) THEN
! 423: ISTARTM = 1
! 424: ISTOPM = N
! 425: ELSE
! 426: ISTARTM = ILO
! 427: ISTOPM = IHI
! 428: END IF
! 429:
! 430: IF ( ISTOPM-IHI > 0 ) THEN
! 431: CALL ZGEMM( 'C', 'N', JW, ISTOPM-IHI, JW, CONE, QC, LDQC,
! 432: $ A( KWTOP, IHI+1 ), LDA, CZERO, WORK, JW )
! 433: CALL ZLACPY( 'ALL', JW, ISTOPM-IHI, WORK, JW, A( KWTOP,
! 434: $ IHI+1 ), LDA )
! 435: CALL ZGEMM( 'C', 'N', JW, ISTOPM-IHI, JW, CONE, QC, LDQC,
! 436: $ B( KWTOP, IHI+1 ), LDB, CZERO, WORK, JW )
! 437: CALL ZLACPY( 'ALL', JW, ISTOPM-IHI, WORK, JW, B( KWTOP,
! 438: $ IHI+1 ), LDB )
! 439: END IF
! 440: IF ( ILQ ) THEN
! 441: CALL ZGEMM( 'N', 'N', N, JW, JW, CONE, Q( 1, KWTOP ), LDQ, QC,
! 442: $ LDQC, CZERO, WORK, N )
! 443: CALL ZLACPY( 'ALL', N, JW, WORK, N, Q( 1, KWTOP ), LDQ )
! 444: END IF
! 445:
! 446: IF ( KWTOP-1-ISTARTM+1 > 0 ) THEN
! 447: CALL ZGEMM( 'N', 'N', KWTOP-ISTARTM, JW, JW, CONE, A( ISTARTM,
! 448: $ KWTOP ), LDA, ZC, LDZC, CZERO, WORK,
! 449: $ KWTOP-ISTARTM )
! 450: CALL ZLACPY( 'ALL', KWTOP-ISTARTM, JW, WORK, KWTOP-ISTARTM,
! 451: $ A( ISTARTM, KWTOP ), LDA )
! 452: CALL ZGEMM( 'N', 'N', KWTOP-ISTARTM, JW, JW, CONE, B( ISTARTM,
! 453: $ KWTOP ), LDB, ZC, LDZC, CZERO, WORK,
! 454: $ KWTOP-ISTARTM )
! 455: CALL ZLACPY( 'ALL', KWTOP-ISTARTM, JW, WORK, KWTOP-ISTARTM,
! 456: $ B( ISTARTM, KWTOP ), LDB )
! 457: END IF
! 458: IF ( ILZ ) THEN
! 459: CALL ZGEMM( 'N', 'N', N, JW, JW, CONE, Z( 1, KWTOP ), LDZ, ZC,
! 460: $ LDZC, CZERO, WORK, N )
! 461: CALL ZLACPY( 'ALL', N, JW, WORK, N, Z( 1, KWTOP ), LDZ )
! 462: END IF
! 463:
! 464: END SUBROUTINE
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaqz2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack