001:       SUBROUTINE ZUPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK,
002:      $                   INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          SIDE, TRANS, UPLO
011:       INTEGER            INFO, LDC, M, N
012: *     ..
013: *     .. Array Arguments ..
014:       COMPLEX*16         AP( * ), C( LDC, * ), TAU( * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  ZUPMTR overwrites the general complex M-by-N matrix C with
021: *
022: *                  SIDE = 'L'     SIDE = 'R'
023: *  TRANS = 'N':      Q * C          C * Q
024: *  TRANS = 'C':      Q**H * C       C * Q**H
025: *
026: *  where Q is a complex unitary matrix of order nq, with nq = m if
027: *  SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
028: *  nq-1 elementary reflectors, as returned by ZHPTRD using packed
029: *  storage:
030: *
031: *  if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
032: *
033: *  if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
034: *
035: *  Arguments
036: *  =========
037: *
038: *  SIDE    (input) CHARACTER*1
039: *          = 'L': apply Q or Q**H from the Left;
040: *          = 'R': apply Q or Q**H from the Right.
041: *
042: *  UPLO    (input) CHARACTER*1
043: *          = 'U': Upper triangular packed storage used in previous
044: *                 call to ZHPTRD;
045: *          = 'L': Lower triangular packed storage used in previous
046: *                 call to ZHPTRD.
047: *
048: *  TRANS   (input) CHARACTER*1
049: *          = 'N':  No transpose, apply Q;
050: *          = 'C':  Conjugate transpose, apply Q**H.
051: *
052: *  M       (input) INTEGER
053: *          The number of rows of the matrix C. M >= 0.
054: *
055: *  N       (input) INTEGER
056: *          The number of columns of the matrix C. N >= 0.
057: *
058: *  AP      (input) COMPLEX*16 array, dimension
059: *                               (M*(M+1)/2) if SIDE = 'L'
060: *                               (N*(N+1)/2) if SIDE = 'R'
061: *          The vectors which define the elementary reflectors, as
062: *          returned by ZHPTRD.  AP is modified by the routine but
063: *          restored on exit.
064: *
065: *  TAU     (input) COMPLEX*16 array, dimension (M-1) if SIDE = 'L'
066: *                                     or (N-1) if SIDE = 'R'
067: *          TAU(i) must contain the scalar factor of the elementary
068: *          reflector H(i), as returned by ZHPTRD.
069: *
070: *  C       (input/output) COMPLEX*16 array, dimension (LDC,N)
071: *          On entry, the M-by-N matrix C.
072: *          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
073: *
074: *  LDC     (input) INTEGER
075: *          The leading dimension of the array C. LDC >= max(1,M).
076: *
077: *  WORK    (workspace) COMPLEX*16 array, dimension
078: *                                   (N) if SIDE = 'L'
079: *                                   (M) if SIDE = 'R'
080: *
081: *  INFO    (output) INTEGER
082: *          = 0:  successful exit
083: *          < 0:  if INFO = -i, the i-th argument had an illegal value
084: *
085: *  =====================================================================
086: *
087: *     .. Parameters ..
088:       COMPLEX*16         ONE
089:       PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ) )
090: *     ..
091: *     .. Local Scalars ..
092:       LOGICAL            FORWRD, LEFT, NOTRAN, UPPER
093:       INTEGER            I, I1, I2, I3, IC, II, JC, MI, NI, NQ
094:       COMPLEX*16         AII, TAUI
095: *     ..
096: *     .. External Functions ..
097:       LOGICAL            LSAME
098:       EXTERNAL           LSAME
099: *     ..
100: *     .. External Subroutines ..
101:       EXTERNAL           XERBLA, ZLARF
102: *     ..
103: *     .. Intrinsic Functions ..
104:       INTRINSIC          DCONJG, MAX
105: *     ..
106: *     .. Executable Statements ..
107: *
108: *     Test the input arguments
109: *
110:       INFO = 0
111:       LEFT = LSAME( SIDE, 'L' )
112:       NOTRAN = LSAME( TRANS, 'N' )
113:       UPPER = LSAME( UPLO, 'U' )
114: *
115: *     NQ is the order of Q
116: *
117:       IF( LEFT ) THEN
118:          NQ = M
119:       ELSE
120:          NQ = N
121:       END IF
122:       IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
123:          INFO = -1
124:       ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
125:          INFO = -2
126:       ELSE IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
127:          INFO = -3
128:       ELSE IF( M.LT.0 ) THEN
129:          INFO = -4
130:       ELSE IF( N.LT.0 ) THEN
131:          INFO = -5
132:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
133:          INFO = -9
134:       END IF
135:       IF( INFO.NE.0 ) THEN
136:          CALL XERBLA( 'ZUPMTR', -INFO )
137:          RETURN
138:       END IF
139: *
140: *     Quick return if possible
141: *
142:       IF( M.EQ.0 .OR. N.EQ.0 )
143:      $   RETURN
144: *
145:       IF( UPPER ) THEN
146: *
147: *        Q was determined by a call to ZHPTRD with UPLO = 'U'
148: *
149:          FORWRD = ( LEFT .AND. NOTRAN ) .OR.
150:      $            ( .NOT.LEFT .AND. .NOT.NOTRAN )
151: *
152:          IF( FORWRD ) THEN
153:             I1 = 1
154:             I2 = NQ - 1
155:             I3 = 1
156:             II = 2
157:          ELSE
158:             I1 = NQ - 1
159:             I2 = 1
160:             I3 = -1
161:             II = NQ*( NQ+1 ) / 2 - 1
162:          END IF
163: *
164:          IF( LEFT ) THEN
165:             NI = N
166:          ELSE
167:             MI = M
168:          END IF
169: *
170:          DO 10 I = I1, I2, I3
171:             IF( LEFT ) THEN
172: *
173: *              H(i) or H(i)' is applied to C(1:i,1:n)
174: *
175:                MI = I
176:             ELSE
177: *
178: *              H(i) or H(i)' is applied to C(1:m,1:i)
179: *
180:                NI = I
181:             END IF
182: *
183: *           Apply H(i) or H(i)'
184: *
185:             IF( NOTRAN ) THEN
186:                TAUI = TAU( I )
187:             ELSE
188:                TAUI = DCONJG( TAU( I ) )
189:             END IF
190:             AII = AP( II )
191:             AP( II ) = ONE
192:             CALL ZLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAUI, C, LDC,
193:      $                  WORK )
194:             AP( II ) = AII
195: *
196:             IF( FORWRD ) THEN
197:                II = II + I + 2
198:             ELSE
199:                II = II - I - 1
200:             END IF
201:    10    CONTINUE
202:       ELSE
203: *
204: *        Q was determined by a call to ZHPTRD with UPLO = 'L'.
205: *
206:          FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR.
207:      $            ( .NOT.LEFT .AND. NOTRAN )
208: *
209:          IF( FORWRD ) THEN
210:             I1 = 1
211:             I2 = NQ - 1
212:             I3 = 1
213:             II = 2
214:          ELSE
215:             I1 = NQ - 1
216:             I2 = 1
217:             I3 = -1
218:             II = NQ*( NQ+1 ) / 2 - 1
219:          END IF
220: *
221:          IF( LEFT ) THEN
222:             NI = N
223:             JC = 1
224:          ELSE
225:             MI = M
226:             IC = 1
227:          END IF
228: *
229:          DO 20 I = I1, I2, I3
230:             AII = AP( II )
231:             AP( II ) = ONE
232:             IF( LEFT ) THEN
233: *
234: *              H(i) or H(i)' is applied to C(i+1:m,1:n)
235: *
236:                MI = M - I
237:                IC = I + 1
238:             ELSE
239: *
240: *              H(i) or H(i)' is applied to C(1:m,i+1:n)
241: *
242:                NI = N - I
243:                JC = I + 1
244:             END IF
245: *
246: *           Apply H(i) or H(i)'
247: *
248:             IF( NOTRAN ) THEN
249:                TAUI = TAU( I )
250:             ELSE
251:                TAUI = DCONJG( TAU( I ) )
252:             END IF
253:             CALL ZLARF( SIDE, MI, NI, AP( II ), 1, TAUI, C( IC, JC ),
254:      $                  LDC, WORK )
255:             AP( II ) = AII
256: *
257:             IF( FORWRD ) THEN
258:                II = II + NQ - I + 1
259:             ELSE
260:                II = II - NQ + I - 2
261:             END IF
262:    20    CONTINUE
263:       END IF
264:       RETURN
265: *
266: *     End of ZUPMTR
267: *
268:       END
269: