File:  [local] / rpl / lapack / lapack / dormbr.f
Revision 1.17: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:39:02 2023 UTC (8 months, 3 weeks 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 DORMBR
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *> \htmlonly
    9: *> Download DORMBR + dependencies
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dormbr.f">
   11: *> [TGZ]</a>
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dormbr.f">
   13: *> [ZIP]</a>
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dormbr.f">
   15: *> [TXT]</a>
   16: *> \endhtmlonly
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DORMBR( VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C,
   22: *                          LDC, WORK, LWORK, INFO )
   23: *
   24: *       .. Scalar Arguments ..
   25: *       CHARACTER          SIDE, TRANS, VECT
   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: *> If VECT = 'Q', DORMBR overwrites the general real M-by-N matrix C
   39: *> with
   40: *>                 SIDE = 'L'     SIDE = 'R'
   41: *> TRANS = 'N':      Q * C          C * Q
   42: *> TRANS = 'T':      Q**T * C       C * Q**T
   43: *>
   44: *> If VECT = 'P', DORMBR overwrites the general real M-by-N matrix C
   45: *> with
   46: *>                 SIDE = 'L'     SIDE = 'R'
   47: *> TRANS = 'N':      P * C          C * P
   48: *> TRANS = 'T':      P**T * C       C * P**T
   49: *>
   50: *> Here Q and P**T are the orthogonal matrices determined by DGEBRD when
   51: *> reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and
   52: *> P**T are defined as products of elementary reflectors H(i) and G(i)
   53: *> respectively.
   54: *>
   55: *> Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the
   56: *> order of the orthogonal matrix Q or P**T that is applied.
   57: *>
   58: *> If VECT = 'Q', A is assumed to have been an NQ-by-K matrix:
   59: *> if nq >= k, Q = H(1) H(2) . . . H(k);
   60: *> if nq < k, Q = H(1) H(2) . . . H(nq-1).
   61: *>
   62: *> If VECT = 'P', A is assumed to have been a K-by-NQ matrix:
   63: *> if k < nq, P = G(1) G(2) . . . G(k);
   64: *> if k >= nq, P = G(1) G(2) . . . G(nq-1).
   65: *> \endverbatim
   66: *
   67: *  Arguments:
   68: *  ==========
   69: *
   70: *> \param[in] VECT
   71: *> \verbatim
   72: *>          VECT is CHARACTER*1
   73: *>          = 'Q': apply Q or Q**T;
   74: *>          = 'P': apply P or P**T.
   75: *> \endverbatim
   76: *>
   77: *> \param[in] SIDE
   78: *> \verbatim
   79: *>          SIDE is CHARACTER*1
   80: *>          = 'L': apply Q, Q**T, P or P**T from the Left;
   81: *>          = 'R': apply Q, Q**T, P or P**T from the Right.
   82: *> \endverbatim
   83: *>
   84: *> \param[in] TRANS
   85: *> \verbatim
   86: *>          TRANS is CHARACTER*1
   87: *>          = 'N':  No transpose, apply Q  or P;
   88: *>          = 'T':  Transpose, apply Q**T or P**T.
   89: *> \endverbatim
   90: *>
   91: *> \param[in] M
   92: *> \verbatim
   93: *>          M is INTEGER
   94: *>          The number of rows of the matrix C. M >= 0.
   95: *> \endverbatim
   96: *>
   97: *> \param[in] N
   98: *> \verbatim
   99: *>          N is INTEGER
  100: *>          The number of columns of the matrix C. N >= 0.
  101: *> \endverbatim
  102: *>
  103: *> \param[in] K
  104: *> \verbatim
  105: *>          K is INTEGER
  106: *>          If VECT = 'Q', the number of columns in the original
  107: *>          matrix reduced by DGEBRD.
  108: *>          If VECT = 'P', the number of rows in the original
  109: *>          matrix reduced by DGEBRD.
  110: *>          K >= 0.
  111: *> \endverbatim
  112: *>
  113: *> \param[in] A
  114: *> \verbatim
  115: *>          A is DOUBLE PRECISION array, dimension
  116: *>                                (LDA,min(nq,K)) if VECT = 'Q'
  117: *>                                (LDA,nq)        if VECT = 'P'
  118: *>          The vectors which define the elementary reflectors H(i) and
  119: *>          G(i), whose products determine the matrices Q and P, as
  120: *>          returned by DGEBRD.
  121: *> \endverbatim
  122: *>
  123: *> \param[in] LDA
  124: *> \verbatim
  125: *>          LDA is INTEGER
  126: *>          The leading dimension of the array A.
  127: *>          If VECT = 'Q', LDA >= max(1,nq);
  128: *>          if VECT = 'P', LDA >= max(1,min(nq,K)).
  129: *> \endverbatim
  130: *>
  131: *> \param[in] TAU
  132: *> \verbatim
  133: *>          TAU is DOUBLE PRECISION array, dimension (min(nq,K))
  134: *>          TAU(i) must contain the scalar factor of the elementary
  135: *>          reflector H(i) or G(i) which determines Q or P, as returned
  136: *>          by DGEBRD in the array argument TAUQ or TAUP.
  137: *> \endverbatim
  138: *>
  139: *> \param[in,out] C
  140: *> \verbatim
  141: *>          C is DOUBLE PRECISION array, dimension (LDC,N)
  142: *>          On entry, the M-by-N matrix C.
  143: *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q
  144: *>          or P*C or P**T*C or C*P or C*P**T.
  145: *> \endverbatim
  146: *>
  147: *> \param[in] LDC
  148: *> \verbatim
  149: *>          LDC is INTEGER
  150: *>          The leading dimension of the array C. LDC >= max(1,M).
  151: *> \endverbatim
  152: *>
  153: *> \param[out] WORK
  154: *> \verbatim
  155: *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
  156: *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  157: *> \endverbatim
  158: *>
  159: *> \param[in] LWORK
  160: *> \verbatim
  161: *>          LWORK is INTEGER
  162: *>          The dimension of the array WORK.
  163: *>          If SIDE = 'L', LWORK >= max(1,N);
  164: *>          if SIDE = 'R', LWORK >= max(1,M).
  165: *>          For optimum performance LWORK >= N*NB if SIDE = 'L', and
  166: *>          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
  167: *>          blocksize.
  168: *>
  169: *>          If LWORK = -1, then a workspace query is assumed; the routine
  170: *>          only calculates the optimal size of the WORK array, returns
  171: *>          this value as the first entry of the WORK array, and no error
  172: *>          message related to LWORK is issued by XERBLA.
  173: *> \endverbatim
  174: *>
  175: *> \param[out] INFO
  176: *> \verbatim
  177: *>          INFO is INTEGER
  178: *>          = 0:  successful exit
  179: *>          < 0:  if INFO = -i, the i-th argument had an illegal value
  180: *> \endverbatim
  181: *
  182: *  Authors:
  183: *  ========
  184: *
  185: *> \author Univ. of Tennessee
  186: *> \author Univ. of California Berkeley
  187: *> \author Univ. of Colorado Denver
  188: *> \author NAG Ltd.
  189: *
  190: *> \ingroup doubleOTHERcomputational
  191: *
  192: *  =====================================================================
  193:       SUBROUTINE DORMBR( VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C,
  194:      $                   LDC, WORK, LWORK, INFO )
  195: *
  196: *  -- LAPACK computational routine --
  197: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  198: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  199: *
  200: *     .. Scalar Arguments ..
  201:       CHARACTER          SIDE, TRANS, VECT
  202:       INTEGER            INFO, K, LDA, LDC, LWORK, M, N
  203: *     ..
  204: *     .. Array Arguments ..
  205:       DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
  206: *     ..
  207: *
  208: *  =====================================================================
  209: *
  210: *     .. Local Scalars ..
  211:       LOGICAL            APPLYQ, LEFT, LQUERY, NOTRAN
  212:       CHARACTER          TRANST
  213:       INTEGER            I1, I2, IINFO, LWKOPT, MI, NB, NI, NQ, NW
  214: *     ..
  215: *     .. External Functions ..
  216:       LOGICAL            LSAME
  217:       INTEGER            ILAENV
  218:       EXTERNAL           LSAME, ILAENV
  219: *     ..
  220: *     .. External Subroutines ..
  221:       EXTERNAL           DORMLQ, DORMQR, XERBLA
  222: *     ..
  223: *     .. Intrinsic Functions ..
  224:       INTRINSIC          MAX, MIN
  225: *     ..
  226: *     .. Executable Statements ..
  227: *
  228: *     Test the input arguments
  229: *
  230:       INFO = 0
  231:       APPLYQ = LSAME( VECT, 'Q' )
  232:       LEFT = LSAME( SIDE, 'L' )
  233:       NOTRAN = LSAME( TRANS, 'N' )
  234:       LQUERY = ( LWORK.EQ.-1 )
  235: *
  236: *     NQ is the order of Q or P and NW is the minimum dimension of WORK
  237: *
  238:       IF( LEFT ) THEN
  239:          NQ = M
  240:          NW = MAX( 1, N )
  241:       ELSE
  242:          NQ = N
  243:          NW = MAX( 1, M )
  244:       END IF
  245:       IF( .NOT.APPLYQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN
  246:          INFO = -1
  247:       ELSE IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
  248:          INFO = -2
  249:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
  250:          INFO = -3
  251:       ELSE IF( M.LT.0 ) THEN
  252:          INFO = -4
  253:       ELSE IF( N.LT.0 ) THEN
  254:          INFO = -5
  255:       ELSE IF( K.LT.0 ) THEN
  256:          INFO = -6
  257:       ELSE IF( ( APPLYQ .AND. LDA.LT.MAX( 1, NQ ) ) .OR.
  258:      $         ( .NOT.APPLYQ .AND. LDA.LT.MAX( 1, MIN( NQ, K ) ) ) )
  259:      $          THEN
  260:          INFO = -8
  261:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
  262:          INFO = -11
  263:       ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
  264:          INFO = -13
  265:       END IF
  266: *
  267:       IF( INFO.EQ.0 ) THEN
  268:          IF( APPLYQ ) THEN
  269:             IF( LEFT ) THEN
  270:                NB = ILAENV( 1, 'DORMQR', SIDE // TRANS, M-1, N, M-1,
  271:      $              -1 )
  272:             ELSE
  273:                NB = ILAENV( 1, 'DORMQR', SIDE // TRANS, M, N-1, N-1,
  274:      $              -1 )
  275:             END IF
  276:          ELSE
  277:             IF( LEFT ) THEN
  278:                NB = ILAENV( 1, 'DORMLQ', SIDE // TRANS, M-1, N, M-1,
  279:      $              -1 )
  280:             ELSE
  281:                NB = ILAENV( 1, 'DORMLQ', SIDE // TRANS, M, N-1, N-1,
  282:      $              -1 )
  283:             END IF
  284:          END IF
  285:          LWKOPT = NW*NB
  286:          WORK( 1 ) = LWKOPT
  287:       END IF
  288: *
  289:       IF( INFO.NE.0 ) THEN
  290:          CALL XERBLA( 'DORMBR', -INFO )
  291:          RETURN
  292:       ELSE IF( LQUERY ) THEN
  293:          RETURN
  294:       END IF
  295: *
  296: *     Quick return if possible
  297: *
  298:       WORK( 1 ) = 1
  299:       IF( M.EQ.0 .OR. N.EQ.0 )
  300:      $   RETURN
  301: *
  302:       IF( APPLYQ ) THEN
  303: *
  304: *        Apply Q
  305: *
  306:          IF( NQ.GE.K ) THEN
  307: *
  308: *           Q was determined by a call to DGEBRD with nq >= k
  309: *
  310:             CALL DORMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
  311:      $                   WORK, LWORK, IINFO )
  312:          ELSE IF( NQ.GT.1 ) THEN
  313: *
  314: *           Q was determined by a call to DGEBRD with nq < k
  315: *
  316:             IF( LEFT ) THEN
  317:                MI = M - 1
  318:                NI = N
  319:                I1 = 2
  320:                I2 = 1
  321:             ELSE
  322:                MI = M
  323:                NI = N - 1
  324:                I1 = 1
  325:                I2 = 2
  326:             END IF
  327:             CALL DORMQR( SIDE, TRANS, MI, NI, NQ-1, A( 2, 1 ), LDA, TAU,
  328:      $                   C( I1, I2 ), LDC, WORK, LWORK, IINFO )
  329:          END IF
  330:       ELSE
  331: *
  332: *        Apply P
  333: *
  334:          IF( NOTRAN ) THEN
  335:             TRANST = 'T'
  336:          ELSE
  337:             TRANST = 'N'
  338:          END IF
  339:          IF( NQ.GT.K ) THEN
  340: *
  341: *           P was determined by a call to DGEBRD with nq > k
  342: *
  343:             CALL DORMLQ( SIDE, TRANST, M, N, K, A, LDA, TAU, C, LDC,
  344:      $                   WORK, LWORK, IINFO )
  345:          ELSE IF( NQ.GT.1 ) THEN
  346: *
  347: *           P was determined by a call to DGEBRD with nq <= k
  348: *
  349:             IF( LEFT ) THEN
  350:                MI = M - 1
  351:                NI = N
  352:                I1 = 2
  353:                I2 = 1
  354:             ELSE
  355:                MI = M
  356:                NI = N - 1
  357:                I1 = 1
  358:                I2 = 2
  359:             END IF
  360:             CALL DORMLQ( SIDE, TRANST, MI, NI, NQ-1, A( 1, 2 ), LDA,
  361:      $                   TAU, C( I1, I2 ), LDC, WORK, LWORK, IINFO )
  362:          END IF
  363:       END IF
  364:       WORK( 1 ) = LWKOPT
  365:       RETURN
  366: *
  367: *     End of DORMBR
  368: *
  369:       END

CVSweb interface <joel.bertrand@systella.fr>