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

◆ dorm22()

subroutine dorm22 ( character  side,
character  trans,
integer  m,
integer  n,
integer  n1,
integer  n2,
double precision, dimension( ldq, * )  q,
integer  ldq,
double precision, dimension( ldc, * )  c,
integer  ldc,
double precision, dimension( * )  work,
integer  lwork,
integer  info 
)

DORM22 multiplies a general matrix by a banded orthogonal matrix.

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

Purpose
  DORM22 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 dorm22.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 DOUBLE PRECISION Q( LDQ, * ), C( LDC, * ), WORK( * )
176* ..
177*
178* =====================================================================
179*
180* .. Parameters ..
181 DOUBLE PRECISION ONE
182 parameter( one = 1.0d+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 dgemm, dlacpy, dtrmm, xerbla
194* ..
195* .. Intrinsic Functions ..
196 INTRINSIC dble, 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 ) = dble( lwkopt )
241 END IF
242*
243 IF( info.NE.0 ) THEN
244 CALL xerbla( 'DORM22', -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 DTRMM.
258*
259 IF( n1.EQ.0 ) THEN
260 CALL dtrmm( 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 dtrmm( 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 dlacpy( 'All', n1, len, c( n2+1, i ), ldc, work,
284 $ ldwork )
285 CALL dtrmm( '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 dgemm( '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 dlacpy( 'All', n2, len, c( 1, i ), ldc,
298 $ work( n1+1 ), ldwork )
299 CALL dtrmm( '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 dgemm( '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 dlacpy( '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 dlacpy( 'All', n2, len, c( n1+1, i ), ldc, work,
322 $ ldwork )
323 CALL dtrmm( '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 dgemm( '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 dlacpy( 'All', n1, len, c( 1, i ), ldc,
336 $ work( n2+1 ), ldwork )
337 CALL dtrmm( '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 dgemm( '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 dlacpy( '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 dlacpy( 'All', len, n2, c( i, n1+1 ), ldc, work,
362 $ ldwork )
363 CALL dtrmm( '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 dgemm( '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 dlacpy( 'All', len, n1, c( i, 1 ), ldc,
376 $ work( 1 + n2*ldwork ), ldwork )
377 CALL dtrmm( '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 dgemm( '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 dlacpy( '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 dlacpy( 'All', len, n1, c( i, n2+1 ), ldc, work,
400 $ ldwork )
401 CALL dtrmm( '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 dgemm( '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 dlacpy( 'All', len, n2, c( i, 1 ), ldc,
414 $ work( 1 + n1*ldwork ), ldwork )
415 CALL dtrmm( '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 dgemm( '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 dlacpy( 'All', len, n, work, ldwork, c( i, 1 ),
428 $ ldc )
429 END DO
430 END IF
431 END IF
432*
433 work( 1 ) = dble( lwkopt )
434 RETURN
435*
436* End of DORM22
437*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:103
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dtrmm
subroutine dtrmm(side, uplo, transa, diag, m, n, alpha, a, lda, b, ldb)
DTRMM
Definition dtrmm.f:177
Here is the call graph for this function:
Here is the caller graph for this function: