File:  [local] / rpl / lapack / lapack / dlaqz3.f
Revision 1.1: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:55:29 2023 UTC (9 months, 1 week ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Ajout de fichiers de lapack 3.11

    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>