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: *> \date January 2015
159: *
160: *> \ingroup complexOTHERcomputational
161: *
162: * =====================================================================
163: SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
164: $ WORK, LWORK, INFO )
165: *
166: * -- LAPACK computational routine (version 3.6.0) --
167: * -- LAPACK is a software package provided by Univ. of Tennessee, --
168: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169: * January 2015
170: *
171: IMPLICIT NONE
172: *
173: * .. Scalar Arguments ..
174: CHARACTER SIDE, TRANS
175: INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
176: * ..
177: * .. Array Arguments ..
178: DOUBLE PRECISION Q( LDQ, * ), C( LDC, * ), WORK( * )
179: * ..
180: *
181: * =====================================================================
182: *
183: * .. Parameters ..
184: DOUBLE PRECISION ONE
185: PARAMETER ( ONE = 1.0D+0 )
186: *
187: * .. Local Scalars ..
188: LOGICAL LEFT, LQUERY, NOTRAN
189: INTEGER I, LDWORK, LEN, LWKOPT, NB, NQ, NW
190: * ..
191: * .. External Functions ..
192: LOGICAL LSAME
193: EXTERNAL LSAME
194: * ..
195: * .. External Subroutines ..
196: EXTERNAL DGEMM, DLACPY, DTRMM, XERBLA
197: * ..
198: * .. Intrinsic Functions ..
199: INTRINSIC DBLE, MAX, MIN
200: * ..
201: * .. Executable Statements ..
202: *
203: * Test the input arguments
204: *
205: INFO = 0
206: LEFT = LSAME( SIDE, 'L' )
207: NOTRAN = LSAME( TRANS, 'N' )
208: LQUERY = ( LWORK.EQ.-1 )
209: *
210: * NQ is the order of Q;
211: * NW is the minimum dimension of WORK.
212: *
213: IF( LEFT ) THEN
214: NQ = M
215: ELSE
216: NQ = N
217: END IF
218: NW = NQ
219: IF( N1.EQ.0 .OR. N2.EQ.0 ) NW = 1
220: IF( .NOT.LEFT .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
221: INFO = -1
222: ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.LSAME( TRANS, 'T' ) )
223: $ THEN
224: INFO = -2
225: ELSE IF( M.LT.0 ) THEN
226: INFO = -3
227: ELSE IF( N.LT.0 ) THEN
228: INFO = -4
229: ELSE IF( N1.LT.0 .OR. N1+N2.NE.NQ ) THEN
230: INFO = -5
231: ELSE IF( N2.LT.0 ) THEN
232: INFO = -6
233: ELSE IF( LDQ.LT.MAX( 1, NQ ) ) THEN
234: INFO = -8
235: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
236: INFO = -10
237: ELSE IF( LWORK.LT.NW .AND. .NOT.LQUERY ) THEN
238: INFO = -12
239: END IF
240: *
241: IF( INFO.EQ.0 ) THEN
242: LWKOPT = M*N
243: WORK( 1 ) = DBLE( LWKOPT )
244: END IF
245: *
246: IF( INFO.NE.0 ) THEN
247: CALL XERBLA( 'DORM22', -INFO )
248: RETURN
249: ELSE IF( LQUERY ) THEN
250: RETURN
251: END IF
252: *
253: * Quick return if possible
254: *
255: IF( M.EQ.0 .OR. N.EQ.0 ) THEN
256: WORK( 1 ) = 1
257: RETURN
258: END IF
259: *
260: * Degenerate cases (N1 = 0 or N2 = 0) are handled using DTRMM.
261: *
262: IF( N1.EQ.0 ) THEN
263: CALL DTRMM( SIDE, 'Upper', TRANS, 'Non-Unit', M, N, ONE,
264: $ Q, LDQ, C, LDC )
265: WORK( 1 ) = ONE
266: RETURN
267: ELSE IF( N2.EQ.0 ) THEN
268: CALL DTRMM( SIDE, 'Lower', TRANS, 'Non-Unit', M, N, ONE,
269: $ Q, LDQ, C, LDC )
270: WORK( 1 ) = ONE
271: RETURN
272: END IF
273: *
274: * Compute the largest chunk size available from the workspace.
275: *
276: NB = MAX( 1, MIN( LWORK, LWKOPT ) / NQ )
277: *
278: IF( LEFT ) THEN
279: IF( NOTRAN ) THEN
280: DO I = 1, N, NB
281: LEN = MIN( NB, N-I+1 )
282: LDWORK = M
283: *
284: * Multiply bottom part of C by Q12.
285: *
286: CALL DLACPY( 'All', N1, LEN, C( N2+1, I ), LDC, WORK,
287: $ LDWORK )
288: CALL DTRMM( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
289: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ, WORK,
290: $ LDWORK )
291: *
292: * Multiply top part of C by Q11.
293: *
294: CALL DGEMM( 'No Transpose', 'No Transpose', N1, LEN, N2,
295: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
296: $ LDWORK )
297: *
298: * Multiply top part of C by Q21.
299: *
300: CALL DLACPY( 'All', N2, LEN, C( 1, I ), LDC,
301: $ WORK( N1+1 ), LDWORK )
302: CALL DTRMM( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
303: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ,
304: $ WORK( N1+1 ), LDWORK )
305: *
306: * Multiply bottom part of C by Q22.
307: *
308: CALL DGEMM( 'No Transpose', 'No Transpose', N2, LEN, N1,
309: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N2+1, I ), LDC,
310: $ ONE, WORK( N1+1 ), LDWORK )
311: *
312: * Copy everything back.
313: *
314: CALL DLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
315: $ LDC )
316: END DO
317: ELSE
318: DO I = 1, N, NB
319: LEN = MIN( NB, N-I+1 )
320: LDWORK = M
321: *
322: * Multiply bottom part of C by Q21**T.
323: *
324: CALL DLACPY( 'All', N2, LEN, C( N1+1, I ), LDC, WORK,
325: $ LDWORK )
326: CALL DTRMM( 'Left', 'Upper', 'Transpose', 'Non-Unit',
327: $ N2, LEN, ONE, Q( N1+1, 1 ), LDQ, WORK,
328: $ LDWORK )
329: *
330: * Multiply top part of C by Q11**T.
331: *
332: CALL DGEMM( 'Transpose', 'No Transpose', N2, LEN, N1,
333: $ ONE, Q, LDQ, C( 1, I ), LDC, ONE, WORK,
334: $ LDWORK )
335: *
336: * Multiply top part of C by Q12**T.
337: *
338: CALL DLACPY( 'All', N1, LEN, C( 1, I ), LDC,
339: $ WORK( N2+1 ), LDWORK )
340: CALL DTRMM( 'Left', 'Lower', 'Transpose', 'Non-Unit',
341: $ N1, LEN, ONE, Q( 1, N2+1 ), LDQ,
342: $ WORK( N2+1 ), LDWORK )
343: *
344: * Multiply bottom part of C by Q22**T.
345: *
346: CALL DGEMM( 'Transpose', 'No Transpose', N1, LEN, N2,
347: $ ONE, Q( N1+1, N2+1 ), LDQ, C( N1+1, I ), LDC,
348: $ ONE, WORK( N2+1 ), LDWORK )
349: *
350: * Copy everything back.
351: *
352: CALL DLACPY( 'All', M, LEN, WORK, LDWORK, C( 1, I ),
353: $ LDC )
354: END DO
355: END IF
356: ELSE
357: IF( NOTRAN ) THEN
358: DO I = 1, M, NB
359: LEN = MIN( NB, M-I+1 )
360: LDWORK = LEN
361: *
362: * Multiply right part of C by Q21.
363: *
364: CALL DLACPY( 'All', LEN, N2, C( I, N1+1 ), LDC, WORK,
365: $ LDWORK )
366: CALL DTRMM( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
367: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ, WORK,
368: $ LDWORK )
369: *
370: * Multiply left part of C by Q11.
371: *
372: CALL DGEMM( 'No Transpose', 'No Transpose', LEN, N2, N1,
373: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
374: $ LDWORK )
375: *
376: * Multiply left part of C by Q12.
377: *
378: CALL DLACPY( 'All', LEN, N1, C( I, 1 ), LDC,
379: $ WORK( 1 + N2*LDWORK ), LDWORK )
380: CALL DTRMM( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
381: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ,
382: $ WORK( 1 + N2*LDWORK ), LDWORK )
383: *
384: * Multiply right part of C by Q22.
385: *
386: CALL DGEMM( 'No Transpose', 'No Transpose', LEN, N1, N2,
387: $ ONE, C( I, N1+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
388: $ ONE, WORK( 1 + N2*LDWORK ), LDWORK )
389: *
390: * Copy everything back.
391: *
392: CALL DLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
393: $ LDC )
394: END DO
395: ELSE
396: DO I = 1, M, NB
397: LEN = MIN( NB, M-I+1 )
398: LDWORK = LEN
399: *
400: * Multiply right part of C by Q12**T.
401: *
402: CALL DLACPY( 'All', LEN, N1, C( I, N2+1 ), LDC, WORK,
403: $ LDWORK )
404: CALL DTRMM( 'Right', 'Lower', 'Transpose', 'Non-Unit',
405: $ LEN, N1, ONE, Q( 1, N2+1 ), LDQ, WORK,
406: $ LDWORK )
407: *
408: * Multiply left part of C by Q11**T.
409: *
410: CALL DGEMM( 'No Transpose', 'Transpose', LEN, N1, N2,
411: $ ONE, C( I, 1 ), LDC, Q, LDQ, ONE, WORK,
412: $ LDWORK )
413: *
414: * Multiply left part of C by Q21**T.
415: *
416: CALL DLACPY( 'All', LEN, N2, C( I, 1 ), LDC,
417: $ WORK( 1 + N1*LDWORK ), LDWORK )
418: CALL DTRMM( 'Right', 'Upper', 'Transpose', 'Non-Unit',
419: $ LEN, N2, ONE, Q( N1+1, 1 ), LDQ,
420: $ WORK( 1 + N1*LDWORK ), LDWORK )
421: *
422: * Multiply right part of C by Q22**T.
423: *
424: CALL DGEMM( 'No Transpose', 'Transpose', LEN, N2, N1,
425: $ ONE, C( I, N2+1 ), LDC, Q( N1+1, N2+1 ), LDQ,
426: $ ONE, WORK( 1 + N1*LDWORK ), LDWORK )
427: *
428: * Copy everything back.
429: *
430: CALL DLACPY( 'All', LEN, N, WORK, LDWORK, C( I, 1 ),
431: $ LDC )
432: END DO
433: END IF
434: END IF
435: *
436: WORK( 1 ) = DBLE( LWKOPT )
437: RETURN
438: *
439: * End of DORM22
440: *
441: 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