Return to zlamswlq.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *
    2: *  Definition:
    3: *  ===========
    4: *
    5: *      SUBROUTINE ZLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T,
    6: *     $                LDT, C, LDC, WORK, LWORK, INFO )
    7: *
    8: *
    9: *     .. Scalar Arguments ..
   10: *      CHARACTER         SIDE, TRANS
   11: *      INTEGER           INFO, LDA, M, N, K, MB, NB, LDT, LWORK, LDC
   12: *     ..
   13: *     .. Array Arguments ..
   14: *      COMPLEX*16        A( LDA, * ), WORK( * ), C(LDC, * ),
   15: *     $                  T( LDT, * )
   16: *> \par Purpose:
   17: *  =============
   18: *>
   19: *> \verbatim
   20: *>
   21: *>    ZLAMQRTS overwrites the general real M-by-N matrix C with
   22: *>
   23: *>
   24: *>                    SIDE = 'L'     SIDE = 'R'
   25: *>    TRANS = 'N':      Q * C          C * Q
   26: *>    TRANS = 'T':      Q**T * C       C * Q**T
   27: *>    where Q is a real orthogonal matrix defined as the product of blocked
   28: *>    elementary reflectors computed by short wide LQ
   29: *>    factorization (ZLASWLQ)
   30: *> \endverbatim
   31: *
   32: *  Arguments:
   33: *  ==========
   34: *
   35: *> \param[in] SIDE
   36: *> \verbatim
   37: *>          SIDE is CHARACTER*1
   38: *>          = 'L': apply Q or Q**T from the Left;
   39: *>          = 'R': apply Q or Q**T from the Right.
   40: *> \endverbatim
   41: *>
   42: *> \param[in] TRANS
   43: *> \verbatim
   44: *>          TRANS is CHARACTER*1
   45: *>          = 'N':  No transpose, apply Q;
   46: *>          = 'T':  Transpose, apply Q**T.
   47: *> \endverbatim
   48: *>
   49: *> \param[in] M
   50: *> \verbatim
   51: *>          M is INTEGER
   52: *>          The number of rows of the matrix A.  M >=0.
   53: *> \endverbatim
   54: *>
   55: *> \param[in] N
   56: *> \verbatim
   57: *>          N is INTEGER
   58: *>          The number of columns of the matrix C. N >= M.
   59: *> \endverbatim
   60: *>
   61: *> \param[in] K
   62: *> \verbatim
   63: *>          K is INTEGER
   64: *>          The number of elementary reflectors whose product defines
   65: *>          the matrix Q.
   66: *>          M >= K >= 0;
   67: *>
   68: *> \endverbatim
   69: *> \param[in] MB
   70: *> \verbatim
   71: *>          MB is INTEGER
   72: *>          The row block size to be used in the blocked QR.
   73: *>          M >= MB >= 1
   74: *> \endverbatim
   75: *>
   76: *> \param[in] NB
   77: *> \verbatim
   78: *>          NB is INTEGER
   79: *>          The column block size to be used in the blocked QR.
   80: *>          NB > M.
   81: *> \endverbatim
   82: *>
   83: *> \param[in] NB
   84: *> \verbatim
   85: *>          NB is INTEGER
   86: *>          The block size to be used in the blocked QR.
   87: *>                MB > M.
   88: *>
   89: *> \endverbatim
   90: *>
   91: *> \param[in,out] A
   92: *> \verbatim
   93: *>          A is COMPLEX*16 array, dimension (LDA,K)
   94: *>          The i-th row must contain the vector which defines the blocked
   95: *>          elementary reflector H(i), for i = 1,2,...,k, as returned by
   96: *>          DLASWLQ in the first k rows of its array argument A.
   97: *> \endverbatim
   98: *>
   99: *> \param[in] LDA
  100: *> \verbatim
  101: *>          LDA is INTEGER
  102: *>          The leading dimension of the array A.
  103: *>          If SIDE = 'L', LDA >= max(1,M);
  104: *>          if SIDE = 'R', LDA >= max(1,N).
  105: *> \endverbatim
  106: *>
  107: *> \param[in] T
  108: *> \verbatim
  109: *>          T is COMPLEX*16 array, dimension
  110: *>          ( M * Number of blocks(CEIL(N-K/NB-K)),
  111: *>          The blocked upper triangular block reflectors stored in compact form
  112: *>          as a sequence of upper triangular blocks.  See below
  113: *>          for further details.
  114: *> \endverbatim
  115: *>
  116: *> \param[in] LDT
  117: *> \verbatim
  118: *>          LDT is INTEGER
  119: *>          The leading dimension of the array T.  LDT >= MB.
  120: *> \endverbatim
  121: *>
  122: *> \param[in,out] C
  123: *> \verbatim
  124: *>          C is COMPLEX*16 array, dimension (LDC,N)
  125: *>          On entry, the M-by-N matrix C.
  126: *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
  127: *> \endverbatim
  128: *>
  129: *> \param[in] LDC
  130: *> \verbatim
  131: *>          LDC is INTEGER
  132: *>          The leading dimension of the array C. LDC >= max(1,M).
  133: *> \endverbatim
  134: *>
  135: *> \param[out] WORK
  136: *> \verbatim
  137: *>         (workspace) COMPLEX*16 array, dimension (MAX(1,LWORK))
  138: *> \endverbatim
  139: *>
  140: *> \param[in] LWORK
  141: *> \verbatim
  142: *>          LWORK is INTEGER
  143: *>          The dimension of the array WORK.
  144: *>          If SIDE = 'L', LWORK >= max(1,NB) * MB;
  145: *>          if SIDE = 'R', LWORK >= max(1,M) * MB.
  146: *>          If LWORK = -1, then a workspace query is assumed; the routine
  147: *>          only calculates the optimal size of the WORK array, returns
  148: *>          this value as the first entry of the WORK array, and no error
  149: *>          message related to LWORK is issued by XERBLA.
  150: *> \endverbatim
  151: *>
  152: *> \param[out] INFO
  153: *> \verbatim
  154: *>          INFO is INTEGER
  155: *>          = 0:  successful exit
  156: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  157: *> \endverbatim
  158: *
  159: *  Authors:
  160: *  ========
  161: *
  162: *> \author Univ. of Tennessee
  163: *> \author Univ. of California Berkeley
  164: *> \author Univ. of Colorado Denver
  165: *> \author NAG Ltd.
  166: *
  167: *> \par Further Details:
  168: *  =====================
  169: *>
  170: *> \verbatim
  171: *> Short-Wide LQ (SWLQ) performs LQ by a sequence of orthogonal transformations,
  172: *> representing Q as a product of other orthogonal matrices
  173: *>   Q = Q(1) * Q(2) * . . . * Q(k)
  174: *> where each Q(i) zeros out upper diagonal entries of a block of NB rows of A:
  175: *>   Q(1) zeros out the upper diagonal entries of rows 1:NB of A
  176: *>   Q(2) zeros out the bottom MB-N rows of rows [1:M,NB+1:2*NB-M] of A
  177: *>   Q(3) zeros out the bottom MB-N rows of rows [1:M,2*NB-M+1:3*NB-2*M] of A
  178: *>   . . .
  179: *>
  180: *> Q(1) is computed by GELQT, which represents Q(1) by Householder vectors
  181: *> stored under the diagonal of rows 1:MB of A, and by upper triangular
  182: *> block reflectors, stored in array T(1:LDT,1:N).
  183: *> For more information see Further Details in GELQT.
  184: *>
  185: *> Q(i) for i>1 is computed by TPLQT, which represents Q(i) by Householder vectors
  186: *> stored in columns [(i-1)*(NB-M)+M+1:i*(NB-M)+M] of A, and by upper triangular
  187: *> block reflectors, stored in array T(1:LDT,(i-1)*M+1:i*M).
  188: *> The last Q(k) may use fewer rows.
  189: *> For more information see Further Details in TPQRT.
  190: *>
  191: *> For more details of the overall algorithm, see the description of
  192: *> Sequential TSQR in Section 2.2 of [1].
  193: *>
  194: *> [1] “Communication-Optimal Parallel and Sequential QR and LU Factorizations,”
  195: *>     J. Demmel, L. Grigori, M. Hoemmen, J. Langou,
  196: *>     SIAM J. Sci. Comput, vol. 34, no. 1, 2012
  197: *> \endverbatim
  198: *>
  199: *  =====================================================================
  200:       SUBROUTINE ZLAMSWLQ( SIDE, TRANS, M, N, K, MB, NB, A, LDA, T,
  201:      $    LDT, C, LDC, WORK, LWORK, INFO )
  202: *
  203: *  -- LAPACK computational routine (version 3.7.0) --
  204: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  205: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  206: *     December 2016
  207: *
  208: *     .. Scalar Arguments ..
  209:       CHARACTER         SIDE, TRANS
  210:       INTEGER           INFO, LDA, M, N, K, MB, NB, LDT, LWORK, LDC, LW
  211: *     ..
  212: *     .. Array Arguments ..
  213:       COMPLEX*16        A( LDA, * ), WORK( * ), C(LDC, * ),
  214:      $      T( LDT, * )
  215: *     ..
  216: *
  217: * =====================================================================
  218: *
  219: *     ..
  220: *     .. Local Scalars ..
  221:       LOGICAL    LEFT, RIGHT, TRAN, NOTRAN, LQUERY
  222:       INTEGER    I, II, KK, CTR
  223: *     ..
  224: *     .. External Functions ..
  225:       LOGICAL            LSAME
  226:       EXTERNAL           LSAME
  227: *     .. External Subroutines ..
  228:       EXTERNAL    ZTPMLQT, ZGEMLQT, XERBLA
  229: *     ..
  230: *     .. Executable Statements ..
  231: *
  232: *     Test the input arguments
  233: *
  234:       LQUERY  = LWORK.LT.0
  235:       NOTRAN  = LSAME( TRANS, 'N' )
  236:       TRAN    = LSAME( TRANS, 'C' )
  237:       LEFT    = LSAME( SIDE, 'L' )
  238:       RIGHT   = LSAME( SIDE, 'R' )
  239:       IF (LEFT) THEN
  240:         LW = N * MB
  241:       ELSE
  242:         LW = M * MB
  243:       END IF
  244: *
  245:       INFO = 0
  246:       IF( .NOT.LEFT .AND. .NOT.RIGHT ) THEN
  247:          INFO = -1
  248:       ELSE IF( .NOT.TRAN .AND. .NOT.NOTRAN ) THEN
  249:          INFO = -2
  250:       ELSE IF( M.LT.0 ) THEN
  251:         INFO = -3
  252:       ELSE IF( N.LT.0 ) THEN
  253:         INFO = -4
  254:       ELSE IF( K.LT.0 ) THEN
  255:         INFO = -5
  256:       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
  257:         INFO = -9
  258:       ELSE IF( LDT.LT.MAX( 1, MB) ) THEN
  259:         INFO = -11
  260:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
  261:          INFO = -13
  262:       ELSE IF(( LWORK.LT.MAX(1,LW)).AND.(.NOT.LQUERY)) THEN
  263:         INFO = -15
  264:       END IF
  265: *
  266:       IF( INFO.NE.0 ) THEN
  267:         CALL XERBLA( 'ZLAMSWLQ', -INFO )
  268:         WORK(1) = LW
  269:         RETURN
  270:       ELSE IF (LQUERY) THEN
  271:         WORK(1) = LW
  272:         RETURN
  273:       END IF
  274: *
  275: *     Quick return if possible
  276: *
  277:       IF( MIN(M,N,K).EQ.0 ) THEN
  278:         RETURN
  279:       END IF
  280: *
  281:       IF((NB.LE.K).OR.(NB.GE.MAX(M,N,K))) THEN
  282:         CALL ZGEMLQT( SIDE, TRANS, M, N, K, MB, A, LDA,
  283:      $        T, LDT, C, LDC, WORK, INFO)
  284:         RETURN
  285:       END IF
  286: *
  287:       IF(LEFT.AND.TRAN) THEN
  288: *
  289: *         Multiply Q to the last block of C
  290: *
  291:           KK = MOD((M-K),(NB-K))
  292:           CTR = (M-K)/(NB-K)
  293: *
  294:           IF (KK.GT.0) THEN
  295:             II=M-KK+1
  296:             CALL ZTPMLQT('L','C',KK , N, K, 0, MB, A(1,II), LDA,
  297:      $        T(1,CTR*K+1), LDT, C(1,1), LDC,
  298:      $        C(II,1), LDC, WORK, INFO )
  299:           ELSE
  300:             II=M+1
  301:           END IF
  302: *
  303:           DO I=II-(NB-K),NB+1,-(NB-K)
  304: *
  305: *         Multiply Q to the current block of C (1:M,I:I+NB)
  306: *
  307:             CTR = CTR - 1
  308:             CALL ZTPMLQT('L','C',NB-K , N, K, 0,MB, A(1,I), LDA,
  309:      $          T(1,CTR*K+1),LDT, C(1,1), LDC,
  310:      $          C(I,1), LDC, WORK, INFO )
  311: 
  312:           END DO
  313: *
  314: *         Multiply Q to the first block of C (1:M,1:NB)
  315: *
  316:           CALL ZGEMLQT('L','C',NB , N, K, MB, A(1,1), LDA, T
  317:      $              ,LDT ,C(1,1), LDC, WORK, INFO )
  318: *
  319:       ELSE IF (LEFT.AND.NOTRAN) THEN
  320: *
  321: *         Multiply Q to the first block of C
  322: *
  323:          KK = MOD((M-K),(NB-K))
  324:          II=M-KK+1
  325:          CTR = 1
  326:          CALL ZGEMLQT('L','N',NB , N, K, MB, A(1,1), LDA, T
  327:      $              ,LDT ,C(1,1), LDC, WORK, INFO )
  328: *
  329:          DO I=NB+1,II-NB+K,(NB-K)
  330: *
  331: *         Multiply Q to the current block of C (I:I+NB,1:N)
  332: *
  333:           CALL ZTPMLQT('L','N',NB-K , N, K, 0,MB, A(1,I), LDA,
  334:      $         T(1, CTR * K + 1), LDT, C(1,1), LDC,
  335:      $         C(I,1), LDC, WORK, INFO )
  336:           CTR = CTR + 1
  337: *
  338:          END DO
  339:          IF(II.LE.M) THEN
  340: *
  341: *         Multiply Q to the last block of C
  342: *
  343:           CALL ZTPMLQT('L','N',KK , N, K, 0, MB, A(1,II), LDA,
  344:      $        T(1, CTR * K + 1), LDT, C(1,1), LDC,
  345:      $        C(II,1), LDC, WORK, INFO )
  346: *
  347:          END IF
  348: *
  349:       ELSE IF(RIGHT.AND.NOTRAN) THEN
  350: *
  351: *         Multiply Q to the last block of C
  352: *
  353:           KK = MOD((N-K),(NB-K))
  354:           CTR = (N-K)/(NB-K)
  355:           IF (KK.GT.0) THEN
  356:             II=N-KK+1
  357:             CALL ZTPMLQT('R','N',M , KK, K, 0, MB, A(1, II), LDA,
  358:      $        T(1, CTR * K + 1), LDT, C(1,1), LDC,
  359:      $        C(1,II), LDC, WORK, INFO )
  360:           ELSE
  361:             II=N+1
  362:           END IF
  363: *
  364:           DO I=II-(NB-K),NB+1,-(NB-K)
  365: *
  366: *         Multiply Q to the current block of C (1:M,I:I+MB)
  367: *
  368:           CTR = CTR - 1
  369:           CALL ZTPMLQT('R','N', M, NB-K, K, 0, MB, A(1, I), LDA,
  370:      $        T(1, CTR * K + 1), LDT, C(1,1), LDC,
  371:      $        C(1,I), LDC, WORK, INFO )
  372: 
  373:           END DO
  374: *
  375: *         Multiply Q to the first block of C (1:M,1:MB)
  376: *
  377:           CALL ZGEMLQT('R','N',M , NB, K, MB, A(1,1), LDA, T
  378:      $            ,LDT ,C(1,1), LDC, WORK, INFO )
  379: *
  380:       ELSE IF (RIGHT.AND.TRAN) THEN
  381: *
  382: *       Multiply Q to the first block of C
  383: *
  384:          KK = MOD((N-K),(NB-K))
  385:          II=N-KK+1
  386:          CALL ZGEMLQT('R','C',M , NB, K, MB, A(1,1), LDA, T
  387:      $            ,LDT ,C(1,1), LDC, WORK, INFO )
  388:          CTR = 1
  389: *
  390:          DO I=NB+1,II-NB+K,(NB-K)
  391: *
  392: *         Multiply Q to the current block of C (1:M,I:I+MB)
  393: *
  394:           CALL ZTPMLQT('R','C',M , NB-K, K, 0,MB, A(1,I), LDA,
  395:      $       T(1,CTR *K+1), LDT, C(1,1), LDC,
  396:      $       C(1,I), LDC, WORK, INFO )
  397:           CTR = CTR + 1
  398: *
  399:          END DO
  400:          IF(II.LE.N) THEN
  401: *
  402: *       Multiply Q to the last block of C
  403: *
  404:           CALL ZTPMLQT('R','C',M , KK, K, 0,MB, A(1,II), LDA,
  405:      $      T(1, CTR * K + 1),LDT, C(1,1), LDC,
  406:      $      C(1,II), LDC, WORK, INFO )
  407: *
  408:          END IF
  409: *
  410:       END IF
  411: *
  412:       WORK(1) = LW
  413:       RETURN
  414: *
  415: *     End of ZLAMSWLQ
  416: *
  417:       END