001:       SUBROUTINE ZUNMBR( VECT, SIDE, TRANS, M, N, K, A, LDA, TAU, C,
002:      $                   LDC, WORK, LWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
006: *     November 2006
007: *
008: *     .. Scalar Arguments ..
009:       CHARACTER          SIDE, TRANS, VECT
010:       INTEGER            INFO, K, LDA, LDC, LWORK, M, N
011: *     ..
012: *     .. Array Arguments ..
013:       COMPLEX*16         A( LDA, * ), C( LDC, * ), TAU( * ), WORK( * )
014: *     ..
015: *
016: *  Purpose
017: *  =======
018: *
019: *  If VECT = 'Q', ZUNMBR overwrites the general complex M-by-N matrix C
020: *  with
021: *                  SIDE = 'L'     SIDE = 'R'
022: *  TRANS = 'N':      Q * C          C * Q
023: *  TRANS = 'C':      Q**H * C       C * Q**H
024: *
025: *  If VECT = 'P', ZUNMBR overwrites the general complex M-by-N matrix C
026: *  with
027: *                  SIDE = 'L'     SIDE = 'R'
028: *  TRANS = 'N':      P * C          C * P
029: *  TRANS = 'C':      P**H * C       C * P**H
030: *
031: *  Here Q and P**H are the unitary matrices determined by ZGEBRD when
032: *  reducing a complex matrix A to bidiagonal form: A = Q * B * P**H. Q
033: *  and P**H are defined as products of elementary reflectors H(i) and
034: *  G(i) respectively.
035: *
036: *  Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the
037: *  order of the unitary matrix Q or P**H that is applied.
038: *
039: *  If VECT = 'Q', A is assumed to have been an NQ-by-K matrix:
040: *  if nq >= k, Q = H(1) H(2) . . . H(k);
041: *  if nq < k, Q = H(1) H(2) . . . H(nq-1).
042: *
043: *  If VECT = 'P', A is assumed to have been a K-by-NQ matrix:
044: *  if k < nq, P = G(1) G(2) . . . G(k);
045: *  if k >= nq, P = G(1) G(2) . . . G(nq-1).
046: *
047: *  Arguments
048: *  =========
049: *
050: *  VECT    (input) CHARACTER*1
051: *          = 'Q': apply Q or Q**H;
052: *          = 'P': apply P or P**H.
053: *
054: *  SIDE    (input) CHARACTER*1
055: *          = 'L': apply Q, Q**H, P or P**H from the Left;
056: *          = 'R': apply Q, Q**H, P or P**H from the Right.
057: *
058: *  TRANS   (input) CHARACTER*1
059: *          = 'N':  No transpose, apply Q or P;
060: *          = 'C':  Conjugate transpose, apply Q**H or P**H.
061: *
062: *  M       (input) INTEGER
063: *          The number of rows of the matrix C. M >= 0.
064: *
065: *  N       (input) INTEGER
066: *          The number of columns of the matrix C. N >= 0.
067: *
068: *  K       (input) INTEGER
069: *          If VECT = 'Q', the number of columns in the original
070: *          matrix reduced by ZGEBRD.
071: *          If VECT = 'P', the number of rows in the original
072: *          matrix reduced by ZGEBRD.
073: *          K >= 0.
074: *
075: *  A       (input) COMPLEX*16 array, dimension
076: *                                (LDA,min(nq,K)) if VECT = 'Q'
077: *                                (LDA,nq)        if VECT = 'P'
078: *          The vectors which define the elementary reflectors H(i) and
079: *          G(i), whose products determine the matrices Q and P, as
080: *          returned by ZGEBRD.
081: *
082: *  LDA     (input) INTEGER
083: *          The leading dimension of the array A.
084: *          If VECT = 'Q', LDA >= max(1,nq);
085: *          if VECT = 'P', LDA >= max(1,min(nq,K)).
086: *
087: *  TAU     (input) COMPLEX*16 array, dimension (min(nq,K))
088: *          TAU(i) must contain the scalar factor of the elementary
089: *          reflector H(i) or G(i) which determines Q or P, as returned
090: *          by ZGEBRD in the array argument TAUQ or TAUP.
091: *
092: *  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
093: *          On entry, the M-by-N matrix C.
094: *          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q
095: *          or P*C or P**H*C or C*P or C*P**H.
096: *
097: *  LDC     (input) INTEGER
098: *          The leading dimension of the array C. LDC >= max(1,M).
099: *
100: *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
101: *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
102: *
103: *  LWORK   (input) INTEGER
104: *          The dimension of the array WORK.
105: *          If SIDE = 'L', LWORK >= max(1,N);
106: *          if SIDE = 'R', LWORK >= max(1,M);
107: *          if N = 0 or M = 0, LWORK >= 1.
108: *          For optimum performance LWORK >= max(1,N*NB) if SIDE = 'L',
109: *          and LWORK >= max(1,M*NB) if SIDE = 'R', where NB is the
110: *          optimal blocksize. (NB = 0 if M = 0 or N = 0.)
111: *
112: *          If LWORK = -1, then a workspace query is assumed; the routine
113: *          only calculates the optimal size of the WORK array, returns
114: *          this value as the first entry of the WORK array, and no error
115: *          message related to LWORK is issued by XERBLA.
116: *
117: *  INFO    (output) INTEGER
118: *          = 0:  successful exit
119: *          < 0:  if INFO = -i, the i-th argument had an illegal value
120: *
121: *  =====================================================================
122: *
123: *     .. Local Scalars ..
124:       LOGICAL            APPLYQ, LEFT, LQUERY, NOTRAN
125:       CHARACTER          TRANST
126:       INTEGER            I1, I2, IINFO, LWKOPT, MI, NB, NI, NQ, NW
127: *     ..
128: *     .. External Functions ..
129:       LOGICAL            LSAME
130:       INTEGER            ILAENV
131:       EXTERNAL           LSAME, ILAENV
132: *     ..
133: *     .. External Subroutines ..
134:       EXTERNAL           XERBLA, ZUNMLQ, ZUNMQR
135: *     ..
136: *     .. Intrinsic Functions ..
137:       INTRINSIC          MAX, MIN
138: *     ..
139: *     .. Executable Statements ..
140: *
141: *     Test the input arguments
142: *
143:       INFO = 0
144:       APPLYQ = LSAME( VECT, 'Q' )
145:       LEFT = LSAME( SIDE, 'L' )
146:       NOTRAN = LSAME( TRANS, 'N' )
147:       LQUERY = ( LWORK.EQ.-1 )
148: *
149: *     NQ is the order of Q or P and NW is the minimum dimension of WORK
150: *
151:       IF( LEFT ) THEN
152:          NQ = M
153:          NW = N
154:       ELSE
155:          NQ = N
156:          NW = M
157:       END IF
158:       IF( M.EQ.0 .OR. N.EQ.0 ) THEN
159:          NW = 0
160:       END IF
161:       IF( .NOT.APPLYQ .AND. .NOT.LSAME( VECT, 'P' ) ) THEN
162:          INFO = -1
163:       ELSE IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
164:          INFO = -2
165:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
166:          INFO = -3
167:       ELSE IF( M.LT.0 ) THEN
168:          INFO = -4
169:       ELSE IF( N.LT.0 ) THEN
170:          INFO = -5
171:       ELSE IF( K.LT.0 ) THEN
172:          INFO = -6
173:       ELSE IF( ( APPLYQ .AND. LDA.LT.MAX( 1, NQ ) ) .OR.
174:      $         ( .NOT.APPLYQ .AND. LDA.LT.MAX( 1, MIN( NQ, K ) ) ) )
175:      $          THEN
176:          INFO = -8
177:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
178:          INFO = -11
179:       ELSE IF( LWORK.LT.MAX( 1, NW ) .AND. .NOT.LQUERY ) THEN
180:          INFO = -13
181:       END IF
182: *
183:       IF( INFO.EQ.0 ) THEN
184:          IF( NW.GT.0 ) THEN
185:             IF( APPLYQ ) THEN
186:                IF( LEFT ) THEN
187:                   NB = ILAENV( 1, 'ZUNMQR', SIDE // TRANS, M-1, N, M-1,
188:      $                 -1 )
189:                ELSE
190:                   NB = ILAENV( 1, 'ZUNMQR', SIDE // TRANS, M, N-1, N-1,
191:      $                 -1 )
192:                END IF
193:             ELSE
194:                IF( LEFT ) THEN
195:                   NB = ILAENV( 1, 'ZUNMLQ', SIDE // TRANS, M-1, N, M-1,
196:      $                 -1 )
197:                ELSE
198:                   NB = ILAENV( 1, 'ZUNMLQ', SIDE // TRANS, M, N-1, N-1,
199:      $                 -1 )
200:                END IF
201:             END IF
202:             LWKOPT = MAX( 1, NW*NB )
203:          ELSE
204:             LWKOPT = 1
205:          END IF
206:          WORK( 1 ) = LWKOPT
207:       END IF
208: *
209:       IF( INFO.NE.0 ) THEN
210:          CALL XERBLA( 'ZUNMBR', -INFO )
211:          RETURN
212:       ELSE IF( LQUERY ) THEN
213:          RETURN
214:       END IF
215: *
216: *     Quick return if possible
217: *
218:       IF( M.EQ.0 .OR. N.EQ.0 )
219:      $   RETURN
220: *
221:       IF( APPLYQ ) THEN
222: *
223: *        Apply Q
224: *
225:          IF( NQ.GE.K ) THEN
226: *
227: *           Q was determined by a call to ZGEBRD with nq >= k
228: *
229:             CALL ZUNMQR( SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC,
230:      $                   WORK, LWORK, IINFO )
231:          ELSE IF( NQ.GT.1 ) THEN
232: *
233: *           Q was determined by a call to ZGEBRD with nq < k
234: *
235:             IF( LEFT ) THEN
236:                MI = M - 1
237:                NI = N
238:                I1 = 2
239:                I2 = 1
240:             ELSE
241:                MI = M
242:                NI = N - 1
243:                I1 = 1
244:                I2 = 2
245:             END IF
246:             CALL ZUNMQR( SIDE, TRANS, MI, NI, NQ-1, A( 2, 1 ), LDA, TAU,
247:      $                   C( I1, I2 ), LDC, WORK, LWORK, IINFO )
248:          END IF
249:       ELSE
250: *
251: *        Apply P
252: *
253:          IF( NOTRAN ) THEN
254:             TRANST = 'C'
255:          ELSE
256:             TRANST = 'N'
257:          END IF
258:          IF( NQ.GT.K ) THEN
259: *
260: *           P was determined by a call to ZGEBRD with nq > k
261: *
262:             CALL ZUNMLQ( SIDE, TRANST, M, N, K, A, LDA, TAU, C, LDC,
263:      $                   WORK, LWORK, IINFO )
264:          ELSE IF( NQ.GT.1 ) THEN
265: *
266: *           P was determined by a call to ZGEBRD with nq <= k
267: *
268:             IF( LEFT ) THEN
269:                MI = M - 1
270:                NI = N
271:                I1 = 2
272:                I2 = 1
273:             ELSE
274:                MI = M
275:                NI = N - 1
276:                I1 = 1
277:                I2 = 2
278:             END IF
279:             CALL ZUNMLQ( SIDE, TRANST, MI, NI, NQ-1, A( 1, 2 ), LDA,
280:      $                   TAU, C( I1, I2 ), LDC, WORK, LWORK, IINFO )
281:          END IF
282:       END IF
283:       WORK( 1 ) = LWKOPT
284:       RETURN
285: *
286: *     End of ZUNMBR
287: *
288:       END
289: