LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ sorm22()

subroutine sorm22 ( character  SIDE,
character  TRANS,
integer  M,
integer  N,
integer  N1,
integer  N2,
real, dimension( ldq, * )  Q,
integer  LDQ,
real, dimension( ldc, * )  C,
integer  LDC,
real, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

SORM22 multiplies a general matrix by a banded orthogonal matrix.

Download SORM22 + dependencies [TGZ] [ZIP] [TXT]

Purpose
  SORM22 overwrites the general real M-by-N matrix C with

                  SIDE = 'L'     SIDE = 'R'
  TRANS = 'N':      Q * C          C * Q
  TRANS = 'T':      Q**T * C       C * Q**T

  where Q is a real orthogonal matrix of order NQ, with NQ = M if
  SIDE = 'L' and NQ = N if SIDE = 'R'.
  The orthogonal matrix Q processes a 2-by-2 block structure

         [  Q11  Q12  ]
     Q = [            ]
         [  Q21  Q22  ],

  where Q12 is an N1-by-N1 lower triangular matrix and Q21 is an
  N2-by-N2 upper triangular matrix.
Parameters
[in]SIDE
          SIDE is CHARACTER*1
          = 'L': apply Q or Q**T from the Left;
          = 'R': apply Q or Q**T from the Right.
[in]TRANS
          TRANS is CHARACTER*1
          = 'N':  apply Q (No transpose);
          = 'C':  apply Q**T (Conjugate transpose).
[in]M
          M is INTEGER
          The number of rows of the matrix C. M >= 0.
[in]N
          N is INTEGER
          The number of columns of the matrix C. N >= 0.
[in]N1
[in]N2
          N1 is INTEGER
          N2 is INTEGER
          The dimension of Q12 and Q21, respectively. N1, N2 >= 0.
          The following requirement must be satisfied:
          N1 + N2 = M if SIDE = 'L' and N1 + N2 = N if SIDE = 'R'.
[in]Q
          Q is REAL array, dimension
                              (LDQ,M) if SIDE = 'L'
                              (LDQ,N) if SIDE = 'R'
[in]LDQ
          LDQ is INTEGER
          The leading dimension of the array Q.
          LDQ >= max(1,M) if SIDE = 'L'; LDQ >= max(1,N) if SIDE = 'R'.
[in,out]C
          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N matrix C.
          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M).
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.
          If SIDE = 'L', LWORK >= max(1,N);
          if SIDE = 'R', LWORK >= max(1,M).
          For optimum performance LWORK >= M*N.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 161 of file sorm22.f.

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*
slacpy
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
strmm
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:177
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
Here is the call graph for this function:
Here is the caller graph for this function: