LAPACK 3.12.1
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*> Download SORM22 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorm22.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorm22.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorm22.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE SORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
20* $ WORK, LWORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER SIDE, TRANS
24* INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
25* ..
26* .. Array Arguments ..
27* REAL Q( LDQ, * ), C( LDC, * ), WORK( * )
28* ..
29*
30*> \par Purpose
31* ============
32*>
33*> \verbatim
34*>
35*>
36*> SORM22 overwrites the general real M-by-N matrix C with
37*>
38*> SIDE = 'L' SIDE = 'R'
39*> TRANS = 'N': Q * C C * Q
40*> TRANS = 'T': Q**T * C C * Q**T
41*>
42*> where Q is a real orthogonal matrix of order NQ, with NQ = M if
43*> SIDE = 'L' and NQ = N if SIDE = 'R'.
44*> The orthogonal matrix Q processes a 2-by-2 block structure
45*>
46*> [ Q11 Q12 ]
47*> Q = [ ]
48*> [ Q21 Q22 ],
49*>
50*> where Q12 is an N1-by-N1 lower triangular matrix and Q21 is an
51*> N2-by-N2 upper triangular matrix.
52*> \endverbatim
53*
54* Arguments:
55* ==========
56*
57*> \param[in] SIDE
58*> \verbatim
59*> SIDE is CHARACTER*1
60*> = 'L': apply Q or Q**T from the Left;
61*> = 'R': apply Q or Q**T from the Right.
62*> \endverbatim
63*>
64*> \param[in] TRANS
65*> \verbatim
66*> TRANS is CHARACTER*1
67*> = 'N': apply Q (No transpose);
68*> = 'C': apply Q**T (Conjugate transpose).
69*> \endverbatim
70*>
71*> \param[in] M
72*> \verbatim
73*> M is INTEGER
74*> The number of rows of the matrix C. M >= 0.
75*> \endverbatim
76*>
77*> \param[in] N
78*> \verbatim
79*> N is INTEGER
80*> The number of columns of the matrix C. N >= 0.
81*> \endverbatim
82*>
83*> \param[in] N1
84*> \param[in] N2
85*> \verbatim
86*> N1 is INTEGER
87*> N2 is INTEGER
88*> The dimension of Q12 and Q21, respectively. N1, N2 >= 0.
89*> The following requirement must be satisfied:
90*> N1 + N2 = M if SIDE = 'L' and N1 + N2 = N if SIDE = 'R'.
91*> \endverbatim
92*>
93*> \param[in] Q
94*> \verbatim
95*> Q is REAL array, dimension
96*> (LDQ,M) if SIDE = 'L'
97*> (LDQ,N) if SIDE = 'R'
98*> \endverbatim
99*>
100*> \param[in] LDQ
101*> \verbatim
102*> LDQ is INTEGER
103*> The leading dimension of the array Q.
104*> LDQ >= max(1,M) if SIDE = 'L'; LDQ >= max(1,N) if SIDE = 'R'.
105*> \endverbatim
106*>
107*> \param[in,out] C
108*> \verbatim
109*> C is REAL array, dimension (LDC,N)
110*> On entry, the M-by-N matrix C.
111*> On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
112*> \endverbatim
113*>
114*> \param[in] LDC
115*> \verbatim
116*> LDC is INTEGER
117*> The leading dimension of the array C. LDC >= max(1,M).
118*> \endverbatim
119*>
120*> \param[out] WORK
121*> \verbatim
122*> WORK is REAL array, dimension (MAX(1,LWORK))
123*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
124*> \endverbatim
125*>
126*> \param[in] LWORK
127*> \verbatim
128*> LWORK is INTEGER
129*> The dimension of the array WORK.
130*> If SIDE = 'L', LWORK >= max(1,N);
131*> if SIDE = 'R', LWORK >= max(1,M).
132*> For optimum performance LWORK >= M*N.
133*>
134*> If LWORK = -1, then a workspace query is assumed; the routine
135*> only calculates the optimal size of the WORK array, returns
136*> this value as the first entry of the WORK array, and no error
137*> message related to LWORK is issued by XERBLA.
138*> \endverbatim
139*>
140*> \param[out] INFO
141*> \verbatim
142*> INFO is INTEGER
143*> = 0: successful exit
144*> < 0: if INFO = -i, the i-th argument had an illegal value
145*> \endverbatim
146*
147*
148* Authors:
149* ========
150*
151*> \author Univ. of Tennessee
152*> \author Univ. of California Berkeley
153*> \author Univ. of Colorado Denver
154*> \author NAG Ltd.
155*
156*> \ingroup unm22
157*
158* =====================================================================
159 SUBROUTINE sorm22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
160 $ WORK, LWORK, INFO )
161*
162* -- LAPACK computational routine --
163* -- LAPACK is a software package provided by Univ. of Tennessee, --
164* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
165*
166 IMPLICIT NONE
167*
168* .. Scalar Arguments ..
169 CHARACTER SIDE, TRANS
170 INTEGER M, N, N1, N2, LDQ, LDC, LWORK, INFO
171* ..
172* .. Array Arguments ..
173 REAL Q( LDQ, * ), C( LDC, * ), WORK( * )
174* ..
175*
176* =====================================================================
177*
178* .. Parameters ..
179 REAL ONE
180 parameter( one = 1.0e+0 )
181*
182* .. Local Scalars ..
183 LOGICAL LEFT, LQUERY, NOTRAN
184 INTEGER I, LDWORK, LEN, LWKOPT, NB, NQ, NW
185* ..
186* .. External Functions ..
187 LOGICAL LSAME
188 REAL SROUNDUP_LWORK
189 EXTERNAL lsame, sroundup_lwork
190* ..
191* .. External Subroutines ..
192 EXTERNAL sgemm, slacpy, strmm, xerbla
193* ..
194* .. Intrinsic Functions ..
195 INTRINSIC max, min
196* ..
197* .. Executable Statements ..
198*
199* Test the input arguments
200*
201 info = 0
202 left = lsame( side, 'L' )
203 notran = lsame( trans, 'N' )
204 lquery = ( lwork.EQ.-1 )
205*
206* NQ is the order of Q;
207* NW is the minimum dimension of WORK.
208*
209 IF( left ) THEN
210 nq = m
211 ELSE
212 nq = n
213 END IF
214 nw = nq
215 IF( n1.EQ.0 .OR. n2.EQ.0 ) nw = 1
216 IF( .NOT.left .AND. .NOT.lsame( side, 'R' ) ) THEN
217 info = -1
218 ELSE IF( .NOT.lsame( trans, 'N' ) .AND.
219 $ .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 ) = sroundup_lwork( lwkopt )
241 END IF
242*
243 IF( info.NE.0 ) THEN
244 CALL xerbla( 'SORM22', -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 STRMM.
258*
259 IF( n1.EQ.0 ) THEN
260 CALL strmm( 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 strmm( 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 slacpy( 'All', n1, len, c( n2+1, i ), ldc, work,
284 $ ldwork )
285 CALL strmm( 'Left', 'Lower', 'No Transpose',
286 $ '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,
293 $ n2,
294 $ one, q, ldq, c( 1, i ), ldc, one, work,
295 $ ldwork )
296*
297* Multiply top part of C by Q21.
298*
299 CALL slacpy( 'All', n2, len, c( 1, i ), ldc,
300 $ work( n1+1 ), ldwork )
301 CALL strmm( 'Left', 'Upper', 'No Transpose',
302 $ '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 sgemm( 'No Transpose', 'No Transpose', n2, len,
309 $ n1,
310 $ one, q( n1+1, n2+1 ), ldq, c( n2+1, i ), ldc,
311 $ one, work( n1+1 ), ldwork )
312*
313* Copy everything back.
314*
315 CALL slacpy( 'All', m, len, work, ldwork, c( 1, i ),
316 $ ldc )
317 END DO
318 ELSE
319 DO i = 1, n, nb
320 len = min( nb, n-i+1 )
321 ldwork = m
322*
323* Multiply bottom part of C by Q21**T.
324*
325 CALL slacpy( 'All', n2, len, c( n1+1, i ), ldc, work,
326 $ ldwork )
327 CALL strmm( 'Left', 'Upper', 'Transpose', 'Non-Unit',
328 $ n2, len, one, q( n1+1, 1 ), ldq, work,
329 $ ldwork )
330*
331* Multiply top part of C by Q11**T.
332*
333 CALL sgemm( 'Transpose', 'No Transpose', n2, len, n1,
334 $ one, q, ldq, c( 1, i ), ldc, one, work,
335 $ ldwork )
336*
337* Multiply top part of C by Q12**T.
338*
339 CALL slacpy( 'All', n1, len, c( 1, i ), ldc,
340 $ work( n2+1 ), ldwork )
341 CALL strmm( 'Left', 'Lower', 'Transpose', 'Non-Unit',
342 $ n1, len, one, q( 1, n2+1 ), ldq,
343 $ work( n2+1 ), ldwork )
344*
345* Multiply bottom part of C by Q22**T.
346*
347 CALL sgemm( 'Transpose', 'No Transpose', n1, len, n2,
348 $ one, q( n1+1, n2+1 ), ldq, c( n1+1, i ), ldc,
349 $ one, work( n2+1 ), ldwork )
350*
351* Copy everything back.
352*
353 CALL slacpy( 'All', m, len, work, ldwork, c( 1, i ),
354 $ ldc )
355 END DO
356 END IF
357 ELSE
358 IF( notran ) THEN
359 DO i = 1, m, nb
360 len = min( nb, m-i+1 )
361 ldwork = len
362*
363* Multiply right part of C by Q21.
364*
365 CALL slacpy( 'All', len, n2, c( i, n1+1 ), ldc, work,
366 $ ldwork )
367 CALL strmm( 'Right', 'Upper', 'No Transpose',
368 $ 'Non-Unit',
369 $ len, n2, one, q( n1+1, 1 ), ldq, work,
370 $ ldwork )
371*
372* Multiply left part of C by Q11.
373*
374 CALL sgemm( 'No Transpose', 'No Transpose', len, n2,
375 $ n1,
376 $ one, c( i, 1 ), ldc, q, ldq, one, work,
377 $ ldwork )
378*
379* Multiply left part of C by Q12.
380*
381 CALL slacpy( 'All', len, n1, c( i, 1 ), ldc,
382 $ work( 1 + n2*ldwork ), ldwork )
383 CALL strmm( 'Right', 'Lower', 'No Transpose',
384 $ 'Non-Unit',
385 $ len, n1, one, q( 1, n2+1 ), ldq,
386 $ work( 1 + n2*ldwork ), ldwork )
387*
388* Multiply right part of C by Q22.
389*
390 CALL sgemm( 'No Transpose', 'No Transpose', len, n1,
391 $ n2,
392 $ one, c( i, n1+1 ), ldc, q( n1+1, n2+1 ), ldq,
393 $ one, work( 1 + n2*ldwork ), ldwork )
394*
395* Copy everything back.
396*
397 CALL slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
398 $ ldc )
399 END DO
400 ELSE
401 DO i = 1, m, nb
402 len = min( nb, m-i+1 )
403 ldwork = len
404*
405* Multiply right part of C by Q12**T.
406*
407 CALL slacpy( 'All', len, n1, c( i, n2+1 ), ldc, work,
408 $ ldwork )
409 CALL strmm( 'Right', 'Lower', 'Transpose', 'Non-Unit',
410 $ len, n1, one, q( 1, n2+1 ), ldq, work,
411 $ ldwork )
412*
413* Multiply left part of C by Q11**T.
414*
415 CALL sgemm( 'No Transpose', 'Transpose', len, n1, n2,
416 $ one, c( i, 1 ), ldc, q, ldq, one, work,
417 $ ldwork )
418*
419* Multiply left part of C by Q21**T.
420*
421 CALL slacpy( 'All', len, n2, c( i, 1 ), ldc,
422 $ work( 1 + n1*ldwork ), ldwork )
423 CALL strmm( 'Right', 'Upper', 'Transpose', 'Non-Unit',
424 $ len, n2, one, q( n1+1, 1 ), ldq,
425 $ work( 1 + n1*ldwork ), ldwork )
426*
427* Multiply right part of C by Q22**T.
428*
429 CALL sgemm( 'No Transpose', 'Transpose', len, n2, n1,
430 $ one, c( i, n2+1 ), ldc, q( n1+1, n2+1 ), ldq,
431 $ one, work( 1 + n1*ldwork ), ldwork )
432*
433* Copy everything back.
434*
435 CALL slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
436 $ ldc )
437 END DO
438 END IF
439 END IF
440*
441 work( 1 ) = sroundup_lwork( lwkopt )
442 RETURN
443*
444* End of SORM22
445*
446 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:101
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:161