LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 complexOTHERcomputational
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  EXTERNAL lsame
191 * ..
192 * .. External Subroutines ..
193  EXTERNAL sgemm, slacpy, strmm, xerbla
194 * ..
195 * .. Intrinsic Functions ..
196  INTRINSIC real, 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 ) = real( 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', '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 sgemm( '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 slacpy( 'All', n2, len, c( 1, i ), ldc,
298  $ work( n1+1 ), ldwork )
299  CALL strmm( '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 sgemm( '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 slacpy( '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 slacpy( 'All', n2, len, c( n1+1, i ), ldc, work,
322  $ ldwork )
323  CALL strmm( '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 sgemm( '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 slacpy( 'All', n1, len, c( 1, i ), ldc,
336  $ work( n2+1 ), ldwork )
337  CALL strmm( '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 sgemm( '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 slacpy( '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 slacpy( 'All', len, n2, c( i, n1+1 ), ldc, work,
362  $ ldwork )
363  CALL strmm( '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 sgemm( '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 slacpy( 'All', len, n1, c( i, 1 ), ldc,
376  $ work( 1 + n2*ldwork ), ldwork )
377  CALL strmm( '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 sgemm( '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 slacpy( '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 slacpy( 'All', len, n1, c( i, n2+1 ), ldc, work,
400  $ ldwork )
401  CALL strmm( '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 sgemm( '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 slacpy( 'All', len, n2, c( i, 1 ), ldc,
414  $ work( 1 + n1*ldwork ), ldwork )
415  CALL strmm( '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 sgemm( '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 slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
428  $ ldc )
429  END DO
430  END IF
431  END IF
432 *
433  work( 1 ) = real( lwkopt )
434  RETURN
435 *
436 * End of SORM22
437 *
438  END
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 xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:177
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187