Return to dorgqr.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE DORGQR( M, N, K, A, LDA, TAU, WORK, LWORK, INFO )
    2: *
    3: *  -- LAPACK routine (version 3.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     November 2006
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            INFO, K, LDA, LWORK, M, N
   10: *     ..
   11: *     .. Array Arguments ..
   12:       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
   13: *     ..
   14: *
   15: *  Purpose
   16: *  =======
   17: *
   18: *  DORGQR generates an M-by-N real matrix Q with orthonormal columns,
   19: *  which is defined as the first N columns of a product of K elementary
   20: *  reflectors of order M
   21: *
   22: *        Q  =  H(1) H(2) . . . H(k)
   23: *
   24: *  as returned by DGEQRF.
   25: *
   26: *  Arguments
   27: *  =========
   28: *
   29: *  M       (input) INTEGER
   30: *          The number of rows of the matrix Q. M >= 0.
   31: *
   32: *  N       (input) INTEGER
   33: *          The number of columns of the matrix Q. M >= N >= 0.
   34: *
   35: *  K       (input) INTEGER
   36: *          The number of elementary reflectors whose product defines the
   37: *          matrix Q. N >= K >= 0.
   38: *
   39: *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
   40: *          On entry, the i-th column must contain the vector which
   41: *          defines the elementary reflector H(i), for i = 1,2,...,k, as
   42: *          returned by DGEQRF in the first k columns of its array
   43: *          argument A.
   44: *          On exit, the M-by-N matrix Q.
   45: *
   46: *  LDA     (input) INTEGER
   47: *          The first dimension of the array A. LDA >= max(1,M).
   48: *
   49: *  TAU     (input) DOUBLE PRECISION array, dimension (K)
   50: *          TAU(i) must contain the scalar factor of the elementary
   51: *          reflector H(i), as returned by DGEQRF.
   52: *
   53: *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
   54: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
   55: *
   56: *  LWORK   (input) INTEGER
   57: *          The dimension of the array WORK. LWORK >= max(1,N).
   58: *          For optimum performance LWORK >= N*NB, where NB is the
   59: *          optimal blocksize.
   60: *
   61: *          If LWORK = -1, then a workspace query is assumed; the routine
   62: *          only calculates the optimal size of the WORK array, returns
   63: *          this value as the first entry of the WORK array, and no error
   64: *          message related to LWORK is issued by XERBLA.
   65: *
   66: *  INFO    (output) INTEGER
   67: *          = 0:  successful exit
   68: *          < 0:  if INFO = -i, the i-th argument has an illegal value
   69: *
   70: *  =====================================================================
   71: *
   72: *     .. Parameters ..
   73:       DOUBLE PRECISION   ZERO
   74:       PARAMETER          ( ZERO = 0.0D+0 )
   75: *     ..
   76: *     .. Local Scalars ..
   77:       LOGICAL            LQUERY
   78:       INTEGER            I, IB, IINFO, IWS, J, KI, KK, L, LDWORK,
   79:      $                   LWKOPT, NB, NBMIN, NX
   80: *     ..
   81: *     .. External Subroutines ..
   82:       EXTERNAL           DLARFB, DLARFT, DORG2R, XERBLA
   83: *     ..
   84: *     .. Intrinsic Functions ..
   85:       INTRINSIC          MAX, MIN
   86: *     ..
   87: *     .. External Functions ..
   88:       INTEGER            ILAENV
   89:       EXTERNAL           ILAENV
   90: *     ..
   91: *     .. Executable Statements ..
   92: *
   93: *     Test the input arguments
   94: *
   95:       INFO = 0
   96:       NB = ILAENV( 1, 'DORGQR', ' ', M, N, K, -1 )
   97:       LWKOPT = MAX( 1, N )*NB
   98:       WORK( 1 ) = LWKOPT
   99:       LQUERY = ( LWORK.EQ.-1 )
  100:       IF( M.LT.0 ) THEN
  101:          INFO = -1
  102:       ELSE IF( N.LT.0 .OR. N.GT.M ) THEN
  103:          INFO = -2
  104:       ELSE IF( K.LT.0 .OR. K.GT.N ) THEN
  105:          INFO = -3
  106:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  107:          INFO = -5
  108:       ELSE IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
  109:          INFO = -8
  110:       END IF
  111:       IF( INFO.NE.0 ) THEN
  112:          CALL XERBLA( 'DORGQR', -INFO )
  113:          RETURN
  114:       ELSE IF( LQUERY ) THEN
  115:          RETURN
  116:       END IF
  117: *
  118: *     Quick return if possible
  119: *
  120:       IF( N.LE.0 ) THEN
  121:          WORK( 1 ) = 1
  122:          RETURN
  123:       END IF
  124: *
  125:       NBMIN = 2
  126:       NX = 0
  127:       IWS = N
  128:       IF( NB.GT.1 .AND. NB.LT.K ) THEN
  129: *
  130: *        Determine when to cross over from blocked to unblocked code.
  131: *
  132:          NX = MAX( 0, ILAENV( 3, 'DORGQR', ' ', M, N, K, -1 ) )
  133:          IF( NX.LT.K ) THEN
  134: *
  135: *           Determine if workspace is large enough for blocked code.
  136: *
  137:             LDWORK = N
  138:             IWS = LDWORK*NB
  139:             IF( LWORK.LT.IWS ) THEN
  140: *
  141: *              Not enough workspace to use optimal NB:  reduce NB and
  142: *              determine the minimum value of NB.
  143: *
  144:                NB = LWORK / LDWORK
  145:                NBMIN = MAX( 2, ILAENV( 2, 'DORGQR', ' ', M, N, K, -1 ) )
  146:             END IF
  147:          END IF
  148:       END IF
  149: *
  150:       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
  151: *
  152: *        Use blocked code after the last block.
  153: *        The first kk columns are handled by the block method.
  154: *
  155:          KI = ( ( K-NX-1 ) / NB )*NB
  156:          KK = MIN( K, KI+NB )
  157: *
  158: *        Set A(1:kk,kk+1:n) to zero.
  159: *
  160:          DO 20 J = KK + 1, N
  161:             DO 10 I = 1, KK
  162:                A( I, J ) = ZERO
  163:    10       CONTINUE
  164:    20    CONTINUE
  165:       ELSE
  166:          KK = 0
  167:       END IF
  168: *
  169: *     Use unblocked code for the last or only block.
  170: *
  171:       IF( KK.LT.N )
  172:      $   CALL DORG2R( M-KK, N-KK, K-KK, A( KK+1, KK+1 ), LDA,
  173:      $                TAU( KK+1 ), WORK, IINFO )
  174: *
  175:       IF( KK.GT.0 ) THEN
  176: *
  177: *        Use blocked code
  178: *
  179:          DO 50 I = KI + 1, 1, -NB
  180:             IB = MIN( NB, K-I+1 )
  181:             IF( I+IB.LE.N ) THEN
  182: *
  183: *              Form the triangular factor of the block reflector
  184: *              H = H(i) H(i+1) . . . H(i+ib-1)
  185: *
  186:                CALL DLARFT( 'Forward', 'Columnwise', M-I+1, IB,
  187:      $                      A( I, I ), LDA, TAU( I ), WORK, LDWORK )
  188: *
  189: *              Apply H to A(i:m,i+ib:n) from the left
  190: *
  191:                CALL DLARFB( 'Left', 'No transpose', 'Forward',
  192:      $                      'Columnwise', M-I+1, N-I-IB+1, IB,
  193:      $                      A( I, I ), LDA, WORK, LDWORK, A( I, I+IB ),
  194:      $                      LDA, WORK( IB+1 ), LDWORK )
  195:             END IF
  196: *
  197: *           Apply H to rows i:m of current block
  198: *
  199:             CALL DORG2R( M-I+1, IB, IB, A( I, I ), LDA, TAU( I ), WORK,
  200:      $                   IINFO )
  201: *
  202: *           Set rows 1:i-1 of current block to zero
  203: *
  204:             DO 40 J = I, I + IB - 1
  205:                DO 30 L = 1, I - 1
  206:                   A( L, J ) = ZERO
  207:    30          CONTINUE
  208:    40       CONTINUE
  209:    50    CONTINUE
  210:       END IF
  211: *
  212:       WORK( 1 ) = IWS
  213:       RETURN
  214: *
  215: *     End of DORGQR
  216: *
  217:       END