![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dopmtr.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dopmtr.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DOPMTR( SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK,
2: $ 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, UPLO
11: INTEGER INFO, LDC, M, N
12: * ..
13: * .. Array Arguments ..
14: DOUBLE PRECISION AP( * ), C( LDC, * ), TAU( * ), WORK( * )
15: * ..
16: *
17: * Purpose
18: * =======
19: *
20: * DOPMTR overwrites the general real M-by-N matrix C with
21: *
22: * SIDE = 'L' SIDE = 'R'
23: * TRANS = 'N': Q * C C * Q
24: * TRANS = 'T': Q**T * C C * Q**T
25: *
26: * where Q is a real orthogonal matrix of order nq, with nq = m if
27: * SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
28: * nq-1 elementary reflectors, as returned by DSPTRD using packed
29: * storage:
30: *
31: * if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
32: *
33: * if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
34: *
35: * Arguments
36: * =========
37: *
38: * SIDE (input) CHARACTER*1
39: * = 'L': apply Q or Q**T from the Left;
40: * = 'R': apply Q or Q**T from the Right.
41: *
42: * UPLO (input) CHARACTER*1
43: * = 'U': Upper triangular packed storage used in previous
44: * call to DSPTRD;
45: * = 'L': Lower triangular packed storage used in previous
46: * call to DSPTRD.
47: *
48: * TRANS (input) CHARACTER*1
49: * = 'N': No transpose, apply Q;
50: * = 'T': Transpose, apply Q**T.
51: *
52: * M (input) INTEGER
53: * The number of rows of the matrix C. M >= 0.
54: *
55: * N (input) INTEGER
56: * The number of columns of the matrix C. N >= 0.
57: *
58: * AP (input) DOUBLE PRECISION array, dimension
59: * (M*(M+1)/2) if SIDE = 'L'
60: * (N*(N+1)/2) if SIDE = 'R'
61: * The vectors which define the elementary reflectors, as
62: * returned by DSPTRD. AP is modified by the routine but
63: * restored on exit.
64: *
65: * TAU (input) DOUBLE PRECISION array, dimension (M-1) if SIDE = 'L'
66: * or (N-1) if SIDE = 'R'
67: * TAU(i) must contain the scalar factor of the elementary
68: * reflector H(i), as returned by DSPTRD.
69: *
70: * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
71: * On entry, the M-by-N matrix C.
72: * On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
73: *
74: * LDC (input) INTEGER
75: * The leading dimension of the array C. LDC >= max(1,M).
76: *
77: * WORK (workspace) DOUBLE PRECISION array, dimension
78: * (N) if SIDE = 'L'
79: * (M) if SIDE = 'R'
80: *
81: * INFO (output) INTEGER
82: * = 0: successful exit
83: * < 0: if INFO = -i, the i-th argument had an illegal value
84: *
85: * =====================================================================
86: *
87: * .. Parameters ..
88: DOUBLE PRECISION ONE
89: PARAMETER ( ONE = 1.0D+0 )
90: * ..
91: * .. Local Scalars ..
92: LOGICAL FORWRD, LEFT, NOTRAN, UPPER
93: INTEGER I, I1, I2, I3, IC, II, JC, MI, NI, NQ
94: DOUBLE PRECISION AII
95: * ..
96: * .. External Functions ..
97: LOGICAL LSAME
98: EXTERNAL LSAME
99: * ..
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