LAPACK 3.12.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 REAL SROUNDUP_LWORK
191 EXTERNAL lsame, sroundup_lwork
192* ..
193* .. External Subroutines ..
194 EXTERNAL sgemm, slacpy, strmm, xerbla
195* ..
196* .. Intrinsic Functions ..
197 INTRINSIC max, min
198* ..
199* .. Executable Statements ..
200*
201* Test the input arguments
202*
203 info = 0
204 left = lsame( side, 'L' )
205 notran = lsame( trans, 'N' )
206 lquery = ( lwork.EQ.-1 )
207*
208* NQ is the order of Q;
209* NW is the minimum dimension of WORK.
210*
211 IF( left ) THEN
212 nq = m
213 ELSE
214 nq = n
215 END IF
216 nw = nq
217 IF( n1.EQ.0 .OR. n2.EQ.0 ) nw = 1
218 IF( .NOT.left .AND. .NOT.lsame( side, 'R' ) ) THEN
219 info = -1
220 ELSE IF( .NOT.lsame( trans, 'N' ) .AND. .NOT.lsame( trans, 'T' ) )
221 $ THEN
222 info = -2
223 ELSE IF( m.LT.0 ) THEN
224 info = -3
225 ELSE IF( n.LT.0 ) THEN
226 info = -4
227 ELSE IF( n1.LT.0 .OR. n1+n2.NE.nq ) THEN
228 info = -5
229 ELSE IF( n2.LT.0 ) THEN
230 info = -6
231 ELSE IF( ldq.LT.max( 1, nq ) ) THEN
232 info = -8
233 ELSE IF( ldc.LT.max( 1, m ) ) THEN
234 info = -10
235 ELSE IF( lwork.LT.nw .AND. .NOT.lquery ) THEN
236 info = -12
237 END IF
238*
239 IF( info.EQ.0 ) THEN
240 lwkopt = m*n
241 work( 1 ) = sroundup_lwork( lwkopt )
242 END IF
243*
244 IF( info.NE.0 ) THEN
245 CALL xerbla( 'SORM22', -info )
246 RETURN
247 ELSE IF( lquery ) THEN
248 RETURN
249 END IF
250*
251* Quick return if possible
252*
253 IF( m.EQ.0 .OR. n.EQ.0 ) THEN
254 work( 1 ) = 1
255 RETURN
256 END IF
257*
258* Degenerate cases (N1 = 0 or N2 = 0) are handled using STRMM.
259*
260 IF( n1.EQ.0 ) THEN
261 CALL strmm( side, 'Upper', trans, 'Non-Unit', m, n, one,
262 $ q, ldq, c, ldc )
263 work( 1 ) = one
264 RETURN
265 ELSE IF( n2.EQ.0 ) THEN
266 CALL strmm( side, 'Lower', trans, 'Non-Unit', m, n, one,
267 $ q, ldq, c, ldc )
268 work( 1 ) = one
269 RETURN
270 END IF
271*
272* Compute the largest chunk size available from the workspace.
273*
274 nb = max( 1, min( lwork, lwkopt ) / nq )
275*
276 IF( left ) THEN
277 IF( notran ) THEN
278 DO i = 1, n, nb
279 len = min( nb, n-i+1 )
280 ldwork = m
281*
282* Multiply bottom part of C by Q12.
283*
284 CALL slacpy( 'All', n1, len, c( n2+1, i ), ldc, work,
285 $ ldwork )
286 CALL strmm( 'Left', 'Lower', 'No Transpose', '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, n2,
293 $ one, q, ldq, c( 1, i ), ldc, one, work,
294 $ ldwork )
295*
296* Multiply top part of C by Q21.
297*
298 CALL slacpy( 'All', n2, len, c( 1, i ), ldc,
299 $ work( n1+1 ), ldwork )
300 CALL strmm( 'Left', 'Upper', 'No Transpose', 'Non-Unit',
301 $ n2, len, one, q( n1+1, 1 ), ldq,
302 $ work( n1+1 ), ldwork )
303*
304* Multiply bottom part of C by Q22.
305*
306 CALL sgemm( 'No Transpose', 'No Transpose', n2, len, n1,
307 $ one, q( n1+1, n2+1 ), ldq, c( n2+1, i ), ldc,
308 $ one, work( n1+1 ), ldwork )
309*
310* Copy everything back.
311*
312 CALL slacpy( 'All', m, len, work, ldwork, c( 1, i ),
313 $ ldc )
314 END DO
315 ELSE
316 DO i = 1, n, nb
317 len = min( nb, n-i+1 )
318 ldwork = m
319*
320* Multiply bottom part of C by Q21**T.
321*
322 CALL slacpy( 'All', n2, len, c( n1+1, i ), ldc, work,
323 $ ldwork )
324 CALL strmm( 'Left', 'Upper', 'Transpose', 'Non-Unit',
325 $ n2, len, one, q( n1+1, 1 ), ldq, work,
326 $ ldwork )
327*
328* Multiply top part of C by Q11**T.
329*
330 CALL sgemm( 'Transpose', 'No Transpose', n2, len, n1,
331 $ one, q, ldq, c( 1, i ), ldc, one, work,
332 $ ldwork )
333*
334* Multiply top part of C by Q12**T.
335*
336 CALL slacpy( 'All', n1, len, c( 1, i ), ldc,
337 $ work( n2+1 ), ldwork )
338 CALL strmm( 'Left', 'Lower', 'Transpose', 'Non-Unit',
339 $ n1, len, one, q( 1, n2+1 ), ldq,
340 $ work( n2+1 ), ldwork )
341*
342* Multiply bottom part of C by Q22**T.
343*
344 CALL sgemm( 'Transpose', 'No Transpose', n1, len, n2,
345 $ one, q( n1+1, n2+1 ), ldq, c( n1+1, i ), ldc,
346 $ one, work( n2+1 ), ldwork )
347*
348* Copy everything back.
349*
350 CALL slacpy( 'All', m, len, work, ldwork, c( 1, i ),
351 $ ldc )
352 END DO
353 END IF
354 ELSE
355 IF( notran ) THEN
356 DO i = 1, m, nb
357 len = min( nb, m-i+1 )
358 ldwork = len
359*
360* Multiply right part of C by Q21.
361*
362 CALL slacpy( 'All', len, n2, c( i, n1+1 ), ldc, work,
363 $ ldwork )
364 CALL strmm( 'Right', 'Upper', 'No Transpose', 'Non-Unit',
365 $ len, n2, one, q( n1+1, 1 ), ldq, work,
366 $ ldwork )
367*
368* Multiply left part of C by Q11.
369*
370 CALL sgemm( 'No Transpose', 'No Transpose', len, n2, n1,
371 $ one, c( i, 1 ), ldc, q, ldq, one, work,
372 $ ldwork )
373*
374* Multiply left part of C by Q12.
375*
376 CALL slacpy( 'All', len, n1, c( i, 1 ), ldc,
377 $ work( 1 + n2*ldwork ), ldwork )
378 CALL strmm( 'Right', 'Lower', 'No Transpose', 'Non-Unit',
379 $ len, n1, one, q( 1, n2+1 ), ldq,
380 $ work( 1 + n2*ldwork ), ldwork )
381*
382* Multiply right part of C by Q22.
383*
384 CALL sgemm( 'No Transpose', 'No Transpose', len, n1, n2,
385 $ one, c( i, n1+1 ), ldc, q( n1+1, n2+1 ), ldq,
386 $ one, work( 1 + n2*ldwork ), ldwork )
387*
388* Copy everything back.
389*
390 CALL slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
391 $ ldc )
392 END DO
393 ELSE
394 DO i = 1, m, nb
395 len = min( nb, m-i+1 )
396 ldwork = len
397*
398* Multiply right part of C by Q12**T.
399*
400 CALL slacpy( 'All', len, n1, c( i, n2+1 ), ldc, work,
401 $ ldwork )
402 CALL strmm( 'Right', 'Lower', 'Transpose', 'Non-Unit',
403 $ len, n1, one, q( 1, n2+1 ), ldq, work,
404 $ ldwork )
405*
406* Multiply left part of C by Q11**T.
407*
408 CALL sgemm( 'No Transpose', 'Transpose', len, n1, n2,
409 $ one, c( i, 1 ), ldc, q, ldq, one, work,
410 $ ldwork )
411*
412* Multiply left part of C by Q21**T.
413*
414 CALL slacpy( 'All', len, n2, c( i, 1 ), ldc,
415 $ work( 1 + n1*ldwork ), ldwork )
416 CALL strmm( 'Right', 'Upper', 'Transpose', 'Non-Unit',
417 $ len, n2, one, q( n1+1, 1 ), ldq,
418 $ work( 1 + n1*ldwork ), ldwork )
419*
420* Multiply right part of C by Q22**T.
421*
422 CALL sgemm( 'No Transpose', 'Transpose', len, n2, n1,
423 $ one, c( i, n2+1 ), ldc, q( n1+1, n2+1 ), ldq,
424 $ one, work( 1 + n1*ldwork ), ldwork )
425*
426* Copy everything back.
427*
428 CALL slacpy( 'All', len, n, work, ldwork, c( i, 1 ),
429 $ ldc )
430 END DO
431 END IF
432 END IF
433*
434 work( 1 ) = sroundup_lwork( lwkopt )
435 RETURN
436*
437* End of SORM22
438*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
sgemm
subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
SGEMM
Definition sgemm.f:188
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
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
strmm
subroutine strmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
STRMM
Definition strmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: