Annotation of rpl/lapack/lapack/zunmqr.f, revision 1.4

1.1       bertrand    1:       SUBROUTINE ZUNMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
                      2:      $                   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
                     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: *  ZUNMQR overwrites the general complex M-by-N matrix C with
                     21: *
                     22: *                  SIDE = 'L'     SIDE = 'R'
                     23: *  TRANS = 'N':      Q * C          C * Q
                     24: *  TRANS = 'C':      Q**H * C       C * Q**H
                     25: *
                     26: *  where Q is a complex unitary matrix defined as the product of k
                     27: *  elementary reflectors
                     28: *
                     29: *        Q = H(1) H(2) . . . H(k)
                     30: *
                     31: *  as returned by ZGEQRF. Q is of order M if SIDE = 'L' and of order N
                     32: *  if SIDE = 'R'.
                     33: *
                     34: *  Arguments
                     35: *  =========
                     36: *
                     37: *  SIDE    (input) CHARACTER*1
                     38: *          = 'L': apply Q or Q**H from the Left;
                     39: *          = 'R': apply Q or Q**H from the Right.
                     40: *
                     41: *  TRANS   (input) CHARACTER*1
                     42: *          = 'N':  No transpose, apply Q;
                     43: *          = 'C':  Conjugate transpose, apply Q**H.
                     44: *
                     45: *  M       (input) INTEGER
                     46: *          The number of rows of the matrix C. M >= 0.
                     47: *
                     48: *  N       (input) INTEGER
                     49: *          The number of columns of the matrix C. N >= 0.
                     50: *
                     51: *  K       (input) INTEGER
                     52: *          The number of elementary reflectors whose product defines
                     53: *          the matrix Q.
                     54: *          If SIDE = 'L', M >= K >= 0;
                     55: *          if SIDE = 'R', N >= K >= 0.
                     56: *
                     57: *  A       (input) COMPLEX*16 array, dimension (LDA,K)
                     58: *          The i-th column must contain the vector which defines the
                     59: *          elementary reflector H(i), for i = 1,2,...,k, as returned by
                     60: *          ZGEQRF in the first k columns of its array argument A.
                     61: *          A is modified by the routine but restored on exit.
                     62: *
                     63: *  LDA     (input) INTEGER
                     64: *          The leading dimension of the array A.
                     65: *          If SIDE = 'L', LDA >= max(1,M);
                     66: *          if SIDE = 'R', LDA >= max(1,N).
                     67: *
                     68: *  TAU     (input) COMPLEX*16 array, dimension (K)
                     69: *          TAU(i) must contain the scalar factor of the elementary
                     70: *          reflector H(i), as returned by ZGEQRF.
                     71: *
                     72: *  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
                     73: *          On entry, the M-by-N matrix C.
                     74: *          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
                     75: *
                     76: *  LDC     (input) INTEGER
                     77: *          The leading dimension of the array C. LDC >= max(1,M).
                     78: *
                     79: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
                     80: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
                     81: *
                     82: *  LWORK   (input) INTEGER
                     83: *          The dimension of the array WORK.
                     84: *          If SIDE = 'L', LWORK >= max(1,N);
                     85: *          if SIDE = 'R', LWORK >= max(1,M).
                     86: *          For optimum performance LWORK >= N*NB if SIDE = 'L', and
                     87: *          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
                     88: *          blocksize.
                     89: *
                     90: *          If LWORK = -1, then a workspace query is assumed; the routine
                     91: *          only calculates the optimal size of the WORK array, returns
                     92: *          this value as the first entry of the WORK array, and no error
                     93: *          message related to LWORK is issued by XERBLA.
                     94: *
                     95: *  INFO    (output) INTEGER
                     96: *          = 0:  successful exit
                     97: *          < 0:  if INFO = -i, the i-th argument had an illegal value
                     98: *
                     99: *  =====================================================================
                    100: *
                    101: *     .. Parameters ..
                    102:       INTEGER            NBMAX, LDT
                    103:       PARAMETER          ( NBMAX = 64, LDT = NBMAX+1 )
                    104: *     ..
                    105: *     .. Local Scalars ..
                    106:       LOGICAL            LEFT, LQUERY, NOTRAN
                    107:       INTEGER            I, I1, I2, I3, IB, IC, IINFO, IWS, JC, LDWORK,
                    108:      $                   LWKOPT, MI, NB, NBMIN, NI, NQ, NW
                    109: *     ..
                    110: *     .. Local Arrays ..
                    111:       COMPLEX*16         T( LDT, NBMAX )
                    112: *     ..
                    113: *     .. External Functions ..
                    114:       LOGICAL            LSAME
                    115:       INTEGER            ILAENV
                    116:       EXTERNAL           LSAME, ILAENV
                    117: *     ..
                    118: *     .. External Subroutines ..
                    119:       EXTERNAL           XERBLA, ZLARFB, ZLARFT, ZUNM2R
                    120: *     ..
                    121: *     .. Intrinsic Functions ..
                    122:       INTRINSIC          MAX, MIN
                    123: *     ..
                    124: *     .. Executable Statements ..
                    125: *
                    126: *     Test the input arguments
                    127: *
                    128:       INFO = 0
                    129:       LEFT = LSAME( SIDE, 'L' )
                    130:       NOTRAN = LSAME( TRANS, 'N' )
                    131:       LQUERY = ( LWORK.EQ.-1 )
                    132: *
                    133: *     NQ is the order of Q and NW is the minimum dimension of WORK
                    134: *
                    135:       IF( LEFT ) THEN
                    136:          NQ = M
                    137:          NW = N
                    138:       ELSE
                    139:          NQ = N
                    140:          NW = M
                    141:       END IF
                    142:       IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
                    143:          INFO = -1
                    144:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
                    145:          INFO = -2
                    146:       ELSE IF( M.LT.0 ) THEN
                    147:          INFO = -3
                    148:       ELSE IF( N.LT.0 ) THEN
                    149:          INFO = -4
                    150:       ELSE IF( K.LT.0 .OR. K.GT.NQ ) THEN
                    151:          INFO = -5
                    152:       ELSE IF( LDA.LT.MAX( 1, NQ ) ) THEN
                    153:          INFO = -7
                    154:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
                    155:          INFO = -10
                    156:       ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
                    157:          INFO = -12
                    158:       END IF
                    159: *
                    160:       IF( INFO.EQ.0 ) THEN
                    161: *
                    162: *        Determine the block size.  NB may be at most NBMAX, where NBMAX
                    163: *        is used to define the local array T.
                    164: *
                    165:          NB = MIN( NBMAX, ILAENV( 1, 'ZUNMQR', SIDE // TRANS, M, N, K,
                    166:      $        -1 ) )
                    167:          LWKOPT = MAX( 1, NW )*NB
                    168:          WORK( 1 ) = LWKOPT
                    169:       END IF
                    170: *
                    171:       IF( INFO.NE.0 ) THEN
                    172:          CALL XERBLA( 'ZUNMQR', -INFO )
                    173:          RETURN
                    174:       ELSE IF( LQUERY ) THEN
                    175:          RETURN
                    176:       END IF
                    177: *
                    178: *     Quick return if possible
                    179: *
                    180:       IF( M.EQ.0 .OR. N.EQ.0 .OR. K.EQ.0 ) THEN
                    181:          WORK( 1 ) = 1
                    182:          RETURN
                    183:       END IF
                    184: *
                    185:       NBMIN = 2
                    186:       LDWORK = NW
                    187:       IF( NB.GT.1 .AND. NB.LT.K ) THEN
                    188:          IWS = NW*NB
                    189:          IF( LWORK.LT.IWS ) THEN
                    190:             NB = LWORK / LDWORK
                    191:             NBMIN = MAX( 2, ILAENV( 2, 'ZUNMQR', SIDE // TRANS, M, N, K,
                    192:      $              -1 ) )
                    193:          END IF
                    194:       ELSE
                    195:          IWS = NW
                    196:       END IF
                    197: *
                    198:       IF( NB.LT.NBMIN .OR. NB.GE.K ) THEN
                    199: *
                    200: *        Use unblocked code
                    201: *
                    202:          CALL ZUNM2R( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK,
                    203:      $                IINFO )
                    204:       ELSE
                    205: *
                    206: *        Use blocked code
                    207: *
                    208:          IF( ( LEFT .AND. .NOT.NOTRAN ) .OR.
                    209:      $       ( .NOT.LEFT .AND. NOTRAN ) ) THEN
                    210:             I1 = 1
                    211:             I2 = K
                    212:             I3 = NB
                    213:          ELSE
                    214:             I1 = ( ( K-1 ) / NB )*NB + 1
                    215:             I2 = 1
                    216:             I3 = -NB
                    217:          END IF
                    218: *
                    219:          IF( LEFT ) THEN
                    220:             NI = N
                    221:             JC = 1
                    222:          ELSE
                    223:             MI = M
                    224:             IC = 1
                    225:          END IF
                    226: *
                    227:          DO 10 I = I1, I2, I3
                    228:             IB = MIN( NB, K-I+1 )
                    229: *
                    230: *           Form the triangular factor of the block reflector
                    231: *           H = H(i) H(i+1) . . . H(i+ib-1)
                    232: *
                    233:             CALL ZLARFT( 'Forward', 'Columnwise', NQ-I+1, IB, A( I, I ),
                    234:      $                   LDA, TAU( I ), T, LDT )
                    235:             IF( LEFT ) THEN
                    236: *
                    237: *              H or H' is applied to C(i:m,1:n)
                    238: *
                    239:                MI = M - I + 1
                    240:                IC = I
                    241:             ELSE
                    242: *
                    243: *              H or H' is applied to C(1:m,i:n)
                    244: *
                    245:                NI = N - I + 1
                    246:                JC = I
                    247:             END IF
                    248: *
                    249: *           Apply H or H'
                    250: *
                    251:             CALL ZLARFB( SIDE, TRANS, 'Forward', 'Columnwise', MI, NI,
                    252:      $                   IB, A( I, I ), LDA, T, LDT, C( IC, JC ), LDC,
                    253:      $                   WORK, LDWORK )
                    254:    10    CONTINUE
                    255:       END IF
                    256:       WORK( 1 ) = LWKOPT
                    257:       RETURN
                    258: *
                    259: *     End of ZUNMQR
                    260: *
                    261:       END

CVSweb interface <joel.bertrand@systella.fr>