LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dorm22.f
Go to the documentation of this file.
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 unm22
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
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
subroutine dorm22(side, trans, m, n, n1, n2, q, ldq, c, ldc, work, lwork, info)
DORM22 multiplies a general matrix by a banded orthogonal matrix.
Definition dorm22.f:163