Return to zungql.f CVS log

Up to [local] / rpl / lapack / lapack

Mise à jour globale de Lapack 3.2.2.

    1:       SUBROUTINE ZUNGQL( 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:       COMPLEX*16         A( LDA, * ), TAU( * ), WORK( * )
   13: *     ..
   14: *
   15: *  Purpose
   16: *  =======
   17: *
   18: *  ZUNGQL generates an M-by-N complex matrix Q with orthonormal columns,
   19: *  which is defined as the last N columns of a product of K elementary
   20: *  reflectors of order M
   21: *
   22: *        Q  =  H(k) . . . H(2) H(1)
   23: *
   24: *  as returned by ZGEQLF.
   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) COMPLEX*16 array, dimension (LDA,N)
   40: *          On entry, the (n-k+i)-th column must contain the vector which
   41: *          defines the elementary reflector H(i), for i = 1,2,...,k, as
   42: *          returned by ZGEQLF in the last 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) COMPLEX*16 array, dimension (K)
   50: *          TAU(i) must contain the scalar factor of the elementary
   51: *          reflector H(i), as returned by ZGEQLF.
   52: *
   53: *  WORK    (workspace/output) COMPLEX*16 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:       COMPLEX*16         ZERO
   74:       PARAMETER          ( ZERO = ( 0.0D+0, 0.0D+0 ) )
   75: *     ..
   76: *     .. Local Scalars ..
   77:       LOGICAL            LQUERY
   78:       INTEGER            I, IB, IINFO, IWS, J, KK, L, LDWORK, LWKOPT,
   79:      $                   NB, NBMIN, NX
   80: *     ..
   81: *     .. External Subroutines ..
   82:       EXTERNAL           XERBLA, ZLARFB, ZLARFT, ZUNG2L
   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:       LQUERY = ( LWORK.EQ.-1 )
   97:       IF( M.LT.0 ) THEN
   98:          INFO = -1
   99:       ELSE IF( N.LT.0 .OR. N.GT.M ) THEN
  100:          INFO = -2
  101:       ELSE IF( K.LT.0 .OR. K.GT.N ) THEN
  102:          INFO = -3
  103:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
  104:          INFO = -5
  105:       END IF
  106: *
  107:       IF( INFO.EQ.0 ) THEN
  108:          IF( N.EQ.0 ) THEN
  109:             LWKOPT = 1
  110:          ELSE
  111:             NB = ILAENV( 1, 'ZUNGQL', ' ', M, N, K, -1 )
  112:             LWKOPT = N*NB
  113:          END IF
  114:          WORK( 1 ) = LWKOPT
  115: *
  116:          IF( LWORK.LT.MAX( 1, N ) .AND. .NOT.LQUERY ) THEN
  117:             INFO = -8
  118:          END IF
  119:       END IF
  120: *
  121:       IF( INFO.NE.0 ) THEN
  122:          CALL XERBLA( 'ZUNGQL', -INFO )
  123:          RETURN
  124:       ELSE IF( LQUERY ) THEN
  125:          RETURN
  126:       END IF
  127: *
  128: *     Quick return if possible
  129: *
  130:       IF( N.LE.0 ) THEN
  131:          RETURN
  132:       END IF
  133: *
  134:       NBMIN = 2
  135:       NX = 0
  136:       IWS = N
  137:       IF( NB.GT.1 .AND. NB.LT.K ) THEN
  138: *
  139: *        Determine when to cross over from blocked to unblocked code.
  140: *
  141:          NX = MAX( 0, ILAENV( 3, 'ZUNGQL', ' ', M, N, K, -1 ) )
  142:          IF( NX.LT.K ) THEN
  143: *
  144: *           Determine if workspace is large enough for blocked code.
  145: *
  146:             LDWORK = N
  147:             IWS = LDWORK*NB
  148:             IF( LWORK.LT.IWS ) THEN
  149: *
  150: *              Not enough workspace to use optimal NB:  reduce NB and
  151: *              determine the minimum value of NB.
  152: *
  153:                NB = LWORK / LDWORK
  154:                NBMIN = MAX( 2, ILAENV( 2, 'ZUNGQL', ' ', M, N, K, -1 ) )
  155:             END IF
  156:          END IF
  157:       END IF
  158: *
  159:       IF( NB.GE.NBMIN .AND. NB.LT.K .AND. NX.LT.K ) THEN
  160: *
  161: *        Use blocked code after the first block.
  162: *        The last kk columns are handled by the block method.
  163: *
  164:          KK = MIN( K, ( ( K-NX+NB-1 ) / NB )*NB )
  165: *
  166: *        Set A(m-kk+1:m,1:n-kk) to zero.
  167: *
  168:          DO 20 J = 1, N - KK
  169:             DO 10 I = M - KK + 1, M
  170:                A( I, J ) = ZERO
  171:    10       CONTINUE
  172:    20    CONTINUE
  173:       ELSE
  174:          KK = 0
  175:       END IF
  176: *
  177: *     Use unblocked code for the first or only block.
  178: *
  179:       CALL ZUNG2L( M-KK, N-KK, K-KK, A, LDA, TAU, WORK, IINFO )
  180: *
  181:       IF( KK.GT.0 ) THEN
  182: *
  183: *        Use blocked code
  184: *
  185:          DO 50 I = K - KK + 1, K, NB
  186:             IB = MIN( NB, K-I+1 )
  187:             IF( N-K+I.GT.1 ) THEN
  188: *
  189: *              Form the triangular factor of the block reflector
  190: *              H = H(i+ib-1) . . . H(i+1) H(i)
  191: *
  192:                CALL ZLARFT( 'Backward', 'Columnwise', M-K+I+IB-1, IB,
  193:      $                      A( 1, N-K+I ), LDA, TAU( I ), WORK, LDWORK )
  194: *
  195: *              Apply H to A(1:m-k+i+ib-1,1:n-k+i-1) from the left
  196: *
  197:                CALL ZLARFB( 'Left', 'No transpose', 'Backward',
  198:      $                      'Columnwise', M-K+I+IB-1, N-K+I-1, IB,
  199:      $                      A( 1, N-K+I ), LDA, WORK, LDWORK, A, LDA,
  200:      $                      WORK( IB+1 ), LDWORK )
  201:             END IF
  202: *
  203: *           Apply H to rows 1:m-k+i+ib-1 of current block
  204: *
  205:             CALL ZUNG2L( M-K+I+IB-1, IB, IB, A( 1, N-K+I ), LDA,
  206:      $                   TAU( I ), WORK, IINFO )
  207: *
  208: *           Set rows m-k+i+ib:m of current block to zero
  209: *
  210:             DO 40 J = N - K + I, N - K + I + IB - 1
  211:                DO 30 L = M - K + I + IB, M
  212:                   A( L, J ) = ZERO
  213:    30          CONTINUE
  214:    40       CONTINUE
  215:    50    CONTINUE
  216:       END IF
  217: *
  218:       WORK( 1 ) = IWS
  219:       RETURN
  220: *
  221: *     End of ZUNGQL
  222: *
  223:       END