LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
sorm22.f
Go to the documentation of this file.
1*> \brief \b SORM22 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 SORM22 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorm22.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorm22.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorm22.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE SORM22( 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* REAL Q( LDQ, * ), C( LDC, * ), WORK( * )
30* ..
31*
32*> \par Purpose
33* ============
34*>
35*> \verbatim
36*>
37*>
38*> SORM22 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 REAL 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 REAL 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 REAL 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 sorm22( 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 REAL Q( LDQ, * ), C( LDC, * ), WORK( * )
176* ..
177*
178* =====================================================================
179*
180* .. Parameters ..
181 REAL ONE
182 parameter( one = 1.0e+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 REAL SROUNDUP_LWORK
191 EXTERNAL lsame, sroundup_lwork
192* ..
193* .. External Subroutines ..
194 EXTERNAL sgemm, slacpy, strmm, xerbla
195* ..
196* .. Intrinsic Functions ..
197 INTRINSIC max, min
198* ..
199* .. Executable Statements ..
200*
201* Test the input arguments
202*
203 info = 0
204 left = lsame( side, 'L' )
205 notran = lsame( trans, 'N' )
206 lquery = ( lwork.EQ.-1 )
207*
208* NQ is the order of Q;
209* NW is the minimum dimension of WORK.
210*
211 IF( left ) THEN
212 nq = m
213 ELSE
214 nq = n
215 END IF
216 nw = nq
217 IF( n1.EQ.0 .OR. n2.EQ.0 ) nw = 1
218 IF( .NOT.left .AND. .NOT.lsame( side, 'R' ) ) THEN
219 info = -1
220 ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'T' ) )
221 $ THEN
222 info = -2
223 ELSE IF( m.LT.0 ) THEN
224 info = -3
225 ELSE IF( n.LT.0 ) THEN
226 info = -4
227 ELSE IF( n1.LT.0 .OR. n1+n2.NE.nq ) THEN
228 info = -5
229 ELSE IF( n2.LT.0 ) THEN
230 info = -6
231 ELSE IF( ldq.LT.max( 1, nq ) ) THEN
232 info = -8
233 ELSE IF( ldc.LT.max( 1, m ) ) THEN
234 info = -10
235 ELSE IF( lwork.LT.nw .AND. .NOT.lquery ) THEN
236 info = -12
237 END IF
238*
239 IF( info.EQ.0 ) THEN
240 lwkopt = m*n
241 work( 1 ) = sroundup_lwork( lwkopt )
242 END IF
243*
244 IF( info.NE.0 ) THEN
245 CALL xerbla( 'SORM22', -info )
246 RETURN
247 ELSE IF( lquery ) THEN
248 RETURN
249 END IF
250*
251* Quick return if possible
252*
253 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
254 work( 1 ) = 1
255 RETURN
256 END IF
257*
258* Degenerate cases (N1 = 0 or N2 = 0) are handled using STRMM.
259*
260 IF( n1.EQ.0 ) THEN
261 CALL strmm( side, 'Upper', trans, 'Non-Unit', m, n, one,
262 $ q, ldq, c, ldc )
263 work( 1 ) = one
264 RETURN
265 ELSE IF( n2.EQ.0 ) THEN
266 CALL strmm( side, 'Lower', trans, 'Non-Unit', m, n, one,
267 $ q, ldq, c, ldc )
268 work( 1 ) = one
269 RETURN
270 END IF
271*
272* Compute the largest chunk size available from the workspace.
273*
274 nb = max( 1, min( lwork, lwkopt ) / nq )
275*
276 IF( left ) THEN
277 IF( notran ) THEN
278 DO i = 1, n, nb
279 len = min( nb, n-i+1 )
280 ldwork = m
281*
282* Multiply bottom part of C by Q12.
283*
284 CALL slacpy( 'All', n1, len, c( n2+1, i ), ldc, work,
285 $ ldwork )
286 CALL strmm( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
287 $ n1, len, one, q( 1, n2+1 ), ldq, work,
288 $ ldwork )
289*
290* Multiply top part of C by Q11.
291*
292 CALL sgemm( 'No Transpose', 'No Transpose', n1, len, n2,
293 $ one, q, ldq, c( 1, i ), ldc, one, work,
294 $ ldwork )
295*
296* Multiply top part of C by Q21.
297*
298 CALL slacpy( 'All', n2, len, c( 1, i ), ldc,
299 $ work( n1+1 ), ldwork )
300 CALL strmm( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
301 $ n2, len, one, q( n1+1, 1 ), ldq,
302 $ work( n1+1 ), ldwork )
303*
304* Multiply bottom part of C by Q22.
305*
306 CALL sgemm( 'No Transpose', 'No Transpose', n2, len, n1,
307 $ one, q( n1+1, n2+1 ), ldq, c( n2+1, i ), ldc,
308 $ one, work( n1+1 ), ldwork )
309*
310* Copy everything back.
311*
312 CALL slacpy( 'All', m, len, work, ldwork, c( 1, i ),
313 $ ldc )
314 END DO
315 ELSE
316 DO i = 1, n, nb
317 len = min( nb, n-i+1 )
318 ldwork = m
319*
320* Multiply bottom part of C by Q21**T.
321*
322 CALL slacpy( 'All', n2, len, c( n1+1, i ), ldc, work,
323 $ ldwork )
324 CALL strmm( 'Left', 'Upper', 'Transpose', 'Non-Unit',
325 $ n2, len, one, q( n1+1, 1 ), ldq, work,
326 $ ldwork )
327*
328* Multiply top part of C by Q11**T.
329*
330 CALL sgemm( 'Transpose', 'No Transpose', n2, len, n1,
331 $ one, q, ldq, c( 1, i ), ldc, one, work,
332 $ ldwork )
333*
334* Multiply top part of C by Q12**T.
335*
336 CALL slacpy( 'All', n1, len, c( 1, i ), ldc,
337 $ work( n2+1 ), ldwork )
338 CALL strmm( 'Left', 'Lower', 'Transpose', 'Non-Unit',
339 $ n1, len, one, q( 1, n2+1 ), ldq,
340 $ work( n2+1 ), ldwork )
341*
342* Multiply bottom part of C by Q22**T.
343*
344 CALL sgemm( 'Transpose', 'No Transpose', n1, len, n2,
345 $ one, q( n1+1, n2+1 ), ldq, c( n1+1, i ), ldc,
346 $ one, work( n2+1 ), ldwork )
347*
348* Copy everything back.
349*
350 CALL slacpy( 'All', m, len, work, ldwork, c( 1, i ),
351 $ ldc )
352 END DO
353 END IF
354 ELSE
355 IF( notran ) THEN
356 DO i = 1, m, nb
357 len = min( nb, m-i+1 )
358 ldwork = len
359*
360* Multiply right part of C by Q21.
361*
362 CALL slacpy( 'All', len, n2, c( i, n1+1 ), ldc, work,
363 $ ldwork )
364 CALL strmm( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
365 $ len, n2, one, q( n1+1, 1 ), ldq, work,
366 $ ldwork )
367*
368* Multiply left part of C by Q11.
369*
370 CALL sgemm( 'No Transpose', 'No Transpose', len, n2, n1,
371 $ one, c( i, 1 ), ldc, q, ldq, one, work,
372 $ ldwork )
373*
374* Multiply left part of C by Q12.
375*
376 CALL slacpy( 'All', len, n1, c( i, 1 ), ldc,
377 $ work( 1 + n2*ldwork ), ldwork )
378 CALL strmm( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
379 $ len, n1, one, q( 1, n2+1 ), ldq,
380 $ work( 1 + n2*ldwork ), ldwork )
381*
382* Multiply right part of C by Q22.
383*
384 CALL sgemm( 'No Transpose', 'No Transpose', len, n1, n2,
385 $ one, c( i, n1+1 ), ldc, q( n1+1, n2+1 ), ldq,
386 $ one, work( 1 + n2*ldwork ), ldwork )
387*
388* Copy everything back.
389*
390 CALL slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
391 $ ldc )
392 END DO
393 ELSE
394 DO i = 1, m, nb
395 len = min( nb, m-i+1 )
396 ldwork = len
397*
398* Multiply right part of C by Q12**T.
399*
400 CALL slacpy( 'All', len, n1, c( i, n2+1 ), ldc, work,
401 $ ldwork )
402 CALL strmm( 'Right', 'Lower', 'Transpose', 'Non-Unit',
403 $ len, n1, one, q( 1, n2+1 ), ldq, work,
404 $ ldwork )
405*
406* Multiply left part of C by Q11**T.
407*
408 CALL sgemm( 'No Transpose', 'Transpose', len, n1, n2,
409 $ one, c( i, 1 ), ldc, q, ldq, one, work,
410 $ ldwork )
411*
412* Multiply left part of C by Q21**T.
413*
414 CALL slacpy( 'All', len, n2, c( i, 1 ), ldc,
415 $ work( 1 + n1*ldwork ), ldwork )
416 CALL strmm( 'Right', 'Upper', 'Transpose', 'Non-Unit',
417 $ len, n2, one, q( n1+1, 1 ), ldq,
418 $ work( 1 + n1*ldwork ), ldwork )
419*
420* Multiply right part of C by Q22**T.
421*
422 CALL sgemm( 'No Transpose', 'Transpose', len, n2, n1,
423 $ one, c( i, n2+1 ), ldc, q( n1+1, n2+1 ), ldq,
424 $ one, work( 1 + n1*ldwork ), ldwork )
425*
426* Copy everything back.
427*
428 CALL slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
429 $ ldc )
430 END DO
431 END IF
432 END IF
433*
434 work( 1 ) = sroundup_lwork( lwkopt )
435 RETURN
436*
437* End of SORM22
438*
439 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
Definition slacpy.f:103
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177
subroutine sorm22(side, trans, m, n, n1, n2, q, ldq, c, ldc, work, lwork, info)
SORM22 multiplies a general matrix by a banded orthogonal matrix.
Definition sorm22.f:163