Annotation of rpl/lapack/lapack/dlaqz3.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b DLAQZ3
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DLAQZ3 + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqz3.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqz3.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqz3.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B,
! 22: * $ LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
! 23: * $ ZC, LDZC, WORK, LWORK, 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: * DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
! 32: * $ Q( LDQ, * ), Z( LDZ, * ), ALPHAR( * ), ALPHAI( * ), BETA( * )
! 33: * INTEGER, INTENT( OUT ) :: NS, ND, INFO
! 34: * DOUBLE PRECISION :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
! 35: * ..
! 36: *
! 37: *
! 38: *> \par Purpose:
! 39: * =============
! 40: *>
! 41: *> \verbatim
! 42: *>
! 43: *> DLAQZ3 performs AED
! 44: *> \endverbatim
! 45: *
! 46: * Arguments:
! 47: * ==========
! 48: *
! 49: *> \param[in] ILSCHUR
! 50: *> \verbatim
! 51: *> ILSCHUR is LOGICAL
! 52: *> Determines whether or not to update the full Schur form
! 53: *> \endverbatim
! 54: *>
! 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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] ALPHAR
! 147: *> \verbatim
! 148: *> ALPHAR is DOUBLE PRECISION array, dimension (N)
! 149: *> The real parts of each scalar alpha defining an eigenvalue
! 150: *> of GNEP.
! 151: *> \endverbatim
! 152: *>
! 153: *> \param[out] ALPHAI
! 154: *> \verbatim
! 155: *> ALPHAI is DOUBLE PRECISION array, dimension (N)
! 156: *> The imaginary parts of each scalar alpha defining an
! 157: *> eigenvalue of GNEP.
! 158: *> If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
! 159: *> positive, then the j-th and (j+1)-st eigenvalues are a
! 160: *> complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
! 161: *> \endverbatim
! 162: *>
! 163: *> \param[out] BETA
! 164: *> \verbatim
! 165: *> BETA is DOUBLE PRECISION array, dimension (N)
! 166: *> The scalars beta that define the eigenvalues of GNEP.
! 167: *> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
! 168: *> beta = BETA(j) represent the j-th eigenvalue of the matrix
! 169: *> pair (A,B), in one of the forms lambda = alpha/beta or
! 170: *> mu = beta/alpha. Since either lambda or mu may overflow,
! 171: *> they should not, in general, be computed.
! 172: *> \endverbatim
! 173: *>
! 174: *> \param[in,out] QC
! 175: *> \verbatim
! 176: *> QC is DOUBLE PRECISION array, dimension (LDQC, NW)
! 177: *> \endverbatim
! 178: *>
! 179: *> \param[in] LDQC
! 180: *> \verbatim
! 181: *> LDQC is INTEGER
! 182: *> \endverbatim
! 183: *>
! 184: *> \param[in,out] ZC
! 185: *> \verbatim
! 186: *> ZC is DOUBLE PRECISION array, dimension (LDZC, NW)
! 187: *> \endverbatim
! 188: *>
! 189: *> \param[in] LDZC
! 190: *> \verbatim
! 191: *> LDZ is INTEGER
! 192: *> \endverbatim
! 193: *>
! 194: *> \param[out] WORK
! 195: *> \verbatim
! 196: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 197: *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
! 198: *> \endverbatim
! 199: *>
! 200: *> \param[in] LWORK
! 201: *> \verbatim
! 202: *> LWORK is INTEGER
! 203: *> The dimension of the array WORK. LWORK >= max(1,N).
! 204: *>
! 205: *> If LWORK = -1, then a workspace query is assumed; the routine
! 206: *> only calculates the optimal size of the WORK array, returns
! 207: *> this value as the first entry of the WORK array, and no error
! 208: *> message related to LWORK is issued by XERBLA.
! 209: *> \endverbatim
! 210: *>
! 211: *> \param[in] REC
! 212: *> \verbatim
! 213: *> REC is INTEGER
! 214: *> REC indicates the current recursion level. Should be set
! 215: *> to 0 on first call.
! 216: *> \endverbatim
! 217: *>
! 218: *> \param[out] INFO
! 219: *> \verbatim
! 220: *> INFO is INTEGER
! 221: *> = 0: successful exit
! 222: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 223: *> \endverbatim
! 224: *
! 225: * Authors:
! 226: * ========
! 227: *
! 228: *> \author Thijs Steel, KU Leuven
! 229: *
! 230: *> \date May 2020
! 231: *
! 232: *> \ingroup doubleGEcomputational
! 233: *>
! 234: * =====================================================================
! 235: RECURSIVE SUBROUTINE DLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
! 236: $ A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
! 237: $ ND, ALPHAR, ALPHAI, BETA, QC, LDQC,
! 238: $ ZC, LDZC, WORK, LWORK, REC, INFO )
! 239: IMPLICIT NONE
! 240:
! 241: * Arguments
! 242: LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
! 243: INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
! 244: $ LDQC, LDZC, LWORK, REC
! 245:
! 246: DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
! 247: $ Q( LDQ, * ), Z( LDZ, * ), ALPHAR( * ),
! 248: $ ALPHAI( * ), BETA( * )
! 249: INTEGER, INTENT( OUT ) :: NS, ND, INFO
! 250: DOUBLE PRECISION :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
! 251:
! 252: * Parameters
! 253: DOUBLE PRECISION :: ZERO, ONE, HALF
! 254: PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 )
! 255:
! 256: * Local Scalars
! 257: LOGICAL :: BULGE
! 258: INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, DTGEXC_INFO,
! 259: $ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
! 260: DOUBLE PRECISION :: S, SMLNUM, ULP, SAFMIN, SAFMAX, C1, S1, TEMP
! 261:
! 262: * External Functions
! 263: EXTERNAL :: XERBLA, DTGEXC, DLABAD, DLAQZ0, DLACPY, DLASET,
! 264: $ DLAQZ2, DROT, DLARTG, DLAG2, DGEMM
! 265: DOUBLE PRECISION, EXTERNAL :: DLAMCH
! 266:
! 267: INFO = 0
! 268:
! 269: * Set up deflation window
! 270: JW = MIN( NW, IHI-ILO+1 )
! 271: KWTOP = IHI-JW+1
! 272: IF ( KWTOP .EQ. ILO ) THEN
! 273: S = ZERO
! 274: ELSE
! 275: S = A( KWTOP, KWTOP-1 )
! 276: END IF
! 277:
! 278: * Determine required workspace
! 279: IFST = 1
! 280: ILST = JW
! 281: CALL DTGEXC( .TRUE., .TRUE., JW, A, LDA, B, LDB, QC, LDQC, ZC,
! 282: $ LDZC, IFST, ILST, WORK, -1, DTGEXC_INFO )
! 283: LWORKREQ = INT( WORK( 1 ) )
! 284: CALL DLAQZ0( 'S', 'V', 'V', JW, 1, JW, A( KWTOP, KWTOP ), LDA,
! 285: $ B( KWTOP, KWTOP ), LDB, ALPHAR, ALPHAI, BETA, QC,
! 286: $ LDQC, ZC, LDZC, WORK, -1, REC+1, QZ_SMALL_INFO )
! 287: LWORKREQ = MAX( LWORKREQ, INT( WORK( 1 ) )+2*JW**2 )
! 288: LWORKREQ = MAX( LWORKREQ, N*NW, 2*NW**2+N )
! 289: IF ( LWORK .EQ.-1 ) THEN
! 290: * workspace query, quick return
! 291: WORK( 1 ) = LWORKREQ
! 292: RETURN
! 293: ELSE IF ( LWORK .LT. LWORKREQ ) THEN
! 294: INFO = -26
! 295: END IF
! 296:
! 297: IF( INFO.NE.0 ) THEN
! 298: CALL XERBLA( 'DLAQZ3', -INFO )
! 299: RETURN
! 300: END IF
! 301:
! 302: * Get machine constants
! 303: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
! 304: SAFMAX = ONE/SAFMIN
! 305: CALL DLABAD( SAFMIN, SAFMAX )
! 306: ULP = DLAMCH( 'PRECISION' )
! 307: SMLNUM = SAFMIN*( DBLE( N )/ULP )
! 308:
! 309: IF ( IHI .EQ. KWTOP ) THEN
! 310: * 1 by 1 deflation window, just try a regular deflation
! 311: ALPHAR( KWTOP ) = A( KWTOP, KWTOP )
! 312: ALPHAI( KWTOP ) = ZERO
! 313: BETA( KWTOP ) = B( KWTOP, KWTOP )
! 314: NS = 1
! 315: ND = 0
! 316: IF ( ABS( S ) .LE. MAX( SMLNUM, ULP*ABS( A( KWTOP,
! 317: $ KWTOP ) ) ) ) THEN
! 318: NS = 0
! 319: ND = 1
! 320: IF ( KWTOP .GT. ILO ) THEN
! 321: A( KWTOP, KWTOP-1 ) = ZERO
! 322: END IF
! 323: END IF
! 324: END IF
! 325:
! 326:
! 327: * Store window in case of convergence failure
! 328: CALL DLACPY( 'ALL', JW, JW, A( KWTOP, KWTOP ), LDA, WORK, JW )
! 329: CALL DLACPY( 'ALL', JW, JW, B( KWTOP, KWTOP ), LDB, WORK( JW**2+
! 330: $ 1 ), JW )
! 331:
! 332: * Transform window to real schur form
! 333: CALL DLASET( 'FULL', JW, JW, ZERO, ONE, QC, LDQC )
! 334: CALL DLASET( 'FULL', JW, JW, ZERO, ONE, ZC, LDZC )
! 335: CALL DLAQZ0( 'S', 'V', 'V', JW, 1, JW, A( KWTOP, KWTOP ), LDA,
! 336: $ B( KWTOP, KWTOP ), LDB, ALPHAR, ALPHAI, BETA, QC,
! 337: $ LDQC, ZC, LDZC, WORK( 2*JW**2+1 ), LWORK-2*JW**2,
! 338: $ REC+1, QZ_SMALL_INFO )
! 339:
! 340: IF( QZ_SMALL_INFO .NE. 0 ) THEN
! 341: * Convergence failure, restore the window and exit
! 342: ND = 0
! 343: NS = JW-QZ_SMALL_INFO
! 344: CALL DLACPY( 'ALL', JW, JW, WORK, JW, A( KWTOP, KWTOP ), LDA )
! 345: CALL DLACPY( 'ALL', JW, JW, WORK( JW**2+1 ), JW, B( KWTOP,
! 346: $ KWTOP ), LDB )
! 347: RETURN
! 348: END IF
! 349:
! 350: * Deflation detection loop
! 351: IF ( KWTOP .EQ. ILO .OR. S .EQ. ZERO ) THEN
! 352: KWBOT = KWTOP-1
! 353: ELSE
! 354: KWBOT = IHI
! 355: K = 1
! 356: K2 = 1
! 357: DO WHILE ( K .LE. JW )
! 358: BULGE = .FALSE.
! 359: IF ( KWBOT-KWTOP+1 .GE. 2 ) THEN
! 360: BULGE = A( KWBOT, KWBOT-1 ) .NE. ZERO
! 361: END IF
! 362: IF ( BULGE ) THEN
! 363:
! 364: * Try to deflate complex conjugate eigenvalue pair
! 365: TEMP = ABS( A( KWBOT, KWBOT ) )+SQRT( ABS( A( KWBOT,
! 366: $ KWBOT-1 ) ) )*SQRT( ABS( A( KWBOT-1, KWBOT ) ) )
! 367: IF( TEMP .EQ. ZERO )THEN
! 368: TEMP = ABS( S )
! 369: END IF
! 370: IF ( MAX( ABS( S*QC( 1, KWBOT-KWTOP ) ), ABS( S*QC( 1,
! 371: $ KWBOT-KWTOP+1 ) ) ) .LE. MAX( SMLNUM,
! 372: $ ULP*TEMP ) ) THEN
! 373: * Deflatable
! 374: KWBOT = KWBOT-2
! 375: ELSE
! 376: * Not deflatable, move out of the way
! 377: IFST = KWBOT-KWTOP+1
! 378: ILST = K2
! 379: CALL DTGEXC( .TRUE., .TRUE., JW, A( KWTOP, KWTOP ),
! 380: $ LDA, B( KWTOP, KWTOP ), LDB, QC, LDQC,
! 381: $ ZC, LDZC, IFST, ILST, WORK, LWORK,
! 382: $ DTGEXC_INFO )
! 383: K2 = K2+2
! 384: END IF
! 385: K = K+2
! 386: ELSE
! 387:
! 388: * Try to deflate real eigenvalue
! 389: TEMP = ABS( A( KWBOT, KWBOT ) )
! 390: IF( TEMP .EQ. ZERO ) THEN
! 391: TEMP = ABS( S )
! 392: END IF
! 393: IF ( ( ABS( S*QC( 1, KWBOT-KWTOP+1 ) ) ) .LE. MAX( ULP*
! 394: $ TEMP, SMLNUM ) ) THEN
! 395: * Deflatable
! 396: KWBOT = KWBOT-1
! 397: ELSE
! 398: * Not deflatable, move out of the way
! 399: IFST = KWBOT-KWTOP+1
! 400: ILST = K2
! 401: CALL DTGEXC( .TRUE., .TRUE., JW, A( KWTOP, KWTOP ),
! 402: $ LDA, B( KWTOP, KWTOP ), LDB, QC, LDQC,
! 403: $ ZC, LDZC, IFST, ILST, WORK, LWORK,
! 404: $ DTGEXC_INFO )
! 405: K2 = K2+1
! 406: END IF
! 407:
! 408: K = K+1
! 409:
! 410: END IF
! 411: END DO
! 412: END IF
! 413:
! 414: * Store eigenvalues
! 415: ND = IHI-KWBOT
! 416: NS = JW-ND
! 417: K = KWTOP
! 418: DO WHILE ( K .LE. IHI )
! 419: BULGE = .FALSE.
! 420: IF ( K .LT. IHI ) THEN
! 421: IF ( A( K+1, K ) .NE. ZERO ) THEN
! 422: BULGE = .TRUE.
! 423: END IF
! 424: END IF
! 425: IF ( BULGE ) THEN
! 426: * 2x2 eigenvalue block
! 427: CALL DLAG2( A( K, K ), LDA, B( K, K ), LDB, SAFMIN,
! 428: $ BETA( K ), BETA( K+1 ), ALPHAR( K ),
! 429: $ ALPHAR( K+1 ), ALPHAI( K ) )
! 430: ALPHAI( K+1 ) = -ALPHAI( K )
! 431: K = K+2
! 432: ELSE
! 433: * 1x1 eigenvalue block
! 434: ALPHAR( K ) = A( K, K )
! 435: ALPHAI( K ) = ZERO
! 436: BETA( K ) = B( K, K )
! 437: K = K+1
! 438: END IF
! 439: END DO
! 440:
! 441: IF ( KWTOP .NE. ILO .AND. S .NE. ZERO ) THEN
! 442: * Reflect spike back, this will create optimally packed bulges
! 443: A( KWTOP:KWBOT, KWTOP-1 ) = A( KWTOP, KWTOP-1 )*QC( 1,
! 444: $ 1:JW-ND )
! 445: DO K = KWBOT-1, KWTOP, -1
! 446: CALL DLARTG( A( K, KWTOP-1 ), A( K+1, KWTOP-1 ), C1, S1,
! 447: $ TEMP )
! 448: A( K, KWTOP-1 ) = TEMP
! 449: A( K+1, KWTOP-1 ) = ZERO
! 450: K2 = MAX( KWTOP, K-1 )
! 451: CALL DROT( IHI-K2+1, A( K, K2 ), LDA, A( K+1, K2 ), LDA, C1,
! 452: $ S1 )
! 453: CALL DROT( IHI-( K-1 )+1, B( K, K-1 ), LDB, B( K+1, K-1 ),
! 454: $ LDB, C1, S1 )
! 455: CALL DROT( JW, QC( 1, K-KWTOP+1 ), 1, QC( 1, K+1-KWTOP+1 ),
! 456: $ 1, C1, S1 )
! 457: END DO
! 458:
! 459: * Chase bulges down
! 460: ISTARTM = KWTOP
! 461: ISTOPM = IHI
! 462: K = KWBOT-1
! 463: DO WHILE ( K .GE. KWTOP )
! 464: IF ( ( K .GE. KWTOP+1 ) .AND. A( K+1, K-1 ) .NE. ZERO ) THEN
! 465:
! 466: * Move double pole block down and remove it
! 467: DO K2 = K-1, KWBOT-2
! 468: CALL DLAQZ2( .TRUE., .TRUE., K2, KWTOP, KWTOP+JW-1,
! 469: $ KWBOT, A, LDA, B, LDB, JW, KWTOP, QC,
! 470: $ LDQC, JW, KWTOP, ZC, LDZC )
! 471: END DO
! 472:
! 473: K = K-2
! 474: ELSE
! 475:
! 476: * k points to single shift
! 477: DO K2 = K, KWBOT-2
! 478:
! 479: * Move shift down
! 480: CALL DLARTG( B( K2+1, K2+1 ), B( K2+1, K2 ), C1, S1,
! 481: $ TEMP )
! 482: B( K2+1, K2+1 ) = TEMP
! 483: B( K2+1, K2 ) = ZERO
! 484: CALL DROT( K2+2-ISTARTM+1, A( ISTARTM, K2+1 ), 1,
! 485: $ A( ISTARTM, K2 ), 1, C1, S1 )
! 486: CALL DROT( K2-ISTARTM+1, B( ISTARTM, K2+1 ), 1,
! 487: $ B( ISTARTM, K2 ), 1, C1, S1 )
! 488: CALL DROT( JW, ZC( 1, K2+1-KWTOP+1 ), 1, ZC( 1,
! 489: $ K2-KWTOP+1 ), 1, C1, S1 )
! 490:
! 491: CALL DLARTG( A( K2+1, K2 ), A( K2+2, K2 ), C1, S1,
! 492: $ TEMP )
! 493: A( K2+1, K2 ) = TEMP
! 494: A( K2+2, K2 ) = ZERO
! 495: CALL DROT( ISTOPM-K2, A( K2+1, K2+1 ), LDA, A( K2+2,
! 496: $ K2+1 ), LDA, C1, S1 )
! 497: CALL DROT( ISTOPM-K2, B( K2+1, K2+1 ), LDB, B( K2+2,
! 498: $ K2+1 ), LDB, C1, S1 )
! 499: CALL DROT( JW, QC( 1, K2+1-KWTOP+1 ), 1, QC( 1,
! 500: $ K2+2-KWTOP+1 ), 1, C1, S1 )
! 501:
! 502: END DO
! 503:
! 504: * Remove the shift
! 505: CALL DLARTG( B( KWBOT, KWBOT ), B( KWBOT, KWBOT-1 ), C1,
! 506: $ S1, TEMP )
! 507: B( KWBOT, KWBOT ) = TEMP
! 508: B( KWBOT, KWBOT-1 ) = ZERO
! 509: CALL DROT( KWBOT-ISTARTM, B( ISTARTM, KWBOT ), 1,
! 510: $ B( ISTARTM, KWBOT-1 ), 1, C1, S1 )
! 511: CALL DROT( KWBOT-ISTARTM+1, A( ISTARTM, KWBOT ), 1,
! 512: $ A( ISTARTM, KWBOT-1 ), 1, C1, S1 )
! 513: CALL DROT( JW, ZC( 1, KWBOT-KWTOP+1 ), 1, ZC( 1,
! 514: $ KWBOT-1-KWTOP+1 ), 1, C1, S1 )
! 515:
! 516: K = K-1
! 517: END IF
! 518: END DO
! 519:
! 520: END IF
! 521:
! 522: * Apply Qc and Zc to rest of the matrix
! 523: IF ( ILSCHUR ) THEN
! 524: ISTARTM = 1
! 525: ISTOPM = N
! 526: ELSE
! 527: ISTARTM = ILO
! 528: ISTOPM = IHI
! 529: END IF
! 530:
! 531: IF ( ISTOPM-IHI > 0 ) THEN
! 532: CALL DGEMM( 'T', 'N', JW, ISTOPM-IHI, JW, ONE, QC, LDQC,
! 533: $ A( KWTOP, IHI+1 ), LDA, ZERO, WORK, JW )
! 534: CALL DLACPY( 'ALL', JW, ISTOPM-IHI, WORK, JW, A( KWTOP,
! 535: $ IHI+1 ), LDA )
! 536: CALL DGEMM( 'T', 'N', JW, ISTOPM-IHI, JW, ONE, QC, LDQC,
! 537: $ B( KWTOP, IHI+1 ), LDB, ZERO, WORK, JW )
! 538: CALL DLACPY( 'ALL', JW, ISTOPM-IHI, WORK, JW, B( KWTOP,
! 539: $ IHI+1 ), LDB )
! 540: END IF
! 541: IF ( ILQ ) THEN
! 542: CALL DGEMM( 'N', 'N', N, JW, JW, ONE, Q( 1, KWTOP ), LDQ, QC,
! 543: $ LDQC, ZERO, WORK, N )
! 544: CALL DLACPY( 'ALL', N, JW, WORK, N, Q( 1, KWTOP ), LDQ )
! 545: END IF
! 546:
! 547: IF ( KWTOP-1-ISTARTM+1 > 0 ) THEN
! 548: CALL DGEMM( 'N', 'N', KWTOP-ISTARTM, JW, JW, ONE, A( ISTARTM,
! 549: $ KWTOP ), LDA, ZC, LDZC, ZERO, WORK,
! 550: $ KWTOP-ISTARTM )
! 551: CALL DLACPY( 'ALL', KWTOP-ISTARTM, JW, WORK, KWTOP-ISTARTM,
! 552: $ A( ISTARTM, KWTOP ), LDA )
! 553: CALL DGEMM( 'N', 'N', KWTOP-ISTARTM, JW, JW, ONE, B( ISTARTM,
! 554: $ KWTOP ), LDB, ZC, LDZC, ZERO, WORK,
! 555: $ KWTOP-ISTARTM )
! 556: CALL DLACPY( 'ALL', KWTOP-ISTARTM, JW, WORK, KWTOP-ISTARTM,
! 557: $ B( ISTARTM, KWTOP ), LDB )
! 558: END IF
! 559: IF ( ILZ ) THEN
! 560: CALL DGEMM( 'N', 'N', N, JW, JW, ONE, Z( 1, KWTOP ), LDZ, ZC,
! 561: $ LDZC, ZERO, WORK, N )
! 562: CALL DLACPY( 'ALL', N, JW, WORK, N, Z( 1, KWTOP ), LDZ )
! 563: END IF
! 564:
! 565: END SUBROUTINE
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dlaqz3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack