Return to zunmbr.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence

    1:       SUBROUTINE ZUNMBR( VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C,
    2:      $                   LDC, WORK, LWORK, INFO )
    3: *
    4: *  -- LAPACK routine (version 3.2) --
    5: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    6: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    7: *     November 2006
    8: *
    9: *     .. Scalar Arguments ..
   10:       CHARACTER          SIDE, TRANS, VECT
   11:       INTEGER            INFO, K, LDA, LDC, LWORK, M, N
   12: *     ..
   13: *     .. Array Arguments ..
   14:       COMPLEX*16         A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
   15: *     ..
   16: *
   17: *  Purpose
   18: *  =======
   19: *
   20: *  If VECT = 'Q', ZUNMBR overwrites the general complex M-by-N matrix C
   21: *  with
   22: *                  SIDE = 'L'     SIDE = 'R'
   23: *  TRANS = 'N':      Q * C          C * Q
   24: *  TRANS = 'C':      Q**H * C       C * Q**H
   25: *
   26: *  If VECT = 'P', ZUNMBR overwrites the general complex M-by-N matrix C
   27: *  with
   28: *                  SIDE = 'L'     SIDE = 'R'
   29: *  TRANS = 'N':      P * C          C * P
   30: *  TRANS = 'C':      P**H * C       C * P**H
   31: *
   32: *  Here Q and P**H are the unitary matrices determined by ZGEBRD when
   33: *  reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q
   34: *  and P**H are defined as products of elementary reflectors H(i) and
   35: *  G(i) respectively.
   36: *
   37: *  Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the
   38: *  order of the unitary matrix Q or P**H that is applied.
   39: *
   40: *  If VECT = 'Q', A is assumed to have been an NQ-by-K matrix:
   41: *  if nq >= k, Q = H(1) H(2) . . . H(k);
   42: *  if nq < k, Q = H(1) H(2) . . . H(nq-1).
   43: *
   44: *  If VECT = 'P', A is assumed to have been a K-by-NQ matrix:
   45: *  if k < nq, P = G(1) G(2) . . . G(k);
   46: *  if k >= nq, P = G(1) G(2) . . . G(nq-1).
   47: *
   48: *  Arguments
   49: *  =========
   50: *
   51: *  VECT    (input) CHARACTER*1
   52: *          = 'Q': apply Q or Q**H;
   53: *          = 'P': apply P or P**H.
   54: *
   55: *  SIDE    (input) CHARACTER*1
   56: *          = 'L': apply Q, Q**H, P or P**H from the Left;
   57: *          = 'R': apply Q, Q**H, P or P**H from the Right.
   58: *
   59: *  TRANS   (input) CHARACTER*1
   60: *          = 'N':  No transpose, apply Q or P;
   61: *          = 'C':  Conjugate transpose, apply Q**H or P**H.
   62: *
   63: *  M       (input) INTEGER
   64: *          The number of rows of the matrix C. M >= 0.
   65: *
   66: *  N       (input) INTEGER
   67: *          The number of columns of the matrix C. N >= 0.
   68: *
   69: *  K       (input) INTEGER
   70: *          If VECT = 'Q', the number of columns in the original
   71: *          matrix reduced by ZGEBRD.
   72: *          If VECT = 'P', the number of rows in the original
   73: *          matrix reduced by ZGEBRD.
   74: *          K >= 0.
   75: *
   76: *  A       (input) COMPLEX*16 array, dimension
   77: *                                (LDA,min(nq,K)) if VECT = 'Q'
   78: *                                (LDA,nq)        if VECT = 'P'
   79: *          The vectors which define the elementary reflectors H(i) and
   80: *          G(i), whose products determine the matrices Q and P, as
   81: *          returned by ZGEBRD.
   82: *
   83: *  LDA     (input) INTEGER
   84: *          The leading dimension of the array A.
   85: *          If VECT = 'Q', LDA >= max(1,nq);
   86: *          if VECT = 'P', LDA >= max(1,min(nq,K)).
   87: *
   88: *  TAU     (input) COMPLEX*16 array, dimension (min(nq,K))
   89: *          TAU(i) must contain the scalar factor of the elementary
   90: *          reflector H(i) or G(i) which determines Q or P, as returned
   91: *          by ZGEBRD in the array argument TAUQ or TAUP.
   92: *
   93: *  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
   94: *          On entry, the M-by-N matrix C.
   95: *          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q
   96: *          or P*C or P**H*C or C*P or C*P**H.
   97: *
   98: *  LDC     (input) INTEGER
   99: *          The leading dimension of the array C. LDC >= max(1,M).
  100: *
  101: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
  102: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  103: *
  104: *  LWORK   (input) INTEGER
  105: *          The dimension of the array WORK.
  106: *          If SIDE = 'L', LWORK >= max(1,N);
  107: *          if SIDE = 'R', LWORK >= max(1,M);
  108: *          if N = 0 or M = 0, LWORK >= 1.
  109: *          For optimum performance LWORK >= max(1,N*NB) if SIDE = 'L',
  110: *          and LWORK >= max(1,M*NB) if SIDE = 'R', where NB is the
  111: *          optimal blocksize. (NB = 0 if M = 0 or N = 0.)
  112: *
  113: *          If LWORK = -1, then a workspace query is assumed; the routine
  114: *          only calculates the optimal size of the WORK array, returns
  115: *          this value as the first entry of the WORK array, and no error
  116: *          message related to LWORK is issued by XERBLA.
  117: *
  118: *  INFO    (output) INTEGER
  119: *          = 0:  successful exit
  120: *          < 0:  if INFO = -i, the i-th argument had an illegal value
  121: *
  122: *  =====================================================================
  123: *
  124: *     .. Local Scalars ..
  125:       LOGICAL            APPLYQ, LEFT, LQUERY, NOTRAN
  126:       CHARACTER          TRANST
  127:       INTEGER            I1, I2, IINFO, LWKOPT, MI, NB, NI, NQ, NW
  128: *     ..
  129: *     .. External Functions ..
  130:       LOGICAL            LSAME
  131:       INTEGER            ILAENV
  132:       EXTERNAL           LSAME, ILAENV
  133: *     ..
  134: *     .. External Subroutines ..
  135:       EXTERNAL           XERBLA, ZUNMLQ, ZUNMQR
  136: *     ..
  137: *     .. Intrinsic Functions ..
  138:       INTRINSIC          MAX, MIN
  139: *     ..
  140: *     .. Executable Statements ..
  141: *
  142: *     Test the input arguments
  143: *
  144:       INFO = 0
  145:       APPLYQ = LSAME( VECT, 'Q' )
  146:       LEFT = LSAME( SIDE, 'L' )
  147:       NOTRAN = LSAME( TRANS, 'N' )
  148:       LQUERY = ( LWORK.EQ.-1 )
  149: *
  150: *     NQ is the order of Q or P and NW is the minimum dimension of WORK
  151: *
  152:       IF( LEFT ) THEN
  153:          NQ = M
  154:          NW = N
  155:       ELSE
  156:          NQ = N
  157:          NW = M
  158:       END IF
  159:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
  160:          NW = 0
  161:       END IF
  162:       IF( .NOT.APPLYQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN
  163:          INFO = -1
  164:       ELSE IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
  165:          INFO = -2
  166:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
  167:          INFO = -3
  168:       ELSE IF( M.LT.0 ) THEN
  169:          INFO = -4
  170:       ELSE IF( N.LT.0 ) THEN
  171:          INFO = -5
  172:       ELSE IF( K.LT.0 ) THEN
  173:          INFO = -6
  174:       ELSE IF( ( APPLYQ .AND. LDA.LT.MAX( 1, NQ ) ) .OR.
  175:      $         ( .NOT.APPLYQ .AND. LDA.LT.MAX( 1, MIN( NQ, K ) ) ) )
  176:      $          THEN
  177:          INFO = -8
  178:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
  179:          INFO = -11
  180:       ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
  181:          INFO = -13
  182:       END IF
  183: *
  184:       IF( INFO.EQ.0 ) THEN
  185:          IF( NW.GT.0 ) THEN
  186:             IF( APPLYQ ) THEN
  187:                IF( LEFT ) THEN
  188:                   NB = ILAENV( 1, 'ZUNMQR', SIDE // TRANS, M-1, N, M-1,
  189:      $                 -1 )
  190:                ELSE
  191:                   NB = ILAENV( 1, 'ZUNMQR', SIDE // TRANS, M, N-1, N-1,
  192:      $                 -1 )
  193:                END IF
  194:             ELSE
  195:                IF( LEFT ) THEN
  196:                   NB = ILAENV( 1, 'ZUNMLQ', SIDE // TRANS, M-1, N, M-1,
  197:      $                 -1 )
  198:                ELSE
  199:                   NB = ILAENV( 1, 'ZUNMLQ', SIDE // TRANS, M, N-1, N-1,
  200:      $                 -1 )
  201:                END IF
  202:             END IF
  203:             LWKOPT = MAX( 1, NW*NB )
  204:          ELSE
  205:             LWKOPT = 1
  206:          END IF
  207:          WORK( 1 ) = LWKOPT
  208:       END IF
  209: *
  210:       IF( INFO.NE.0 ) THEN
  211:          CALL XERBLA( 'ZUNMBR', -INFO )
  212:          RETURN
  213:       ELSE IF( LQUERY ) THEN
  214:          RETURN
  215:       END IF
  216: *
  217: *     Quick return if possible
  218: *
  219:       IF( M.EQ.0 .OR. N.EQ.0 )
  220:      $   RETURN
  221: *
  222:       IF( APPLYQ ) THEN
  223: *
  224: *        Apply Q
  225: *
  226:          IF( NQ.GE.K ) THEN
  227: *
  228: *           Q was determined by a call to ZGEBRD with nq >= k
  229: *
  230:             CALL ZUNMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
  231:      $                   WORK, LWORK, IINFO )
  232:          ELSE IF( NQ.GT.1 ) THEN
  233: *
  234: *           Q was determined by a call to ZGEBRD with nq < k
  235: *
  236:             IF( LEFT ) THEN
  237:                MI = M - 1
  238:                NI = N
  239:                I1 = 2
  240:                I2 = 1
  241:             ELSE
  242:                MI = M
  243:                NI = N - 1
  244:                I1 = 1
  245:                I2 = 2
  246:             END IF
  247:             CALL ZUNMQR( SIDE, TRANS, MI, NI, NQ-1, A( 2, 1 ), LDA, TAU,
  248:      $                   C( I1, I2 ), LDC, WORK, LWORK, IINFO )
  249:          END IF
  250:       ELSE
  251: *
  252: *        Apply P
  253: *
  254:          IF( NOTRAN ) THEN
  255:             TRANST = 'C'
  256:          ELSE
  257:             TRANST = 'N'
  258:          END IF
  259:          IF( NQ.GT.K ) THEN
  260: *
  261: *           P was determined by a call to ZGEBRD with nq > k
  262: *
  263:             CALL ZUNMLQ( SIDE, TRANST, M, N, K, A, LDA, TAU, C, LDC,
  264:      $                   WORK, LWORK, IINFO )
  265:          ELSE IF( NQ.GT.1 ) THEN
  266: *
  267: *           P was determined by a call to ZGEBRD with nq <= k
  268: *
  269:             IF( LEFT ) THEN
  270:                MI = M - 1
  271:                NI = N
  272:                I1 = 2
  273:                I2 = 1
  274:             ELSE
  275:                MI = M
  276:                NI = N - 1
  277:                I1 = 1
  278:                I2 = 2
  279:             END IF
  280:             CALL ZUNMLQ( SIDE, TRANST, MI, NI, NQ-1, A( 1, 2 ), LDA,
  281:      $                   TAU, C( I1, I2 ), LDC, WORK, LWORK, IINFO )
  282:          END IF
  283:       END IF
  284:       WORK( 1 ) = LWKOPT
  285:       RETURN
  286: *
  287: *     End of ZUNMBR
  288: *
  289:       END