File:  [local] / rpl / lapack / lapack / dormrq.f
Revision 1.20: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:03 2023 UTC (9 months, 1 week ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    1: *> \brief \b DORMRQ
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DORMRQ + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dormrq.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dormrq.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dormrq.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DORMRQ( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
   22: *                          WORK, LWORK, INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          SIDE, TRANS
   26: *       INTEGER            INFO, K, LDA, LDC, LWORK, M, N
   27: *       ..
   28: *       .. Array Arguments ..
   29: *       DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
   30: *       ..
   31: *
   32: *
   33: *> \par Purpose:
   34: *  =============
   35: *>
   36: *> \verbatim
   37: *>
   38: *> DORMRQ overwrites the general real M-by-N matrix C with
   39: *>
   40: *>                 SIDE = 'L'     SIDE = 'R'
   41: *> TRANS = 'N':      Q * C          C * Q
   42: *> TRANS = 'T':      Q**T * C       C * Q**T
   43: *>
   44: *> where Q is a real orthogonal matrix defined as the product of k
   45: *> elementary reflectors
   46: *>
   47: *>       Q = H(1) H(2) . . . H(k)
   48: *>
   49: *> as returned by DGERQF. Q is of order M if SIDE = 'L' and of order N
   50: *> if SIDE = 'R'.
   51: *> \endverbatim
   52: *
   53: *  Arguments:
   54: *  ==========
   55: *
   56: *> \param[in] SIDE
   57: *> \verbatim
   58: *>          SIDE is CHARACTER*1
   59: *>          = 'L': apply Q or Q**T from the Left;
   60: *>          = 'R': apply Q or Q**T from the Right.
   61: *> \endverbatim
   62: *>
   63: *> \param[in] TRANS
   64: *> \verbatim
   65: *>          TRANS is CHARACTER*1
   66: *>          = 'N':  No transpose, apply Q;
   67: *>          = 'T':  Transpose, apply Q**T.
   68: *> \endverbatim
   69: *>
   70: *> \param[in] M
   71: *> \verbatim
   72: *>          M is INTEGER
   73: *>          The number of rows of the matrix C. M >= 0.
   74: *> \endverbatim
   75: *>
   76: *> \param[in] N
   77: *> \verbatim
   78: *>          N is INTEGER
   79: *>          The number of columns of the matrix C. N >= 0.
   80: *> \endverbatim
   81: *>
   82: *> \param[in] K
   83: *> \verbatim
   84: *>          K is INTEGER
   85: *>          The number of elementary reflectors whose product defines
   86: *>          the matrix Q.
   87: *>          If SIDE = 'L', M >= K >= 0;
   88: *>          if SIDE = 'R', N >= K >= 0.
   89: *> \endverbatim
   90: *>
   91: *> \param[in] A
   92: *> \verbatim
   93: *>          A is DOUBLE PRECISION array, dimension
   94: *>                               (LDA,M) if SIDE = 'L',
   95: *>                               (LDA,N) if SIDE = 'R'
   96: *>          The i-th row must contain the vector which defines the
   97: *>          elementary reflector H(i), for i = 1,2,...,k, as returned by
   98: *>          DGERQF in the last k rows of its array argument A.
   99: *> \endverbatim
  100: *>
  101: *> \param[in] LDA
  102: *> \verbatim
  103: *>          LDA is INTEGER
  104: *>          The leading dimension of the array A. LDA >= max(1,K).
  105: *> \endverbatim
  106: *>
  107: *> \param[in] TAU
  108: *> \verbatim
  109: *>          TAU is DOUBLE PRECISION array, dimension (K)
  110: *>          TAU(i) must contain the scalar factor of the elementary
  111: *>          reflector H(i), as returned by DGERQF.
  112: *> \endverbatim
  113: *>
  114: *> \param[in,out] C
  115: *> \verbatim
  116: *>          C is DOUBLE PRECISION array, dimension (LDC,N)
  117: *>          On entry, the M-by-N matrix C.
  118: *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
  119: *> \endverbatim
  120: *>
  121: *> \param[in] LDC
  122: *> \verbatim
  123: *>          LDC is INTEGER
  124: *>          The leading dimension of the array C. LDC >= max(1,M).
  125: *> \endverbatim
  126: *>
  127: *> \param[out] WORK
  128: *> \verbatim
  129: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  130: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  131: *> \endverbatim
  132: *>
  133: *> \param[in] LWORK
  134: *> \verbatim
  135: *>          LWORK is INTEGER
  136: *>          The dimension of the array WORK.
  137: *>          If SIDE = 'L', LWORK >= max(1,N);
  138: *>          if SIDE = 'R', LWORK >= max(1,M).
  139: *>          For good performance, LWORK should generally be larger.
  140: *>
  141: *>          If LWORK = -1, then a workspace query is assumed; the routine
  142: *>          only calculates the optimal size of the WORK array, returns
  143: *>          this value as the first entry of the WORK array, and no error
  144: *>          message related to LWORK is issued by XERBLA.
  145: *> \endverbatim
  146: *>
  147: *> \param[out] INFO
  148: *> \verbatim
  149: *>          INFO is INTEGER
  150: *>          = 0:  successful exit
  151: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  152: *> \endverbatim
  153: *
  154: *  Authors:
  155: *  ========
  156: *
  157: *> \author Univ. of Tennessee
  158: *> \author Univ. of California Berkeley
  159: *> \author Univ. of Colorado Denver
  160: *> \author NAG Ltd.
  161: *
  162: *> \ingroup doubleOTHERcomputational
  163: *
  164: *  =====================================================================
  165:       SUBROUTINE DORMRQ( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
  166:      $                   WORK, LWORK, INFO )
  167: *
  168: *  -- LAPACK computational routine --
  169: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  170: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  171: *
  172: *     .. Scalar Arguments ..
  173:       CHARACTER          SIDE, TRANS
  174:       INTEGER            INFO, K, LDA, LDC, LWORK, M, N
  175: *     ..
  176: *     .. Array Arguments ..
  177:       DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
  178: *     ..
  179: *
  180: *  =====================================================================
  181: *
  182: *     .. Parameters ..
  183:       INTEGER            NBMAX, LDT, TSIZE
  184:       PARAMETER          ( NBMAX = 64, LDT = NBMAX+1,
  185:      $                     TSIZE = LDT*NBMAX )
  186: *     ..
  187: *     .. Local Scalars ..
  188:       LOGICAL            LEFT, LQUERY, NOTRAN
  189:       CHARACTER          TRANST
  190:       INTEGER            I, I1, I2, I3, IB, IINFO, IWT, LDWORK, LWKOPT,
  191:      $                   MI, NB, NBMIN, NI, NQ, NW
  192: *     ..
  193: *     .. External Functions ..
  194:       LOGICAL            LSAME
  195:       INTEGER            ILAENV
  196:       EXTERNAL           LSAME, ILAENV
  197: *     ..
  198: *     .. External Subroutines ..
  199:       EXTERNAL           DLARFB, DLARFT, DORMR2, XERBLA
  200: *     ..
  201: *     .. Intrinsic Functions ..
  202:       INTRINSIC          MAX, MIN
  203: *     ..
  204: *     .. Executable Statements ..
  205: *
  206: *     Test the input arguments
  207: *
  208:       INFO = 0
  209:       LEFT = LSAME( SIDE, 'L' )
  210:       NOTRAN = LSAME( TRANS, 'N' )
  211:       LQUERY = ( LWORK.EQ.-1 )
  212: *
  213: *     NQ is the order of Q and NW is the minimum dimension of WORK
  214: *
  215:       IF( LEFT ) THEN
  216:          NQ = M
  217:          NW = MAX( 1, N )
  218:       ELSE
  219:          NQ = N
  220:          NW = MAX( 1, M )
  221:       END IF
  222:       IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
  223:          INFO = -1
  224:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
  225:          INFO = -2
  226:       ELSE IF( M.LT.0 ) THEN
  227:          INFO = -3
  228:       ELSE IF( N.LT.0 ) THEN
  229:          INFO = -4
  230:       ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
  231:          INFO = -5
  232:       ELSE IF( LDA.LT.MAX( 1, K ) ) THEN
  233:          INFO = -7
  234:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
  235:          INFO = -10
  236:       ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
  237:          INFO = -12
  238:       END IF
  239: *
  240:       IF( INFO.EQ.0 ) THEN
  241: *
  242: *        Compute the workspace requirements
  243: *
  244:          IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  245:             LWKOPT = 1
  246:          ELSE
  247:             NB = MIN( NBMAX, ILAENV( 1, 'DORMRQ', SIDE // TRANS, M, N,
  248:      $                               K, -1 ) )
  249:             LWKOPT = NW*NB + TSIZE
  250:          END IF
  251:          WORK( 1 ) = LWKOPT
  252:       END IF
  253: *
  254:       IF( INFO.NE.0 ) THEN
  255:          CALL XERBLA( 'DORMRQ', -INFO )
  256:          RETURN
  257:       ELSE IF( LQUERY ) THEN
  258:          RETURN
  259:       END IF
  260: *
  261: *     Quick return if possible
  262: *
  263:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  264:          RETURN
  265:       END IF
  266: *
  267:       NBMIN = 2
  268:       LDWORK = NW
  269:       IF( NB.GT.1 .AND. NB.LT.K ) THEN
  270:          IF( LWORK.LT.LWKOPT ) THEN
  271:             NB = (LWORK-TSIZE) / LDWORK
  272:             NBMIN = MAX( 2, ILAENV( 2, 'DORMRQ', SIDE // TRANS, M, N, K,
  273:      $              -1 ) )
  274:          END IF
  275:       END IF
  276: *
  277:       IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN
  278: *
  279: *        Use unblocked code
  280: *
  281:          CALL DORMR2( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK,
  282:      $                IINFO )
  283:       ELSE
  284: *
  285: *        Use blocked code
  286: *
  287:          IWT = 1 + NW*NB
  288:          IF( ( LEFT .AND. .NOT.NOTRAN ) .OR.
  289:      $       ( .NOT.LEFT .AND. NOTRAN ) ) THEN
  290:             I1 = 1
  291:             I2 = K
  292:             I3 = NB
  293:          ELSE
  294:             I1 = ( ( K-1 ) / NB )*NB + 1
  295:             I2 = 1
  296:             I3 = -NB
  297:          END IF
  298: *
  299:          IF( LEFT ) THEN
  300:             NI = N
  301:          ELSE
  302:             MI = M
  303:          END IF
  304: *
  305:          IF( NOTRAN ) THEN
  306:             TRANST = 'T'
  307:          ELSE
  308:             TRANST = 'N'
  309:          END IF
  310: *
  311:          DO 10 I = I1, I2, I3
  312:             IB = MIN( NB, K-I+1 )
  313: *
  314: *           Form the triangular factor of the block reflector
  315: *           H = H(i+ib-1) . . . H(i+1) H(i)
  316: *
  317:             CALL DLARFT( 'Backward', 'Rowwise', NQ-K+I+IB-1, IB,
  318:      $                   A( I, 1 ), LDA, TAU( I ), WORK( IWT ), LDT )
  319:             IF( LEFT ) THEN
  320: *
  321: *              H or H**T is applied to C(1:m-k+i+ib-1,1:n)
  322: *
  323:                MI = M - K + I + IB - 1
  324:             ELSE
  325: *
  326: *              H or H**T is applied to C(1:m,1:n-k+i+ib-1)
  327: *
  328:                NI = N - K + I + IB - 1
  329:             END IF
  330: *
  331: *           Apply H or H**T
  332: *
  333:             CALL DLARFB( SIDE, TRANST, 'Backward', 'Rowwise', MI, NI,
  334:      $                   IB, A( I, 1 ), LDA, WORK( IWT ), LDT, C, LDC,
  335:      $                   WORK, LDWORK )
  336:    10    CONTINUE
  337:       END IF
  338:       WORK( 1 ) = LWKOPT
  339:       RETURN
  340: *
  341: *     End of DORMRQ
  342: *
  343:       END

CVSweb interface <joel.bertrand@systella.fr>