[BACK]Return to dormql.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dormql.f between versions 1.6 and 1.18

version 1.6, 2010/08/13 21:03:54 version 1.18, 2017/06/17 11:06:29
Line 1 Line 1
   *> \brief \b DORMQL
   *
   *  =========== DOCUMENTATION ===========
   *
   * Online html documentation available at
   *            http://www.netlib.org/lapack/explore-html/
   *
   *> \htmlonly
   *> Download DORMQL + dependencies
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dormql.f">
   *> [TGZ]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dormql.f">
   *> [ZIP]</a>
   *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dormql.f">
   *> [TXT]</a>
   *> \endhtmlonly
   *
   *  Definition:
   *  ===========
   *
   *       SUBROUTINE DORMQL( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
   *                          WORK, LWORK, INFO )
   *
   *       .. Scalar Arguments ..
   *       CHARACTER          SIDE, TRANS
   *       INTEGER            INFO, K, LDA, LDC, LWORK, M, N
   *       ..
   *       .. Array Arguments ..
   *       DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
   *       ..
   *
   *
   *> \par Purpose:
   *  =============
   *>
   *> \verbatim
   *>
   *> DORMQL overwrites the general real M-by-N matrix C with
   *>
   *>                 SIDE = 'L'     SIDE = 'R'
   *> TRANS = 'N':      Q * C          C * Q
   *> TRANS = 'T':      Q**T * C       C * Q**T
   *>
   *> where Q is a real orthogonal matrix defined as the product of k
   *> elementary reflectors
   *>
   *>       Q = H(k) . . . H(2) H(1)
   *>
   *> as returned by DGEQLF. Q is of order M if SIDE = 'L' and of order N
   *> if SIDE = 'R'.
   *> \endverbatim
   *
   *  Arguments:
   *  ==========
   *
   *> \param[in] SIDE
   *> \verbatim
   *>          SIDE is CHARACTER*1
   *>          = 'L': apply Q or Q**T from the Left;
   *>          = 'R': apply Q or Q**T from the Right.
   *> \endverbatim
   *>
   *> \param[in] TRANS
   *> \verbatim
   *>          TRANS is CHARACTER*1
   *>          = 'N':  No transpose, apply Q;
   *>          = 'T':  Transpose, apply Q**T.
   *> \endverbatim
   *>
   *> \param[in] M
   *> \verbatim
   *>          M is INTEGER
   *>          The number of rows of the matrix C. M >= 0.
   *> \endverbatim
   *>
   *> \param[in] N
   *> \verbatim
   *>          N is INTEGER
   *>          The number of columns of the matrix C. N >= 0.
   *> \endverbatim
   *>
   *> \param[in] K
   *> \verbatim
   *>          K is INTEGER
   *>          The number of elementary reflectors whose product defines
   *>          the matrix Q.
   *>          If SIDE = 'L', M >= K >= 0;
   *>          if SIDE = 'R', N >= K >= 0.
   *> \endverbatim
   *>
   *> \param[in] A
   *> \verbatim
   *>          A is DOUBLE PRECISION array, dimension (LDA,K)
   *>          The i-th column must contain the vector which defines the
   *>          elementary reflector H(i), for i = 1,2,...,k, as returned by
   *>          DGEQLF in the last k columns of its array argument A.
   *> \endverbatim
   *>
   *> \param[in] LDA
   *> \verbatim
   *>          LDA is INTEGER
   *>          The leading dimension of the array A.
   *>          If SIDE = 'L', LDA >= max(1,M);
   *>          if SIDE = 'R', LDA >= max(1,N).
   *> \endverbatim
   *>
   *> \param[in] TAU
   *> \verbatim
   *>          TAU is DOUBLE PRECISION array, dimension (K)
   *>          TAU(i) must contain the scalar factor of the elementary
   *>          reflector H(i), as returned by DGEQLF.
   *> \endverbatim
   *>
   *> \param[in,out] C
   *> \verbatim
   *>          C is DOUBLE PRECISION array, dimension (LDC,N)
   *>          On entry, the M-by-N matrix C.
   *>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
   *> \endverbatim
   *>
   *> \param[in] LDC
   *> \verbatim
   *>          LDC is INTEGER
   *>          The leading dimension of the array C. LDC >= max(1,M).
   *> \endverbatim
   *>
   *> \param[out] WORK
   *> \verbatim
   *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
   *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   *> \endverbatim
   *>
   *> \param[in] LWORK
   *> \verbatim
   *>          LWORK is INTEGER
   *>          The dimension of the array WORK.
   *>          If SIDE = 'L', LWORK >= max(1,N);
   *>          if SIDE = 'R', LWORK >= max(1,M).
   *>          For good performance, LWORK should generally be larger.
   *>
   *>          If LWORK = -1, then a workspace query is assumed; the routine
   *>          only calculates the optimal size of the WORK array, returns
   *>          this value as the first entry of the WORK array, and no error
   *>          message related to LWORK is issued by XERBLA.
   *> \endverbatim
   *>
   *> \param[out] INFO
   *> \verbatim
   *>          INFO is INTEGER
   *>          = 0:  successful exit
   *>          < 0:  if INFO = -i, the i-th argument had an illegal value
   *> \endverbatim
   *
   *  Authors:
   *  ========
   *
   *> \author Univ. of Tennessee
   *> \author Univ. of California Berkeley
   *> \author Univ. of Colorado Denver
   *> \author NAG Ltd.
   *
   *> \date December 2016
   *
   *> \ingroup doubleOTHERcomputational
   *
   *  =====================================================================
       SUBROUTINE DORMQL( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,        SUBROUTINE DORMQL( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
      $                   WORK, LWORK, INFO )       $                   WORK, LWORK, INFO )
 *  *
 *  -- LAPACK routine (version 3.2) --  *  -- LAPACK computational routine (version 3.7.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 *     November 2006  *     December 2016
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          SIDE, TRANS        CHARACTER          SIDE, TRANS
Line 14 Line 180
       DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )        DOUBLE PRECISION   A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
 *     ..  *     ..
 *  *
 *  Purpose  
 *  =======  
 *  
 *  DORMQL overwrites the general real M-by-N matrix C with  
 *  
 *                  SIDE = 'L'     SIDE = 'R'  
 *  TRANS = 'N':      Q * C          C * Q  
 *  TRANS = 'T':      Q**T * C       C * Q**T  
 *  
 *  where Q is a real orthogonal matrix defined as the product of k  
 *  elementary reflectors  
 *  
 *        Q = H(k) . . . H(2) H(1)  
 *  
 *  as returned by DGEQLF. Q is of order M if SIDE = 'L' and of order N  
 *  if SIDE = 'R'.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  SIDE    (input) CHARACTER*1  
 *          = 'L': apply Q or Q**T from the Left;  
 *          = 'R': apply Q or Q**T from the Right.  
 *  
 *  TRANS   (input) CHARACTER*1  
 *          = 'N':  No transpose, apply Q;  
 *          = 'T':  Transpose, apply Q**T.  
 *  
 *  M       (input) INTEGER  
 *          The number of rows of the matrix C. M >= 0.  
 *  
 *  N       (input) INTEGER  
 *          The number of columns of the matrix C. N >= 0.  
 *  
 *  K       (input) INTEGER  
 *          The number of elementary reflectors whose product defines  
 *          the matrix Q.  
 *          If SIDE = 'L', M >= K >= 0;  
 *          if SIDE = 'R', N >= K >= 0.  
 *  
 *  A       (input) DOUBLE PRECISION array, dimension (LDA,K)  
 *          The i-th column must contain the vector which defines the  
 *          elementary reflector H(i), for i = 1,2,...,k, as returned by  
 *          DGEQLF in the last k columns of its array argument A.  
 *          A is modified by the routine but restored on exit.  
 *  
 *  LDA     (input) INTEGER  
 *          The leading dimension of the array A.  
 *          If SIDE = 'L', LDA >= max(1,M);  
 *          if SIDE = 'R', LDA >= max(1,N).  
 *  
 *  TAU     (input) DOUBLE PRECISION array, dimension (K)  
 *          TAU(i) must contain the scalar factor of the elementary  
 *          reflector H(i), as returned by DGEQLF.  
 *  
 *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)  
 *          On entry, the M-by-N matrix C.  
 *          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.  
 *  
 *  LDC     (input) INTEGER  
 *          The leading dimension of the array C. LDC >= max(1,M).  
 *  
 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))  
 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.  
 *  
 *  LWORK   (input) INTEGER  
 *          The dimension of the array WORK.  
 *          If SIDE = 'L', LWORK >= max(1,N);  
 *          if SIDE = 'R', LWORK >= max(1,M).  
 *          For optimum performance LWORK >= N*NB if SIDE = 'L', and  
 *          LWORK >= M*NB if SIDE = 'R', where NB is the optimal  
 *          blocksize.  
 *  
 *          If LWORK = -1, then a workspace query is assumed; the routine  
 *          only calculates the optimal size of the WORK array, returns  
 *          this value as the first entry of the WORK array, and no error  
 *          message related to LWORK is issued by XERBLA.  
 *  
 *  INFO    (output) INTEGER  
 *          = 0:  successful exit  
 *          < 0:  if INFO = -i, the i-th argument had an illegal value  
 *  
 *  =====================================================================  *  =====================================================================
 *  *
 *     .. Parameters ..  *     .. Parameters ..
       INTEGER            NBMAX, LDT        INTEGER            NBMAX, LDT, TSIZE
       PARAMETER          ( NBMAX = 64, LDT = NBMAX+1 )        PARAMETER          ( NBMAX = 64, LDT = NBMAX+1,
        $                     TSIZE = LDT*NBMAX )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       LOGICAL            LEFT, LQUERY, NOTRAN        LOGICAL            LEFT, LQUERY, NOTRAN
       INTEGER            I, I1, I2, I3, IB, IINFO, IWS, LDWORK, LWKOPT,        INTEGER            I, I1, I2, I3, IB, IINFO, IWT, LDWORK, LWKOPT,
      $                   MI, NB, NBMIN, NI, NQ, NW       $                   MI, NB, NBMIN, NI, NQ, NW
 *     ..  *     ..
 *     .. Local Arrays ..  
       DOUBLE PRECISION   T( LDT, NBMAX )  
 *     ..  
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME
       INTEGER            ILAENV        INTEGER            ILAENV
Line 153 Line 235
          INFO = -7           INFO = -7
       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN        ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
          INFO = -10           INFO = -10
         ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
            INFO = -12
       END IF        END IF
 *  *
       IF( INFO.EQ.0 ) THEN        IF( INFO.EQ.0 ) THEN
   *
   *        Compute the workspace requirements
   *
          IF( M.EQ.0 .OR. N.EQ.0 ) THEN           IF( M.EQ.0 .OR. N.EQ.0 ) THEN
             LWKOPT = 1              LWKOPT = 1
          ELSE           ELSE
 *  
 *           Determine the block size.  NB may be at most NBMAX, where  
 *           NBMAX is used to define the local array T.  
 *  
             NB = MIN( NBMAX, ILAENV( 1, 'DORMQL', SIDE // TRANS, M, N,              NB = MIN( NBMAX, ILAENV( 1, 'DORMQL', SIDE // TRANS, M, N,
      $                               K, -1 ) )       $                               K, -1 ) )
             LWKOPT = NW*NB              LWKOPT = NW*NB + TSIZE
          END IF           END IF
          WORK( 1 ) = LWKOPT           WORK( 1 ) = LWKOPT
 *  
          IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN  
             INFO = -12  
          END IF  
       END IF        END IF
 *  *
       IF( INFO.NE.0 ) THEN        IF( INFO.NE.0 ) THEN
Line 190 Line 269
       NBMIN = 2        NBMIN = 2
       LDWORK = NW        LDWORK = NW
       IF( NB.GT.1 .AND. NB.LT.K ) THEN        IF( NB.GT.1 .AND. NB.LT.K ) THEN
          IWS = NW*NB           IF( LWORK.LT.NW*NB+TSIZE ) THEN
          IF( LWORK.LT.IWS ) THEN              NB = (LWORK-TSIZE) / LDWORK
             NB = LWORK / LDWORK  
             NBMIN = MAX( 2, ILAENV( 2, 'DORMQL', SIDE // TRANS, M, N, K,              NBMIN = MAX( 2, ILAENV( 2, 'DORMQL', SIDE // TRANS, M, N, K,
      $              -1 ) )       $              -1 ) )
          END IF           END IF
       ELSE  
          IWS = NW  
       END IF        END IF
 *  *
       IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN        IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN
Line 210 Line 286
 *  *
 *        Use blocked code  *        Use blocked code
 *  *
            IWT = 1 + NW*NB
          IF( ( LEFT .AND. NOTRAN ) .OR.           IF( ( LEFT .AND. NOTRAN ) .OR.
      $       ( .NOT.LEFT .AND. .NOT.NOTRAN ) ) THEN       $       ( .NOT.LEFT .AND. .NOT.NOTRAN ) ) THEN
             I1 = 1              I1 = 1
Line 234 Line 311
 *           H = H(i+ib-1) . . . H(i+1) H(i)  *           H = H(i+ib-1) . . . H(i+1) H(i)
 *  *
             CALL DLARFT( 'Backward', 'Columnwise', NQ-K+I+IB-1, IB,              CALL DLARFT( 'Backward', 'Columnwise', NQ-K+I+IB-1, IB,
      $                   A( 1, I ), LDA, TAU( I ), T, LDT )       $                   A( 1, I ), LDA, TAU( I ), WORK( IWT ), LDT )
             IF( LEFT ) THEN              IF( LEFT ) THEN
 *  *
 *              H or H' is applied to C(1:m-k+i+ib-1,1:n)  *              H or H**T is applied to C(1:m-k+i+ib-1,1:n)
 *  *
                MI = M - K + I + IB - 1                 MI = M - K + I + IB - 1
             ELSE              ELSE
 *  *
 *              H or H' is applied to C(1:m,1:n-k+i+ib-1)  *              H or H**T is applied to C(1:m,1:n-k+i+ib-1)
 *  *
                NI = N - K + I + IB - 1                 NI = N - K + I + IB - 1
             END IF              END IF
 *  *
 *           Apply H or H'  *           Apply H or H**T
 *  *
             CALL DLARFB( SIDE, TRANS, 'Backward', 'Columnwise', MI, NI,              CALL DLARFB( SIDE, TRANS, 'Backward', 'Columnwise', MI, NI,
      $                   IB, A( 1, I ), LDA, T, LDT, C, LDC, WORK,       $                   IB, A( 1, I ), LDA, WORK( IWT ), LDT, C, LDC,
      $                   LDWORK )       $                   WORK, LDWORK )
    10    CONTINUE     10    CONTINUE
       END IF        END IF
       WORK( 1 ) = LWKOPT        WORK( 1 ) = LWKOPT
Removed from v.1.6  
changed lines
  Added in v.1.18

CVSweb interface <joel.bertrand@systella.fr>