LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ 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 *
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
Here is the call graph for this function:
Here is the caller graph for this function: