1: *> \brief \b DORM22 multiplies a general matrix by a banded orthogonal matrix.
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DORM22 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorm22.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorm22.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorm22.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
22: * $ WORK, LWORK, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER SIDE, TRANS
26: * INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
27: * ..
28: * .. Array Arguments ..
29: * DOUBLE PRECISION Q( LDQ, * ), C( LDC, * ), WORK( * )
30: * ..
31: *
32: *> \par Purpose
33: * ============
34: *>
35: *> \verbatim
36: *>
37: *>
38: *> DORM22 overwrites the general real M-by-N matrix C with
39: *>
40: *> SIDE = 'L' SIDE = 'R'
41: *> TRANS = 'N': Q * C C * Q
42: *> TRANS = 'T': Q**T * C C * Q**T
43: *>
44: *> where Q is a real orthogonal matrix of order NQ, with NQ = M if
45: *> SIDE = 'L' and NQ = N if SIDE = 'R'.
46: *> The orthogonal matrix Q processes a 2-by-2 block structure
47: *>
48: *> [ Q11 Q12 ]
49: *> Q = [ ]
50: *> [ Q21 Q22 ],
51: *>
52: *> where Q12 is an N1-by-N1 lower triangular matrix and Q21 is an
53: *> N2-by-N2 upper triangular matrix.
54: *> \endverbatim
55: *
56: * Arguments:
57: * ==========
58: *
59: *> \param[in] SIDE
60: *> \verbatim
61: *> SIDE is CHARACTER*1
62: *> = 'L': apply Q or Q**T from the Left;
63: *> = 'R': apply Q or Q**T from the Right.
64: *> \endverbatim
65: *>
66: *> \param[in] TRANS
67: *> \verbatim
68: *> TRANS is CHARACTER*1
69: *> = 'N': apply Q (No transpose);
70: *> = 'C': apply Q**T (Conjugate transpose).
71: *> \endverbatim
72: *>
73: *> \param[in] M
74: *> \verbatim
75: *> M is INTEGER
76: *> The number of rows of the matrix C. M >= 0.
77: *> \endverbatim
78: *>
79: *> \param[in] N
80: *> \verbatim
81: *> N is INTEGER
82: *> The number of columns of the matrix C. N >= 0.
83: *> \endverbatim
84: *>
85: *> \param[in] N1
86: *> \param[in] N2
87: *> \verbatim
88: *> N1 is INTEGER
89: *> N2 is INTEGER
90: *> The dimension of Q12 and Q21, respectively. N1, N2 >= 0.
91: *> The following requirement must be satisfied:
92: *> N1 + N2 = M if SIDE = 'L' and N1 + N2 = N if SIDE = 'R'.
93: *> \endverbatim
94: *>
95: *> \param[in] Q
96: *> \verbatim
97: *> Q is DOUBLE PRECISION array, dimension
98: *> (LDQ,M) if SIDE = 'L'
99: *> (LDQ,N) if SIDE = 'R'
100: *> \endverbatim
101: *>
102: *> \param[in] LDQ
103: *> \verbatim
104: *> LDQ is INTEGER
105: *> The leading dimension of the array Q.
106: *> LDQ >= max(1,M) if SIDE = 'L'; LDQ >= max(1,N) if SIDE = 'R'.
107: *> \endverbatim
108: *>
109: *> \param[in,out] C
110: *> \verbatim
111: *> C is DOUBLE PRECISION array, dimension (LDC,N)
112: *> On entry, the M-by-N matrix C.
113: *> On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
114: *> \endverbatim
115: *>
116: *> \param[in] LDC
117: *> \verbatim
118: *> LDC is INTEGER
119: *> The leading dimension of the array C. LDC >= max(1,M).
120: *> \endverbatim
121: *>
122: *> \param[out] WORK
123: *> \verbatim
124: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
125: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
126: *> \endverbatim
127: *>
128: *> \param[in] LWORK
129: *> \verbatim
130: *> LWORK is INTEGER
131: *> The dimension of the array WORK.
132: *> If SIDE = 'L', LWORK >= max(1,N);
133: *> if SIDE = 'R', LWORK >= max(1,M).
134: *> For optimum performance LWORK >= M*N.
135: *>
136: *> If LWORK = -1, then a workspace query is assumed; the routine
137: *> only calculates the optimal size of the WORK array, returns
138: *> this value as the first entry of the WORK array, and no error
139: *> message related to LWORK is issued by XERBLA.
140: *> \endverbatim
141: *>
142: *> \param[out] INFO
143: *> \verbatim
144: *> INFO is INTEGER
145: *> = 0: successful exit
146: *> < 0: if INFO = -i, the i-th argument had an illegal value
147: *> \endverbatim
148: *
149: *
150: * Authors:
151: * ========
152: *
153: *> \author Univ. of Tennessee
154: *> \author Univ. of California Berkeley
155: *> \author Univ. of Colorado Denver
156: *> \author NAG Ltd.
157: *
158: *> \ingroup complexOTHERcomputational
159: *
160: * =====================================================================
161: SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
162: $ WORK, LWORK, INFO )
163: *
164: * -- LAPACK computational routine --
165: * -- LAPACK is a software package provided by Univ. of Tennessee, --
166: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
167: *
168: IMPLICIT NONE
169: *
170: * .. Scalar Arguments ..
171: CHARACTER SIDE, TRANS
172: INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
173: * ..
174: * .. Array Arguments ..
175: DOUBLE PRECISION Q( LDQ, * ), C( LDC, * ), WORK( * )
176: * ..
177: *
178: * =====================================================================
179: *
180: * .. Parameters ..
181: DOUBLE PRECISION ONE
182: PARAMETER ( ONE = 1.0D+0 )
183: *
184: * .. Local Scalars ..
185: LOGICAL LEFT, LQUERY, NOTRAN
186: INTEGER I, LDWORK, LEN, LWKOPT, NB, NQ, NW
187: * ..
188: * .. External Functions ..
189: LOGICAL LSAME
190: EXTERNAL LSAME
191: * ..
192: * .. External Subroutines ..
193: EXTERNAL DGEMM, DLACPY, DTRMM, XERBLA
194: * ..
195: * .. Intrinsic Functions ..
196: INTRINSIC DBLE, MAX, MIN
197: * ..
198: * .. Executable Statements ..
199: *
200: * Test the input arguments
201: *
202: INFO = 0
203: LEFT = LSAME( SIDE, 'L' )
204: NOTRAN = LSAME( TRANS, 'N' )
205: LQUERY = ( LWORK.EQ.-1 )
206: *
207: * NQ is the order of Q;
208: * NW is the minimum dimension of WORK.
209: *
210: IF( LEFT ) THEN
211: NQ = M
212: ELSE
213: NQ = N
214: END IF
215: NW = NQ
216: IF( N1.EQ.0 .OR. N2.EQ.0 ) NW = 1
217: IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
218: INFO = -1
219: ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'T' ) )
220: $ THEN
221: INFO = -2
222: ELSE IF( M.LT.0 ) THEN
223: INFO = -3
224: ELSE IF( N.LT.0 ) THEN
225: INFO = -4
226: ELSE IF( N1.LT.0 .OR. N1+N2.NE.NQ ) THEN
227: INFO = -5
228: ELSE IF( N2.LT.0 ) THEN
229: INFO = -6
230: ELSE IF( LDQ.LT.MAX( 1, NQ ) ) THEN
231: INFO = -8
232: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
233: INFO = -10
234: ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
235: INFO = -12
236: END IF
237: *
238: IF( INFO.EQ.0 ) THEN
239: LWKOPT = M*N
240: WORK( 1 ) = DBLE( LWKOPT )
241: END IF
242: *
243: IF( INFO.NE.0 ) THEN
244: CALL XERBLA( 'DORM22', -INFO )
245: RETURN
246: ELSE IF( LQUERY ) THEN
247: RETURN
248: END IF
249: *
250: * Quick return if possible
251: *
252: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
253: WORK( 1 ) = 1
254: RETURN
255: END IF
256: *
257: * Degenerate cases (N1 = 0 or N2 = 0) are handled using DTRMM.
258: *
259: IF( N1.EQ.0 ) THEN
260: CALL DTRMM( SIDE, 'Upper', TRANS, 'Non-Unit', M, N, ONE,
261: $ Q, LDQ, C, LDC )
262: WORK( 1 ) = ONE
263: RETURN
264: ELSE IF( N2.EQ.0 ) THEN
265: CALL DTRMM( SIDE, 'Lower', TRANS, 'Non-Unit', M, N, ONE,
266: $ Q, LDQ, C, LDC )
267: WORK( 1 ) = ONE
268: RETURN
269: END IF
270: *
271: * Compute the largest chunk size available from the workspace.
272: *
273: NB = MAX( 1, MIN( LWORK, LWKOPT ) / NQ )
274: *
275: IF( LEFT ) THEN
276: IF( NOTRAN ) THEN
277: DO I = 1, N, NB
278: LEN = MIN( NB, N-I+1 )
279: LDWORK = M
280: *
281: * Multiply bottom part of C by Q12.
282: *
283: CALL DLACPY( 'All', N1, LEN, C( N2+1, I ), LDC, WORK,
284: $ LDWORK )
285: CALL DTRMM( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
286: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ, WORK,
287: $ LDWORK )
288: *
289: * Multiply top part of C by Q11.
290: *
291: CALL DGEMM( 'No Transpose', 'No Transpose', N1, LEN, N2,
292: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
293: $ LDWORK )
294: *
295: * Multiply top part of C by Q21.
296: *
297: CALL DLACPY( 'All', N2, LEN, C( 1, I ), LDC,
298: $ WORK( N1+1 ), LDWORK )
299: CALL DTRMM( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
300: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ,
301: $ WORK( N1+1 ), LDWORK )
302: *
303: * Multiply bottom part of C by Q22.
304: *
305: CALL DGEMM( 'No Transpose', 'No Transpose', N2, LEN, N1,
306: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N2+1, I ), LDC,
307: $ ONE, WORK( N1+1 ), LDWORK )
308: *
309: * Copy everything back.
310: *
311: CALL DLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
312: $ LDC )
313: END DO
314: ELSE
315: DO I = 1, N, NB
316: LEN = MIN( NB, N-I+1 )
317: LDWORK = M
318: *
319: * Multiply bottom part of C by Q21**T.
320: *
321: CALL DLACPY( 'All', N2, LEN, C( N1+1, I ), LDC, WORK,
322: $ LDWORK )
323: CALL DTRMM( 'Left', 'Upper', 'Transpose', 'Non-Unit',
324: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ, WORK,
325: $ LDWORK )
326: *
327: * Multiply top part of C by Q11**T.
328: *
329: CALL DGEMM( 'Transpose', 'No Transpose', N2, LEN, N1,
330: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
331: $ LDWORK )
332: *
333: * Multiply top part of C by Q12**T.
334: *
335: CALL DLACPY( 'All', N1, LEN, C( 1, I ), LDC,
336: $ WORK( N2+1 ), LDWORK )
337: CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-Unit',
338: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ,
339: $ WORK( N2+1 ), LDWORK )
340: *
341: * Multiply bottom part of C by Q22**T.
342: *
343: CALL DGEMM( 'Transpose', 'No Transpose', N1, LEN, N2,
344: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N1+1, I ), LDC,
345: $ ONE, WORK( N2+1 ), LDWORK )
346: *
347: * Copy everything back.
348: *
349: CALL DLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
350: $ LDC )
351: END DO
352: END IF
353: ELSE
354: IF( NOTRAN ) THEN
355: DO I = 1, M, NB
356: LEN = MIN( NB, M-I+1 )
357: LDWORK = LEN
358: *
359: * Multiply right part of C by Q21.
360: *
361: CALL DLACPY( 'All', LEN, N2, C( I, N1+1 ), LDC, WORK,
362: $ LDWORK )
363: CALL DTRMM( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
364: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ, WORK,
365: $ LDWORK )
366: *
367: * Multiply left part of C by Q11.
368: *
369: CALL DGEMM( 'No Transpose', 'No Transpose', LEN, N2, N1,
370: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
371: $ LDWORK )
372: *
373: * Multiply left part of C by Q12.
374: *
375: CALL DLACPY( 'All', LEN, N1, C( I, 1 ), LDC,
376: $ WORK( 1 + N2*LDWORK ), LDWORK )
377: CALL DTRMM( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
378: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ,
379: $ WORK( 1 + N2*LDWORK ), LDWORK )
380: *
381: * Multiply right part of C by Q22.
382: *
383: CALL DGEMM( 'No Transpose', 'No Transpose', LEN, N1, N2,
384: $ ONE, C( I, N1+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
385: $ ONE, WORK( 1 + N2*LDWORK ), LDWORK )
386: *
387: * Copy everything back.
388: *
389: CALL DLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
390: $ LDC )
391: END DO
392: ELSE
393: DO I = 1, M, NB
394: LEN = MIN( NB, M-I+1 )
395: LDWORK = LEN
396: *
397: * Multiply right part of C by Q12**T.
398: *
399: CALL DLACPY( 'All', LEN, N1, C( I, N2+1 ), LDC, WORK,
400: $ LDWORK )
401: CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Non-Unit',
402: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ, WORK,
403: $ LDWORK )
404: *
405: * Multiply left part of C by Q11**T.
406: *
407: CALL DGEMM( 'No Transpose', 'Transpose', LEN, N1, N2,
408: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
409: $ LDWORK )
410: *
411: * Multiply left part of C by Q21**T.
412: *
413: CALL DLACPY( 'All', LEN, N2, C( I, 1 ), LDC,
414: $ WORK( 1 + N1*LDWORK ), LDWORK )
415: CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-Unit',
416: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ,
417: $ WORK( 1 + N1*LDWORK ), LDWORK )
418: *
419: * Multiply right part of C by Q22**T.
420: *
421: CALL DGEMM( 'No Transpose', 'Transpose', LEN, N2, N1,
422: $ ONE, C( I, N2+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
423: $ ONE, WORK( 1 + N1*LDWORK ), LDWORK )
424: *
425: * Copy everything back.
426: *
427: CALL DLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
428: $ LDC )
429: END DO
430: END IF
431: END IF
432: *
433: WORK( 1 ) = DBLE( LWKOPT )
434: RETURN
435: *
436: * End of DORM22
437: *
438: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dorm22.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack