001:       SUBROUTINE DOPMTR( 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:       DOUBLE PRECISION   AP( * ), C( LDC, * ), TAU( * ), WORK( * )
015: *     ..
016: *
017: *  Purpose
018: *  =======
019: *
020: *  DOPMTR overwrites the general real M-by-N matrix C with
021: *
022: *                  SIDE = 'L'     SIDE = 'R'
023: *  TRANS = 'N':      Q * C          C * Q
024: *  TRANS = 'T':      Q**T * C       C * Q**T
025: *
026: *  where Q is a real orthogonal 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 DSPTRD 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**T from the Left;
040: *          = 'R': apply Q or Q**T from the Right.
041: *
042: *  UPLO    (input) CHARACTER*1
043: *          = 'U': Upper triangular packed storage used in previous
044: *                 call to DSPTRD;
045: *          = 'L': Lower triangular packed storage used in previous
046: *                 call to DSPTRD.
047: *
048: *  TRANS   (input) CHARACTER*1
049: *          = 'N':  No transpose, apply Q;
050: *          = 'T':  Transpose, apply Q**T.
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) DOUBLE PRECISION 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 DSPTRD.  AP is modified by the routine but
063: *          restored on exit.
064: *
065: *  TAU     (input) DOUBLE PRECISION 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 DSPTRD.
069: *
070: *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
071: *          On entry, the M-by-N matrix C.
072: *          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
073: *
074: *  LDC     (input) INTEGER
075: *          The leading dimension of the array C. LDC >= max(1,M).
076: *
077: *  WORK    (workspace) DOUBLE PRECISION 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:       DOUBLE PRECISION   ONE
089:       PARAMETER          ( ONE = 1.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:       DOUBLE PRECISION   AII
095: *     ..
096: *     .. External Functions ..
097:       LOGICAL            LSAME
098:       EXTERNAL           LSAME
099: *     ..
100: *     .. External Subroutines ..
101:       EXTERNAL           DLARF, XERBLA
102: *     ..
103: *     .. Intrinsic Functions ..
104:       INTRINSIC          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, 'T' ) ) 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( 'DOPMTR', -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 DSPTRD 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) is applied to C(1:i,1:n)
174: *
175:                MI = I
176:             ELSE
177: *
178: *              H(i) is applied to C(1:m,1:i)
179: *
180:                NI = I
181:             END IF
182: *
183: *           Apply H(i)
184: *
185:             AII = AP( II )
186:             AP( II ) = ONE
187:             CALL DLARF( SIDE, MI, NI, AP( II-I+1 ), 1, TAU( I ), C, LDC,
188:      $                  WORK )
189:             AP( II ) = AII
190: *
191:             IF( FORWRD ) THEN
192:                II = II + I + 2
193:             ELSE
194:                II = II - I - 1
195:             END IF
196:    10    CONTINUE
197:       ELSE
198: *
199: *        Q was determined by a call to DSPTRD with UPLO = 'L'.
200: *
201:          FORWRD = ( LEFT .AND. .NOT.NOTRAN ) .OR.
202:      $            ( .NOT.LEFT .AND. NOTRAN )
203: *
204:          IF( FORWRD ) THEN
205:             I1 = 1
206:             I2 = NQ - 1
207:             I3 = 1
208:             II = 2
209:          ELSE
210:             I1 = NQ - 1
211:             I2 = 1
212:             I3 = -1
213:             II = NQ*( NQ+1 ) / 2 - 1
214:          END IF
215: *
216:          IF( LEFT ) THEN
217:             NI = N
218:             JC = 1
219:          ELSE
220:             MI = M
221:             IC = 1
222:          END IF
223: *
224:          DO 20 I = I1, I2, I3
225:             AII = AP( II )
226:             AP( II ) = ONE
227:             IF( LEFT ) THEN
228: *
229: *              H(i) is applied to C(i+1:m,1:n)
230: *
231:                MI = M - I
232:                IC = I + 1
233:             ELSE
234: *
235: *              H(i) is applied to C(1:m,i+1:n)
236: *
237:                NI = N - I
238:                JC = I + 1
239:             END IF
240: *
241: *           Apply H(i)
242: *
243:             CALL DLARF( SIDE, MI, NI, AP( II ), 1, TAU( I ),
244:      $                  C( IC, JC ), LDC, WORK )
245:             AP( II ) = AII
246: *
247:             IF( FORWRD ) THEN
248:                II = II + NQ - I + 1
249:             ELSE
250:                II = II - NQ + I - 2
251:             END IF
252:    20    CONTINUE
253:       END IF
254:       RETURN
255: *
256: *     End of DOPMTR
257: *
258:       END
259: