File:  [local] / rpl / lapack / lapack / dlaqz4.f
Revision 1.1: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:55:30 2023 UTC (8 months, 3 weeks 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 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>