Annotation of rpl/lapack/lapack/dlaqz4.f, revision 1.1

1.1     ! bertrand    1: *> \brief \b DLAQZ4
        !             2: *
        !             3: *  =========== DOCUMENTATION ===========
        !             4: *
        !             5: * Online html documentation available at
        !             6: *            http://www.netlib.org/lapack/explore-html/
        !             7: *
        !             8: *> \htmlonly
        !             9: *> Download DLAQZ4 + dependencies
        !            10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqz4.f">
        !            11: *> [TGZ]</a>
        !            12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqz4.f">
        !            13: *> [ZIP]</a>
        !            14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqz4.f">
        !            15: *> [TXT]</a>
        !            16: *> \endhtmlonly
        !            17: *
        !            18: *  Definition:
        !            19: *  ===========
        !            20: *
        !            21: *      SUBROUTINE DLAQZ4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
        !            22: *     $    NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ,
        !            23: *     $    QC, LDQC, ZC, LDZC, WORK, LWORK, INFO )
        !            24: *      IMPLICIT NONE
        !            25: *
        !            26: *      Function arguments
        !            27: *      LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
        !            28: *      INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
        !            29: *     $    NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
        !            30: *
        !            31: *      DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
        !            32: *     $    Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ),
        !            33: *     $    WORK( * ), SR( * ), SI( * ), SS( * )
        !            34: *
        !            35: *      INTEGER, INTENT( OUT ) :: INFO
        !            36: *       ..
        !            37: *
        !            38: *
        !            39: *> \par Purpose:
        !            40: *  =============
        !            41: *>
        !            42: *> \verbatim
        !            43: *>
        !            44: *> DLAQZ4 Executes a single multishift QZ sweep
        !            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: *> \endverbatim
        !            82: *>
        !            83: *> \param[in] NSHIFTS
        !            84: *> \verbatim
        !            85: *>          NSHIFTS is INTEGER
        !            86: *>          The desired number of shifts to use
        !            87: *> \endverbatim
        !            88: *>
        !            89: *> \param[in] NBLOCK_DESIRED
        !            90: *> \verbatim
        !            91: *>          NBLOCK_DESIRED is INTEGER
        !            92: *>          The desired size of the computational windows
        !            93: *> \endverbatim
        !            94: *>
        !            95: *> \param[in] SR
        !            96: *> \verbatim
        !            97: *>          SR is DOUBLE PRECISION array. SR contains
        !            98: *>          the real parts of the shifts to use.
        !            99: *> \endverbatim
        !           100: *>
        !           101: *> \param[in] SI
        !           102: *> \verbatim
        !           103: *>          SI is DOUBLE PRECISION array. SI contains
        !           104: *>          the imaginary parts of the shifts to use.
        !           105: *> \endverbatim
        !           106: *>
        !           107: *> \param[in] SS
        !           108: *> \verbatim
        !           109: *>          SS is DOUBLE PRECISION array. SS contains
        !           110: *>          the scale of the shifts to use.
        !           111: *> \endverbatim
        !           112: *>
        !           113: *> \param[in,out] A
        !           114: *> \verbatim
        !           115: *>          A is DOUBLE PRECISION array, dimension (LDA, N)
        !           116: *> \endverbatim
        !           117: *>
        !           118: *> \param[in] LDA
        !           119: *> \verbatim
        !           120: *>          LDA is INTEGER
        !           121: *>          The leading dimension of the array A.  LDA >= max( 1, N ).
        !           122: *> \endverbatim
        !           123: *>
        !           124: *> \param[in,out] B
        !           125: *> \verbatim
        !           126: *>          B is DOUBLE PRECISION array, dimension (LDB, N)
        !           127: *> \endverbatim
        !           128: *>
        !           129: *> \param[in] LDB
        !           130: *> \verbatim
        !           131: *>          LDB is INTEGER
        !           132: *>          The leading dimension of the array B.  LDB >= max( 1, N ).
        !           133: *> \endverbatim
        !           134: *>
        !           135: *> \param[in,out] Q
        !           136: *> \verbatim
        !           137: *>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
        !           138: *> \endverbatim
        !           139: *>
        !           140: *> \param[in] LDQ
        !           141: *> \verbatim
        !           142: *>          LDQ is INTEGER
        !           143: *> \endverbatim
        !           144: *>
        !           145: *> \param[in,out] Z
        !           146: *> \verbatim
        !           147: *>          Z is DOUBLE PRECISION array, dimension (LDZ, N)
        !           148: *> \endverbatim
        !           149: *>
        !           150: *> \param[in] LDZ
        !           151: *> \verbatim
        !           152: *>          LDZ is INTEGER
        !           153: *> \endverbatim
        !           154: *>
        !           155: *> \param[in,out] QC
        !           156: *> \verbatim
        !           157: *>          QC is DOUBLE PRECISION array, dimension (LDQC, NBLOCK_DESIRED)
        !           158: *> \endverbatim
        !           159: *>
        !           160: *> \param[in] LDQC
        !           161: *> \verbatim
        !           162: *>          LDQC is INTEGER
        !           163: *> \endverbatim
        !           164: *>
        !           165: *> \param[in,out] ZC
        !           166: *> \verbatim
        !           167: *>          ZC is DOUBLE PRECISION array, dimension (LDZC, NBLOCK_DESIRED)
        !           168: *> \endverbatim
        !           169: *>
        !           170: *> \param[in] LDZC
        !           171: *> \verbatim
        !           172: *>          LDZ is INTEGER
        !           173: *> \endverbatim
        !           174: *>
        !           175: *> \param[out] WORK
        !           176: *> \verbatim
        !           177: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
        !           178: *>          On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
        !           179: *> \endverbatim
        !           180: *>
        !           181: *> \param[in] LWORK
        !           182: *> \verbatim
        !           183: *>          LWORK is INTEGER
        !           184: *>          The dimension of the array WORK.  LWORK >= max(1,N).
        !           185: *>
        !           186: *>          If LWORK = -1, then a workspace query is assumed; the routine
        !           187: *>          only calculates the optimal size of the WORK array, returns
        !           188: *>          this value as the first entry of the WORK array, and no error
        !           189: *>          message related to LWORK is issued by XERBLA.
        !           190: *> \endverbatim
        !           191: *>
        !           192: *> \param[out] INFO
        !           193: *> \verbatim
        !           194: *>          INFO is INTEGER
        !           195: *>          = 0: successful exit
        !           196: *>          < 0: if INFO = -i, the i-th argument had an illegal value
        !           197: *> \endverbatim
        !           198: *
        !           199: *  Authors:
        !           200: *  ========
        !           201: *
        !           202: *> \author Thijs Steel, KU Leuven
        !           203: *
        !           204: *> \date May 2020
        !           205: *
        !           206: *> \ingroup doubleGEcomputational
        !           207: *>
        !           208: *  =====================================================================
        !           209:       SUBROUTINE DLAQZ4( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
        !           210:      $                   NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q,
        !           211:      $                   LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK,
        !           212:      $                   INFO )
        !           213:       IMPLICIT NONE
        !           214: 
        !           215: *     Function arguments
        !           216:       LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
        !           217:       INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
        !           218:      $         NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
        !           219: 
        !           220:       DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
        !           221:      $                  Q( LDQ, * ), Z( LDZ, * ), QC( LDQC, * ),
        !           222:      $                  ZC( LDZC, * ), WORK( * ), SR( * ), SI( * ),
        !           223:      $                  SS( * )
        !           224: 
        !           225:       INTEGER, INTENT( OUT ) :: INFO
        !           226: 
        !           227: *     Parameters
        !           228:       DOUBLE PRECISION :: ZERO, ONE, HALF
        !           229:       PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 )
        !           230: 
        !           231: *     Local scalars
        !           232:       INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
        !           233:      $           ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
        !           234:       DOUBLE PRECISION :: TEMP, V( 3 ), C1, S1, C2, S2, SWAP
        !           235: *
        !           236: *     External functions
        !           237:       EXTERNAL :: XERBLA, DGEMM, DLAQZ1, DLAQZ2, DLASET, DLARTG, DROT,
        !           238:      $            DLACPY
        !           239: 
        !           240:       INFO = 0
        !           241:       IF ( NBLOCK_DESIRED .LT. NSHIFTS+1 ) THEN
        !           242:          INFO = -8
        !           243:       END IF
        !           244:       IF ( LWORK .EQ.-1 ) THEN
        !           245: *        workspace query, quick return
        !           246:          WORK( 1 ) = N*NBLOCK_DESIRED
        !           247:          RETURN
        !           248:       ELSE IF ( LWORK .LT. N*NBLOCK_DESIRED ) THEN
        !           249:          INFO = -25
        !           250:       END IF
        !           251: 
        !           252:       IF( INFO.NE.0 ) THEN
        !           253:          CALL XERBLA( 'DLAQZ4', -INFO )
        !           254:          RETURN
        !           255:       END IF
        !           256: 
        !           257: *     Executable statements
        !           258: 
        !           259:       IF ( NSHIFTS .LT. 2 ) THEN
        !           260:          RETURN
        !           261:       END IF
        !           262: 
        !           263:       IF ( ILO .GE. IHI ) THEN
        !           264:          RETURN
        !           265:       END IF
        !           266: 
        !           267:       IF ( ILSCHUR ) THEN
        !           268:          ISTARTM = 1
        !           269:          ISTOPM = N
        !           270:       ELSE
        !           271:          ISTARTM = ILO
        !           272:          ISTOPM = IHI
        !           273:       END IF
        !           274: 
        !           275: *     Shuffle shifts into pairs of real shifts and pairs
        !           276: *     of complex conjugate shifts assuming complex
        !           277: *     conjugate shifts are already adjacent to one
        !           278: *     another
        !           279: 
        !           280:       DO I = 1, NSHIFTS-2, 2
        !           281:          IF( SI( I ).NE.-SI( I+1 ) ) THEN
        !           282: *
        !           283:             SWAP = SR( I )
        !           284:             SR( I ) = SR( I+1 )
        !           285:             SR( I+1 ) = SR( I+2 )
        !           286:             SR( I+2 ) = SWAP
        !           287: 
        !           288:             SWAP = SI( I )
        !           289:             SI( I ) = SI( I+1 )
        !           290:             SI( I+1 ) = SI( I+2 )
        !           291:             SI( I+2 ) = SWAP
        !           292:             
        !           293:             SWAP = SS( I )
        !           294:             SS( I ) = SS( I+1 )
        !           295:             SS( I+1 ) = SS( I+2 )
        !           296:             SS( I+2 ) = SWAP
        !           297:          END IF
        !           298:       END DO
        !           299: 
        !           300: *     NSHFTS is supposed to be even, but if it is odd,
        !           301: *     then simply reduce it by one.  The shuffle above
        !           302: *     ensures that the dropped shift is real and that
        !           303: *     the remaining shifts are paired.
        !           304: 
        !           305:       NS = NSHIFTS-MOD( NSHIFTS, 2 )
        !           306:       NPOS = MAX( NBLOCK_DESIRED-NS, 1 )
        !           307: 
        !           308: *     The following block introduces the shifts and chases
        !           309: *     them down one by one just enough to make space for
        !           310: *     the other shifts. The near-the-diagonal block is
        !           311: *     of size (ns+1) x ns.
        !           312: 
        !           313:       CALL DLASET( 'FULL', NS+1, NS+1, ZERO, ONE, QC, LDQC )
        !           314:       CALL DLASET( 'FULL', NS, NS, ZERO, ONE, ZC, LDZC )
        !           315: 
        !           316:       DO I = 1, NS, 2
        !           317: *        Introduce the shift
        !           318:          CALL DLAQZ1( A( ILO, ILO ), LDA, B( ILO, ILO ), LDB, SR( I ),
        !           319:      $                SR( I+1 ), SI( I ), SS( I ), SS( I+1 ), V )
        !           320: 
        !           321:          TEMP = V( 2 )
        !           322:          CALL DLARTG( TEMP, V( 3 ), C1, S1, V( 2 ) )
        !           323:          CALL DLARTG( V( 1 ), V( 2 ), C2, S2, TEMP )
        !           324: 
        !           325:          CALL DROT( NS, A( ILO+1, ILO ), LDA, A( ILO+2, ILO ), LDA, C1,
        !           326:      $              S1 )
        !           327:          CALL DROT( NS, A( ILO, ILO ), LDA, A( ILO+1, ILO ), LDA, C2,
        !           328:      $              S2 )
        !           329:          CALL DROT( NS, B( ILO+1, ILO ), LDB, B( ILO+2, ILO ), LDB, C1,
        !           330:      $              S1 )
        !           331:          CALL DROT( NS, B( ILO, ILO ), LDB, B( ILO+1, ILO ), LDB, C2,
        !           332:      $              S2 )
        !           333:          CALL DROT( NS+1, QC( 1, 2 ), 1, QC( 1, 3 ), 1, C1, S1 )
        !           334:          CALL DROT( NS+1, QC( 1, 1 ), 1, QC( 1, 2 ), 1, C2, S2 )
        !           335: 
        !           336: *        Chase the shift down
        !           337:          DO J = 1, NS-1-I
        !           338: 
        !           339:             CALL DLAQZ2( .TRUE., .TRUE., J, 1, NS, IHI-ILO+1, A( ILO,
        !           340:      $                   ILO ), LDA, B( ILO, ILO ), LDB, NS+1, 1, QC,
        !           341:      $                   LDQC, NS, 1, ZC, LDZC )
        !           342: 
        !           343:          END DO
        !           344: 
        !           345:       END DO
        !           346: 
        !           347: *     Update the rest of the pencil
        !           348: 
        !           349: *     Update A(ilo:ilo+ns,ilo+ns:istopm) and B(ilo:ilo+ns,ilo+ns:istopm)
        !           350: *     from the left with Qc(1:ns+1,1:ns+1)'
        !           351:       SHEIGHT = NS+1
        !           352:       SWIDTH = ISTOPM-( ILO+NS )+1
        !           353:       IF ( SWIDTH > 0 ) THEN
        !           354:          CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
        !           355:      $               A( ILO, ILO+NS ), LDA, ZERO, WORK, SHEIGHT )
        !           356:          CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ILO,
        !           357:      $                ILO+NS ), LDA )
        !           358:          CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
        !           359:      $               B( ILO, ILO+NS ), LDB, ZERO, WORK, SHEIGHT )
        !           360:          CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ILO,
        !           361:      $                ILO+NS ), LDB )
        !           362:       END IF
        !           363:       IF ( ILQ ) THEN
        !           364:          CALL DGEMM( 'N', 'N', N, SHEIGHT, SHEIGHT, ONE, Q( 1, ILO ),
        !           365:      $               LDQ, QC, LDQC, ZERO, WORK, N )
        !           366:          CALL DLACPY( 'ALL', N, SHEIGHT, WORK, N, Q( 1, ILO ), LDQ )
        !           367:       END IF
        !           368: 
        !           369: *     Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
        !           370: *     from the right with Zc(1:ns,1:ns)
        !           371:       SHEIGHT = ILO-1-ISTARTM+1
        !           372:       SWIDTH = NS
        !           373:       IF ( SHEIGHT > 0 ) THEN
        !           374:          CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, A( ISTARTM,
        !           375:      $               ILO ), LDA, ZC, LDZC, ZERO, WORK, SHEIGHT )
        !           376:          CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ISTARTM,
        !           377:      $                ILO ), LDA )
        !           378:          CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, B( ISTARTM,
        !           379:      $               ILO ), LDB, ZC, LDZC, ZERO, WORK, SHEIGHT )
        !           380:          CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ISTARTM,
        !           381:      $                ILO ), LDB )
        !           382:       END IF
        !           383:       IF ( ILZ ) THEN
        !           384:          CALL DGEMM( 'N', 'N', N, SWIDTH, SWIDTH, ONE, Z( 1, ILO ), LDZ,
        !           385:      $               ZC, LDZC, ZERO, WORK, N )
        !           386:          CALL DLACPY( 'ALL', N, SWIDTH, WORK, N, Z( 1, ILO ), LDZ )
        !           387:       END IF
        !           388: 
        !           389: *     The following block chases the shifts down to the bottom
        !           390: *     right block. If possible, a shift is moved down npos
        !           391: *     positions at a time
        !           392: 
        !           393:       K = ILO
        !           394:       DO WHILE ( K < IHI-NS )
        !           395:          NP = MIN( IHI-NS-K, NPOS )
        !           396: *        Size of the near-the-diagonal block
        !           397:          NBLOCK = NS+NP
        !           398: *        istartb points to the first row we will be updating
        !           399:          ISTARTB = K+1
        !           400: *        istopb points to the last column we will be updating
        !           401:          ISTOPB = K+NBLOCK-1
        !           402: 
        !           403:          CALL DLASET( 'FULL', NS+NP, NS+NP, ZERO, ONE, QC, LDQC )
        !           404:          CALL DLASET( 'FULL', NS+NP, NS+NP, ZERO, ONE, ZC, LDZC )
        !           405: 
        !           406: *        Near the diagonal shift chase
        !           407:          DO I = NS-1, 0, -2
        !           408:             DO J = 0, NP-1
        !           409: *              Move down the block with index k+i+j-1, updating
        !           410: *              the (ns+np x ns+np) block:
        !           411: *              (k:k+ns+np,k:k+ns+np-1)
        !           412:                CALL DLAQZ2( .TRUE., .TRUE., K+I+J-1, ISTARTB, ISTOPB,
        !           413:      $                      IHI, A, LDA, B, LDB, NBLOCK, K+1, QC, LDQC,
        !           414:      $                      NBLOCK, K, ZC, LDZC )
        !           415:             END DO
        !           416:          END DO
        !           417: 
        !           418: *        Update rest of the pencil
        !           419: 
        !           420: *        Update A(k+1:k+ns+np, k+ns+np:istopm) and
        !           421: *        B(k+1:k+ns+np, k+ns+np:istopm)
        !           422: *        from the left with Qc(1:ns+np,1:ns+np)'
        !           423:          SHEIGHT = NS+NP
        !           424:          SWIDTH = ISTOPM-( K+NS+NP )+1
        !           425:          IF ( SWIDTH > 0 ) THEN
        !           426:             CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC,
        !           427:      $                  LDQC, A( K+1, K+NS+NP ), LDA, ZERO, WORK,
        !           428:      $                  SHEIGHT )
        !           429:             CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( K+1,
        !           430:      $                   K+NS+NP ), LDA )
        !           431:             CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC,
        !           432:      $                  LDQC, B( K+1, K+NS+NP ), LDB, ZERO, WORK,
        !           433:      $                  SHEIGHT )
        !           434:             CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( K+1,
        !           435:      $                   K+NS+NP ), LDB )
        !           436:          END IF
        !           437:          IF ( ILQ ) THEN
        !           438:             CALL DGEMM( 'N', 'N', N, NBLOCK, NBLOCK, ONE, Q( 1, K+1 ),
        !           439:      $                  LDQ, QC, LDQC, ZERO, WORK, N )
        !           440:             CALL DLACPY( 'ALL', N, NBLOCK, WORK, N, Q( 1, K+1 ), LDQ )
        !           441:          END IF
        !           442: 
        !           443: *        Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
        !           444: *        from the right with Zc(1:ns+np,1:ns+np)
        !           445:          SHEIGHT = K-ISTARTM+1
        !           446:          SWIDTH = NBLOCK
        !           447:          IF ( SHEIGHT > 0 ) THEN
        !           448:             CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE,
        !           449:      $                  A( ISTARTM, K ), LDA, ZC, LDZC, ZERO, WORK,
        !           450:      $                  SHEIGHT )
        !           451:             CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
        !           452:      $                   A( ISTARTM, K ), LDA )
        !           453:             CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE,
        !           454:      $                  B( ISTARTM, K ), LDB, ZC, LDZC, ZERO, WORK,
        !           455:      $                  SHEIGHT )
        !           456:             CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
        !           457:      $                   B( ISTARTM, K ), LDB )
        !           458:          END IF
        !           459:          IF ( ILZ ) THEN
        !           460:             CALL DGEMM( 'N', 'N', N, NBLOCK, NBLOCK, ONE, Z( 1, K ),
        !           461:      $                  LDZ, ZC, LDZC, ZERO, WORK, N )
        !           462:             CALL DLACPY( 'ALL', N, NBLOCK, WORK, N, Z( 1, K ), LDZ )
        !           463:          END IF
        !           464: 
        !           465:          K = K+NP
        !           466: 
        !           467:       END DO
        !           468: 
        !           469: *     The following block removes the shifts from the bottom right corner
        !           470: *     one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).
        !           471: 
        !           472:       CALL DLASET( 'FULL', NS, NS, ZERO, ONE, QC, LDQC )
        !           473:       CALL DLASET( 'FULL', NS+1, NS+1, ZERO, ONE, ZC, LDZC )
        !           474: 
        !           475: *     istartb points to the first row we will be updating
        !           476:       ISTARTB = IHI-NS+1
        !           477: *     istopb points to the last column we will be updating
        !           478:       ISTOPB = IHI
        !           479: 
        !           480:       DO I = 1, NS, 2
        !           481: *        Chase the shift down to the bottom right corner
        !           482:          DO ISHIFT = IHI-I-1, IHI-2
        !           483:             CALL DLAQZ2( .TRUE., .TRUE., ISHIFT, ISTARTB, ISTOPB, IHI,
        !           484:      $                   A, LDA, B, LDB, NS, IHI-NS+1, QC, LDQC, NS+1,
        !           485:      $                   IHI-NS, ZC, LDZC )
        !           486:          END DO
        !           487:          
        !           488:       END DO
        !           489: 
        !           490: *     Update rest of the pencil
        !           491: 
        !           492: *     Update A(ihi-ns+1:ihi, ihi+1:istopm)
        !           493: *     from the left with Qc(1:ns,1:ns)'
        !           494:       SHEIGHT = NS
        !           495:       SWIDTH = ISTOPM-( IHI+1 )+1
        !           496:       IF ( SWIDTH > 0 ) THEN
        !           497:          CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
        !           498:      $               A( IHI-NS+1, IHI+1 ), LDA, ZERO, WORK, SHEIGHT )
        !           499:          CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
        !           500:      $                A( IHI-NS+1, IHI+1 ), LDA )
        !           501:          CALL DGEMM( 'T', 'N', SHEIGHT, SWIDTH, SHEIGHT, ONE, QC, LDQC,
        !           502:      $               B( IHI-NS+1, IHI+1 ), LDB, ZERO, WORK, SHEIGHT )
        !           503:          CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT,
        !           504:      $                B( IHI-NS+1, IHI+1 ), LDB )
        !           505:       END IF
        !           506:       IF ( ILQ ) THEN
        !           507:          CALL DGEMM( 'N', 'N', N, NS, NS, ONE, Q( 1, IHI-NS+1 ), LDQ,
        !           508:      $               QC, LDQC, ZERO, WORK, N )
        !           509:          CALL DLACPY( 'ALL', N, NS, WORK, N, Q( 1, IHI-NS+1 ), LDQ )
        !           510:       END IF
        !           511: 
        !           512: *     Update A(istartm:ihi-ns,ihi-ns:ihi)
        !           513: *     from the right with Zc(1:ns+1,1:ns+1)
        !           514:       SHEIGHT = IHI-NS-ISTARTM+1
        !           515:       SWIDTH = NS+1
        !           516:       IF ( SHEIGHT > 0 ) THEN
        !           517:          CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, A( ISTARTM,
        !           518:      $               IHI-NS ), LDA, ZC, LDZC, ZERO, WORK, SHEIGHT )
        !           519:          CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, A( ISTARTM,
        !           520:      $                IHI-NS ), LDA )
        !           521:          CALL DGEMM( 'N', 'N', SHEIGHT, SWIDTH, SWIDTH, ONE, B( ISTARTM,
        !           522:      $               IHI-NS ), LDB, ZC, LDZC, ZERO, WORK, SHEIGHT )
        !           523:          CALL DLACPY( 'ALL', SHEIGHT, SWIDTH, WORK, SHEIGHT, B( ISTARTM,
        !           524:      $                IHI-NS ), LDB )
        !           525:       END IF
        !           526:       IF ( ILZ ) THEN
        !           527:          CALL DGEMM( 'N', 'N', N, NS+1, NS+1, ONE, Z( 1, IHI-NS ), LDZ,
        !           528:      $               ZC, LDZC, ZERO, WORK, N )
        !           529:          CALL DLACPY( 'ALL', N, NS+1, WORK, N, Z( 1, IHI-NS ), LDZ )
        !           530:       END IF
        !           531: 
        !           532:       END SUBROUTINE

CVSweb interface <joel.bertrand@systella.fr>