LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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.
Date
January 2015

Definition at line 165 of file sorm22.f.

165 *
166 * -- LAPACK computational routine (version 3.6.0) --
167 * -- LAPACK is a software package provided by Univ. of Tennessee, --
168 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
169 * January 2015
170 *
171  IMPLICIT NONE
172 *
173 * .. Scalar Arguments ..
174  CHARACTER side, trans
175  INTEGER m, n, n1, n2, ldq, ldc, lwork, info
176 * ..
177 * .. Array Arguments ..
178  REAL q( ldq, * ), c( ldc, * ), work( * )
179 * ..
180 *
181 * =====================================================================
182 *
183 * .. Parameters ..
184  REAL one
185  parameter ( one = 1.0e+0 )
186 *
187 * .. Local Scalars ..
188  LOGICAL left, lquery, notran
189  INTEGER i, ldwork, len, lwkopt, nb, nq, nw
190 * ..
191 * .. External Functions ..
192  LOGICAL lsame
193  EXTERNAL lsame
194 * ..
195 * .. External Subroutines ..
196  EXTERNAL sgemm, slacpy, strmm, xerbla
197 * ..
198 * .. Intrinsic Functions ..
199  INTRINSIC REAL, max, min
200 * ..
201 * .. Executable Statements ..
202 *
203 * Test the input arguments
204 *
205  info = 0
206  left = lsame( side, 'L' )
207  notran = lsame( trans, 'N' )
208  lquery = ( lwork.EQ.-1 )
209 *
210 * NQ is the order of Q;
211 * NW is the minimum dimension of WORK.
212 *
213  IF( left ) THEN
214  nq = m
215  ELSE
216  nq = n
217  END IF
218  nw = nq
219  IF( n1.EQ.0 .OR. n2.EQ.0 ) nw = 1
220  IF( .NOT.left .AND. .NOT.lsame( side, 'R' ) ) THEN
221  info = -1
222  ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'T' ) )
223  $ THEN
224  info = -2
225  ELSE IF( m.LT.0 ) THEN
226  info = -3
227  ELSE IF( n.LT.0 ) THEN
228  info = -4
229  ELSE IF( n1.LT.0 .OR. n1+n2.NE.nq ) THEN
230  info = -5
231  ELSE IF( n2.LT.0 ) THEN
232  info = -6
233  ELSE IF( ldq.LT.max( 1, nq ) ) THEN
234  info = -8
235  ELSE IF( ldc.LT.max( 1, m ) ) THEN
236  info = -10
237  ELSE IF( lwork.LT.nw .AND. .NOT.lquery ) THEN
238  info = -12
239  END IF
240 *
241  IF( info.EQ.0 ) THEN
242  lwkopt = m*n
243  work( 1 ) = REAL( lwkopt )
244  END IF
245 *
246  IF( info.NE.0 ) THEN
247  CALL xerbla( 'SORM22', -info )
248  RETURN
249  ELSE IF( lquery ) THEN
250  RETURN
251  END IF
252 *
253 * Quick return if possible
254 *
255  IF( m.EQ.0 .OR. n.EQ.0 ) THEN
256  work( 1 ) = 1
257  RETURN
258  END IF
259 *
260 * Degenerate cases (N1 = 0 or N2 = 0) are handled using STRMM.
261 *
262  IF( n1.EQ.0 ) THEN
263  CALL strmm( side, 'Upper', trans, 'Non-Unit', m, n, one,
264  $ q, ldq, c, ldc )
265  work( 1 ) = one
266  RETURN
267  ELSE IF( n2.EQ.0 ) THEN
268  CALL strmm( side, 'Lower', trans, 'Non-Unit', m, n, one,
269  $ q, ldq, c, ldc )
270  work( 1 ) = one
271  RETURN
272  END IF
273 *
274 * Compute the largest chunk size available from the workspace.
275 *
276  nb = max( 1, min( lwork, lwkopt ) / nq )
277 *
278  IF( left ) THEN
279  IF( notran ) THEN
280  DO i = 1, n, nb
281  len = min( nb, n-i+1 )
282  ldwork = m
283 *
284 * Multiply bottom part of C by Q12.
285 *
286  CALL slacpy( 'All', n1, len, c( n2+1, i ), ldc, work,
287  $ ldwork )
288  CALL strmm( 'Left', 'Lower', 'No Transpose', 'Non-Unit',
289  $ n1, len, one, q( 1, n2+1 ), ldq, work,
290  $ ldwork )
291 *
292 * Multiply top part of C by Q11.
293 *
294  CALL sgemm( 'No Transpose', 'No Transpose', n1, len, n2,
295  $ one, q, ldq, c( 1, i ), ldc, one, work,
296  $ ldwork )
297 *
298 * Multiply top part of C by Q21.
299 *
300  CALL slacpy( 'All', n2, len, c( 1, i ), ldc,
301  $ work( n1+1 ), ldwork )
302  CALL strmm( 'Left', 'Upper', 'No Transpose', '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, n1,
309  $ one, q( n1+1, n2+1 ), ldq, c( n2+1, i ), ldc,
310  $ one, work( n1+1 ), ldwork )
311 *
312 * Copy everything back.
313 *
314  CALL slacpy( 'All', m, len, work, ldwork, c( 1, i ),
315  $ ldc )
316  END DO
317  ELSE
318  DO i = 1, n, nb
319  len = min( nb, n-i+1 )
320  ldwork = m
321 *
322 * Multiply bottom part of C by Q21**T.
323 *
324  CALL slacpy( 'All', n2, len, c( n1+1, i ), ldc, work,
325  $ ldwork )
326  CALL strmm( 'Left', 'Upper', 'Transpose', 'Non-Unit',
327  $ n2, len, one, q( n1+1, 1 ), ldq, work,
328  $ ldwork )
329 *
330 * Multiply top part of C by Q11**T.
331 *
332  CALL sgemm( 'Transpose', 'No Transpose', n2, len, n1,
333  $ one, q, ldq, c( 1, i ), ldc, one, work,
334  $ ldwork )
335 *
336 * Multiply top part of C by Q12**T.
337 *
338  CALL slacpy( 'All', n1, len, c( 1, i ), ldc,
339  $ work( n2+1 ), ldwork )
340  CALL strmm( 'Left', 'Lower', 'Transpose', 'Non-Unit',
341  $ n1, len, one, q( 1, n2+1 ), ldq,
342  $ work( n2+1 ), ldwork )
343 *
344 * Multiply bottom part of C by Q22**T.
345 *
346  CALL sgemm( 'Transpose', 'No Transpose', n1, len, n2,
347  $ one, q( n1+1, n2+1 ), ldq, c( n1+1, i ), ldc,
348  $ one, work( n2+1 ), ldwork )
349 *
350 * Copy everything back.
351 *
352  CALL slacpy( 'All', m, len, work, ldwork, c( 1, i ),
353  $ ldc )
354  END DO
355  END IF
356  ELSE
357  IF( notran ) THEN
358  DO i = 1, m, nb
359  len = min( nb, m-i+1 )
360  ldwork = len
361 *
362 * Multiply right part of C by Q21.
363 *
364  CALL slacpy( 'All', len, n2, c( i, n1+1 ), ldc, work,
365  $ ldwork )
366  CALL strmm( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
367  $ len, n2, one, q( n1+1, 1 ), ldq, work,
368  $ ldwork )
369 *
370 * Multiply left part of C by Q11.
371 *
372  CALL sgemm( 'No Transpose', 'No Transpose', len, n2, n1,
373  $ one, c( i, 1 ), ldc, q, ldq, one, work,
374  $ ldwork )
375 *
376 * Multiply left part of C by Q12.
377 *
378  CALL slacpy( 'All', len, n1, c( i, 1 ), ldc,
379  $ work( 1 + n2*ldwork ), ldwork )
380  CALL strmm( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
381  $ len, n1, one, q( 1, n2+1 ), ldq,
382  $ work( 1 + n2*ldwork ), ldwork )
383 *
384 * Multiply right part of C by Q22.
385 *
386  CALL sgemm( 'No Transpose', 'No Transpose', len, n1, n2,
387  $ one, c( i, n1+1 ), ldc, q( n1+1, n2+1 ), ldq,
388  $ one, work( 1 + n2*ldwork ), ldwork )
389 *
390 * Copy everything back.
391 *
392  CALL slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
393  $ ldc )
394  END DO
395  ELSE
396  DO i = 1, m, nb
397  len = min( nb, m-i+1 )
398  ldwork = len
399 *
400 * Multiply right part of C by Q12**T.
401 *
402  CALL slacpy( 'All', len, n1, c( i, n2+1 ), ldc, work,
403  $ ldwork )
404  CALL strmm( 'Right', 'Lower', 'Transpose', 'Non-Unit',
405  $ len, n1, one, q( 1, n2+1 ), ldq, work,
406  $ ldwork )
407 *
408 * Multiply left part of C by Q11**T.
409 *
410  CALL sgemm( 'No Transpose', 'Transpose', len, n1, n2,
411  $ one, c( i, 1 ), ldc, q, ldq, one, work,
412  $ ldwork )
413 *
414 * Multiply left part of C by Q21**T.
415 *
416  CALL slacpy( 'All', len, n2, c( i, 1 ), ldc,
417  $ work( 1 + n1*ldwork ), ldwork )
418  CALL strmm( 'Right', 'Upper', 'Transpose', 'Non-Unit',
419  $ len, n2, one, q( n1+1, 1 ), ldq,
420  $ work( 1 + n1*ldwork ), ldwork )
421 *
422 * Multiply right part of C by Q22**T.
423 *
424  CALL sgemm( 'No Transpose', 'Transpose', len, n2, n1,
425  $ one, c( i, n2+1 ), ldc, q( n1+1, n2+1 ), ldq,
426  $ one, work( 1 + n1*ldwork ), ldwork )
427 *
428 * Copy everything back.
429 *
430  CALL slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
431  $ ldc )
432  END DO
433  END IF
434  END IF
435 *
436  work( 1 ) = REAL( lwkopt )
437  RETURN
438 *
439 * End of SORM22
440 *
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine strmm(SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA, B, LDB)
STRMM
Definition: strmm.f:179
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: