1*> \brief \b CUNBDB4
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download CUNBDB4 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cunbdb4.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cunbdb4.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cunbdb4.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
22*                           TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK,
23*                           INFO )
24*
25*       .. Scalar Arguments ..
26*       INTEGER            INFO, LWORK, M, P, Q, LDX11, LDX21
27*       ..
28*       .. Array Arguments ..
29*       REAL               PHI(*), THETA(*)
30*       COMPLEX            PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
31*      $                   WORK(*), X11(LDX11,*), X21(LDX21,*)
32*       ..
33*
34*
35*> \par Purpose:
36*  =============
37*>
38*>\verbatim
39*>
40*> CUNBDB4 simultaneously bidiagonalizes the blocks of a tall and skinny
41*> matrix X with orthonomal columns:
42*>
43*>                            [ B11 ]
44*>      [ X11 ]   [ P1 |    ] [  0  ]
45*>      [-----] = [---------] [-----] Q1**T .
46*>      [ X21 ]   [    | P2 ] [ B21 ]
47*>                            [  0  ]
48*>
49*> X11 is P-by-Q, and X21 is (M-P)-by-Q. M-Q must be no larger than P,
50*> M-P, or Q. Routines CUNBDB1, CUNBDB2, and CUNBDB3 handle cases in
51*> which M-Q is not the minimum dimension.
52*>
53*> The unitary matrices P1, P2, and Q1 are P-by-P, (M-P)-by-(M-P),
54*> and (M-Q)-by-(M-Q), respectively. They are represented implicitly by
55*> Householder vectors.
56*>
57*> B11 and B12 are (M-Q)-by-(M-Q) bidiagonal matrices represented
58*> implicitly by angles THETA, PHI.
59*>
60*>\endverbatim
61*
62*  Arguments:
63*  ==========
64*
65*> \param[in] M
66*> \verbatim
67*>          M is INTEGER
68*>           The number of rows X11 plus the number of rows in X21.
69*> \endverbatim
70*>
71*> \param[in] P
72*> \verbatim
73*>          P is INTEGER
74*>           The number of rows in X11. 0 <= P <= M.
75*> \endverbatim
76*>
77*> \param[in] Q
78*> \verbatim
79*>          Q is INTEGER
80*>           The number of columns in X11 and X21. 0 <= Q <= M and
81*>           M-Q <= min(P,M-P,Q).
82*> \endverbatim
83*>
84*> \param[in,out] X11
85*> \verbatim
86*>          X11 is COMPLEX array, dimension (LDX11,Q)
87*>           On entry, the top block of the matrix X to be reduced. On
88*>           exit, the columns of tril(X11) specify reflectors for P1 and
89*>           the rows of triu(X11,1) specify reflectors for Q1.
90*> \endverbatim
91*>
92*> \param[in] LDX11
93*> \verbatim
94*>          LDX11 is INTEGER
95*>           The leading dimension of X11. LDX11 >= P.
96*> \endverbatim
97*>
98*> \param[in,out] X21
99*> \verbatim
100*>          X21 is COMPLEX array, dimension (LDX21,Q)
101*>           On entry, the bottom block of the matrix X to be reduced. On
102*>           exit, the columns of tril(X21) specify reflectors for P2.
103*> \endverbatim
104*>
105*> \param[in] LDX21
106*> \verbatim
107*>          LDX21 is INTEGER
108*>           The leading dimension of X21. LDX21 >= M-P.
109*> \endverbatim
110*>
111*> \param[out] THETA
112*> \verbatim
113*>          THETA is REAL array, dimension (Q)
114*>           The entries of the bidiagonal blocks B11, B21 are defined by
115*>           THETA and PHI. See Further Details.
116*> \endverbatim
117*>
118*> \param[out] PHI
119*> \verbatim
120*>          PHI is REAL array, dimension (Q-1)
121*>           The entries of the bidiagonal blocks B11, B21 are defined by
122*>           THETA and PHI. See Further Details.
123*> \endverbatim
124*>
125*> \param[out] TAUP1
126*> \verbatim
127*>          TAUP1 is COMPLEX array, dimension (P)
128*>           The scalar factors of the elementary reflectors that define
129*>           P1.
130*> \endverbatim
131*>
132*> \param[out] TAUP2
133*> \verbatim
134*>          TAUP2 is COMPLEX array, dimension (M-P)
135*>           The scalar factors of the elementary reflectors that define
136*>           P2.
137*> \endverbatim
138*>
139*> \param[out] TAUQ1
140*> \verbatim
141*>          TAUQ1 is COMPLEX array, dimension (Q)
142*>           The scalar factors of the elementary reflectors that define
143*>           Q1.
144*> \endverbatim
145*>
146*> \param[out] PHANTOM
147*> \verbatim
148*>          PHANTOM is COMPLEX array, dimension (M)
149*>           The routine computes an M-by-1 column vector Y that is
150*>           orthogonal to the columns of [ X11; X21 ]. PHANTOM(1:P) and
151*>           PHANTOM(P+1:M) contain Householder vectors for Y(1:P) and
152*>           Y(P+1:M), respectively.
153*> \endverbatim
154*>
155*> \param[out] WORK
156*> \verbatim
157*>          WORK is COMPLEX array, dimension (LWORK)
158*> \endverbatim
159*>
160*> \param[in] LWORK
161*> \verbatim
162*>          LWORK is INTEGER
163*>           The dimension of the array WORK. LWORK >= M-Q.
164*>
165*>           If LWORK = -1, then a workspace query is assumed; the routine
166*>           only calculates the optimal size of the WORK array, returns
167*>           this value as the first entry of the WORK array, and no error
168*>           message related to LWORK is issued by XERBLA.
169*> \endverbatim
170*>
171*> \param[out] INFO
172*> \verbatim
173*>          INFO is INTEGER
174*>           = 0:  successful exit.
175*>           < 0:  if INFO = -i, the i-th argument had an illegal value.
176*> \endverbatim
177*
178*  Authors:
179*  ========
180*
181*> \author Univ. of Tennessee
182*> \author Univ. of California Berkeley
183*> \author Univ. of Colorado Denver
184*> \author NAG Ltd.
185*
186*> \ingroup complexOTHERcomputational
187*
188*> \par Further Details:
189*  =====================
190
191*> \verbatim
192*>
193*>  The upper-bidiagonal blocks B11, B21 are represented implicitly by
194*>  angles THETA(1), ..., THETA(Q) and PHI(1), ..., PHI(Q-1). Every entry
195*>  in each bidiagonal band is a product of a sine or cosine of a THETA
196*>  with a sine or cosine of a PHI. See [1] or CUNCSD for details.
197*>
198*>  P1, P2, and Q1 are represented as products of elementary reflectors.
199*>  See CUNCSD2BY1 for details on generating P1, P2, and Q1 using CUNGQR
200*>  and CUNGLQ.
201*> \endverbatim
202*
203*> \par References:
204*  ================
205*>
206*>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
207*>      Algorithms, 50(1):33-65, 2009.
208*>
209*  =====================================================================
210      SUBROUTINE CUNBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI,
211     $                    TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK,
212     $                    INFO )
213*
214*  -- LAPACK computational routine --
215*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
216*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217*
218*     .. Scalar Arguments ..
219      INTEGER            INFO, LWORK, M, P, Q, LDX11, LDX21
220*     ..
221*     .. Array Arguments ..
222      REAL               PHI(*), THETA(*)
223      COMPLEX            PHANTOM(*), TAUP1(*), TAUP2(*), TAUQ1(*),
224     $                   WORK(*), X11(LDX11,*), X21(LDX21,*)
225*     ..
226*
227*  ====================================================================
228*
229*     .. Parameters ..
230      COMPLEX            NEGONE, ONE, ZERO
231      PARAMETER          ( NEGONE = (-1.0E0,0.0E0), ONE = (1.0E0,0.0E0),
232     $                     ZERO = (0.0E0,0.0E0) )
233*     ..
234*     .. Local Scalars ..
235      REAL               C, S
236      INTEGER            CHILDINFO, I, ILARF, IORBDB5, J, LLARF,
237     $                   LORBDB5, LWORKMIN, LWORKOPT
238      LOGICAL            LQUERY
239*     ..
240*     .. External Subroutines ..
241      EXTERNAL           CLARF, CLARFGP, CUNBDB5, CSROT, CSCAL, CLACGV,
242     $                   XERBLA
243*     ..
244*     .. External Functions ..
245      REAL               SCNRM2
246      EXTERNAL           SCNRM2
247*     ..
248*     .. Intrinsic Function ..
249      INTRINSIC          ATAN2, COS, MAX, SIN, SQRT
250*     ..
251*     .. Executable Statements ..
252*
253*     Test input arguments
254*
255      INFO = 0
256      LQUERY = LWORK .EQ. -1
257*
258      IF( M .LT. 0 ) THEN
259         INFO = -1
260      ELSE IF( P .LT. M-Q .OR. M-P .LT. M-Q ) THEN
261         INFO = -2
262      ELSE IF( Q .LT. M-Q .OR. Q .GT. M ) THEN
263         INFO = -3
264      ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN
265         INFO = -5
266      ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN
267         INFO = -7
268      END IF
269*
270*     Compute workspace
271*
272      IF( INFO .EQ. 0 ) THEN
273         ILARF = 2
274         LLARF = MAX( Q-1, P-1, M-P-1 )
275         IORBDB5 = 2
276         LORBDB5 = Q
277         LWORKOPT = ILARF + LLARF - 1
278         LWORKOPT = MAX( LWORKOPT, IORBDB5 + LORBDB5 - 1 )
279         LWORKMIN = LWORKOPT
280         WORK(1) = LWORKOPT
281         IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN
282           INFO = -14
283         END IF
284      END IF
285      IF( INFO .NE. 0 ) THEN
286         CALL XERBLA( 'CUNBDB4', -INFO )
287         RETURN
288      ELSE IF( LQUERY ) THEN
289         RETURN
290      END IF
291*
292*     Reduce columns 1, ..., M-Q of X11 and X21
293*
294      DO I = 1, M-Q
295*
296         IF( I .EQ. 1 ) THEN
297            DO J = 1, M
298               PHANTOM(J) = ZERO
299            END DO
300            CALL CUNBDB5( P, M-P, Q, PHANTOM(1), 1, PHANTOM(P+1), 1,
301     $                    X11, LDX11, X21, LDX21, WORK(IORBDB5),
302     $                    LORBDB5, CHILDINFO )
303            CALL CSCAL( P, NEGONE, PHANTOM(1), 1 )
304            CALL CLARFGP( P, PHANTOM(1), PHANTOM(2), 1, TAUP1(1) )
305            CALL CLARFGP( M-P, PHANTOM(P+1), PHANTOM(P+2), 1, TAUP2(1) )
306            THETA(I) = ATAN2( REAL( PHANTOM(1) ), REAL( PHANTOM(P+1) ) )
307            C = COS( THETA(I) )
308            S = SIN( THETA(I) )
309            PHANTOM(1) = ONE
310            PHANTOM(P+1) = ONE
311            CALL CLARF( 'L', P, Q, PHANTOM(1), 1, CONJG(TAUP1(1)), X11,
312     $                  LDX11, WORK(ILARF) )
313            CALL CLARF( 'L', M-P, Q, PHANTOM(P+1), 1, CONJG(TAUP2(1)),
314     $                  X21, LDX21, WORK(ILARF) )
315         ELSE
316            CALL CUNBDB5( P-I+1, M-P-I+1, Q-I+1, X11(I,I-1), 1,
317     $                    X21(I,I-1), 1, X11(I,I), LDX11, X21(I,I),
318     $                    LDX21, WORK(IORBDB5), LORBDB5, CHILDINFO )
319            CALL CSCAL( P-I+1, NEGONE, X11(I,I-1), 1 )
320            CALL CLARFGP( P-I+1, X11(I,I-1), X11(I+1,I-1), 1, TAUP1(I) )
321            CALL CLARFGP( M-P-I+1, X21(I,I-1), X21(I+1,I-1), 1,
322     $                    TAUP2(I) )
323            THETA(I) = ATAN2( REAL( X11(I,I-1) ), REAL( X21(I,I-1) ) )
324            C = COS( THETA(I) )
325            S = SIN( THETA(I) )
326            X11(I,I-1) = ONE
327            X21(I,I-1) = ONE
328            CALL CLARF( 'L', P-I+1, Q-I+1, X11(I,I-1), 1,
329     $                  CONJG(TAUP1(I)), X11(I,I), LDX11, WORK(ILARF) )
330            CALL CLARF( 'L', M-P-I+1, Q-I+1, X21(I,I-1), 1,
331     $                  CONJG(TAUP2(I)), X21(I,I), LDX21, WORK(ILARF) )
332         END IF
333*
334         CALL CSROT( Q-I+1, X11(I,I), LDX11, X21(I,I), LDX21, S, -C )
335         CALL CLACGV( Q-I+1, X21(I,I), LDX21 )
336         CALL CLARFGP( Q-I+1, X21(I,I), X21(I,I+1), LDX21, TAUQ1(I) )
337         C = REAL( X21(I,I) )
338         X21(I,I) = ONE
339         CALL CLARF( 'R', P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I),
340     $               X11(I+1,I), LDX11, WORK(ILARF) )
341         CALL CLARF( 'R', M-P-I, Q-I+1, X21(I,I), LDX21, TAUQ1(I),
342     $               X21(I+1,I), LDX21, WORK(ILARF) )
343         CALL CLACGV( Q-I+1, X21(I,I), LDX21 )
344         IF( I .LT. M-Q ) THEN
345            S = SQRT( SCNRM2( P-I, X11(I+1,I), 1 )**2
346     $              + SCNRM2( M-P-I, X21(I+1,I), 1 )**2 )
347            PHI(I) = ATAN2( S, C )
348         END IF
349*
350      END DO
351*
352*     Reduce the bottom-right portion of X11 to [ I 0 ]
353*
354      DO I = M - Q + 1, P
355         CALL CLACGV( Q-I+1, X11(I,I), LDX11 )
356         CALL CLARFGP( Q-I+1, X11(I,I), X11(I,I+1), LDX11, TAUQ1(I) )
357         X11(I,I) = ONE
358         CALL CLARF( 'R', P-I, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
359     $               X11(I+1,I), LDX11, WORK(ILARF) )
360         CALL CLARF( 'R', Q-P, Q-I+1, X11(I,I), LDX11, TAUQ1(I),
361     $               X21(M-Q+1,I), LDX21, WORK(ILARF) )
362         CALL CLACGV( Q-I+1, X11(I,I), LDX11 )
363      END DO
364*
365*     Reduce the bottom-right portion of X21 to [ 0 I ]
366*
367      DO I = P + 1, Q
368         CALL CLACGV( Q-I+1, X21(M-Q+I-P,I), LDX21 )
369         CALL CLARFGP( Q-I+1, X21(M-Q+I-P,I), X21(M-Q+I-P,I+1), LDX21,
370     $                 TAUQ1(I) )
371         X21(M-Q+I-P,I) = ONE
372         CALL CLARF( 'R', Q-I, Q-I+1, X21(M-Q+I-P,I), LDX21, TAUQ1(I),
373     $               X21(M-Q+I-P+1,I), LDX21, WORK(ILARF) )
374         CALL CLACGV( Q-I+1, X21(M-Q+I-P,I), LDX21 )
375      END DO
376*
377      RETURN
378*
379*     End of CUNBDB4
380*
381      END
382
383