1*> \brief \b DORM22 multiplies a general matrix by a banded orthogonal matrix.
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DORM22 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorm22.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorm22.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorm22.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*     SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
22*    $                   WORK, LWORK, INFO )
23*
24*     .. Scalar Arguments ..
25*     CHARACTER          SIDE, TRANS
26*     INTEGER            M, N, N1, N2, LDQ, LDC, LWORK, INFO
27*     ..
28*     .. Array Arguments ..
29*     DOUBLE PRECISION   Q( LDQ, * ), C( LDC, * ), WORK( * )
30*     ..
31*
32*> \par Purpose
33*  ============
34*>
35*> \verbatim
36*>
37*>
38*>  DORM22 overwrites the general real M-by-N matrix C with
39*>
40*>                  SIDE = 'L'     SIDE = 'R'
41*>  TRANS = 'N':      Q * C          C * Q
42*>  TRANS = 'T':      Q**T * C       C * Q**T
43*>
44*>  where Q is a real orthogonal matrix of order NQ, with NQ = M if
45*>  SIDE = 'L' and NQ = N if SIDE = 'R'.
46*>  The orthogonal matrix Q processes a 2-by-2 block structure
47*>
48*>         [  Q11  Q12  ]
49*>     Q = [            ]
50*>         [  Q21  Q22  ],
51*>
52*>  where Q12 is an N1-by-N1 lower triangular matrix and Q21 is an
53*>  N2-by-N2 upper triangular matrix.
54*> \endverbatim
55*
56*  Arguments:
57*  ==========
58*
59*> \param[in] SIDE
60*> \verbatim
61*>          SIDE is CHARACTER*1
62*>          = 'L': apply Q or Q**T from the Left;
63*>          = 'R': apply Q or Q**T from the Right.
64*> \endverbatim
65*>
66*> \param[in] TRANS
67*> \verbatim
68*>          TRANS is CHARACTER*1
69*>          = 'N':  apply Q (No transpose);
70*>          = 'C':  apply Q**T (Conjugate transpose).
71*> \endverbatim
72*>
73*> \param[in] M
74*> \verbatim
75*>          M is INTEGER
76*>          The number of rows of the matrix C. M >= 0.
77*> \endverbatim
78*>
79*> \param[in] N
80*> \verbatim
81*>          N is INTEGER
82*>          The number of columns of the matrix C. N >= 0.
83*> \endverbatim
84*>
85*> \param[in] N1
86*> \param[in] N2
87*> \verbatim
88*>          N1 is INTEGER
89*>          N2 is INTEGER
90*>          The dimension of Q12 and Q21, respectively. N1, N2 >= 0.
91*>          The following requirement must be satisfied:
92*>          N1 + N2 = M if SIDE = 'L' and N1 + N2 = N if SIDE = 'R'.
93*> \endverbatim
94*>
95*> \param[in] Q
96*> \verbatim
97*>          Q is DOUBLE PRECISION array, dimension
98*>                                       (LDQ,M) if SIDE = 'L'
99*>                                       (LDQ,N) if SIDE = 'R'
100*> \endverbatim
101*>
102*> \param[in] LDQ
103*> \verbatim
104*>          LDQ is INTEGER
105*>          The leading dimension of the array Q.
106*>          LDQ >= max(1,M) if SIDE = 'L'; LDQ >= max(1,N) if SIDE = 'R'.
107*> \endverbatim
108*>
109*> \param[in,out] C
110*> \verbatim
111*>          C is DOUBLE PRECISION array, dimension (LDC,N)
112*>          On entry, the M-by-N matrix C.
113*>          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
114*> \endverbatim
115*>
116*> \param[in] LDC
117*> \verbatim
118*>          LDC is INTEGER
119*>          The leading dimension of the array C. LDC >= max(1,M).
120*> \endverbatim
121*>
122*> \param[out] WORK
123*> \verbatim
124*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
125*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
126*> \endverbatim
127*>
128*> \param[in] LWORK
129*> \verbatim
130*>          LWORK is INTEGER
131*>          The dimension of the array WORK.
132*>          If SIDE = 'L', LWORK >= max(1,N);
133*>          if SIDE = 'R', LWORK >= max(1,M).
134*>          For optimum performance LWORK >= M*N.
135*>
136*>          If LWORK = -1, then a workspace query is assumed; the routine
137*>          only calculates the optimal size of the WORK array, returns
138*>          this value as the first entry of the WORK array, and no error
139*>          message related to LWORK is issued by XERBLA.
140*> \endverbatim
141*>
142*> \param[out] INFO
143*> \verbatim
144*>          INFO is INTEGER
145*>          = 0:  successful exit
146*>          < 0:  if INFO = -i, the i-th argument had an illegal value
147*> \endverbatim
148*
149*
150*  Authors:
151*  ========
152*
153*> \author Univ. of Tennessee
154*> \author Univ. of California Berkeley
155*> \author Univ. of Colorado Denver
156*> \author NAG Ltd.
157*
158*> \ingroup complexOTHERcomputational
159*
160*  =====================================================================
161      SUBROUTINE DORM22( SIDE, TRANS, M, N, N1, N2, Q, LDQ, C, LDC,
162     $                   WORK, LWORK, INFO )
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*
438      END
439