1*> \brief \b DORCSD2BY1
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DORCSD2BY1 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd2by1.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd2by1.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd2by1.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
22*                              X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
23*                              LDV1T, WORK, LWORK, IWORK, INFO )
24*
25*       .. Scalar Arguments ..
26*       CHARACTER          JOBU1, JOBU2, JOBV1T
27*       INTEGER            INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
28*      $                   M, P, Q
29*       ..
30*       .. Array Arguments ..
31*       DOUBLE PRECISION   THETA(*)
32*       DOUBLE PRECISION   U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
33*      $                   X11(LDX11,*), X21(LDX21,*)
34*       INTEGER            IWORK(*)
35*       ..
36*
37*
38*> \par Purpose:
39*  =============
40*>
41*>\verbatim
42*>
43*> DORCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
44*> orthonormal columns that has been partitioned into a 2-by-1 block
45*> structure:
46*>
47*>                                [  I1 0  0 ]
48*>                                [  0  C  0 ]
49*>          [ X11 ]   [ U1 |    ] [  0  0  0 ]
50*>      X = [-----] = [---------] [----------] V1**T .
51*>          [ X21 ]   [    | U2 ] [  0  0  0 ]
52*>                                [  0  S  0 ]
53*>                                [  0  0  I2]
54*>
55*> X11 is P-by-Q. The orthogonal matrices U1, U2, and V1 are P-by-P,
56*> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
57*> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
58*> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a
59*> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
60*> \endverbatim
61*
62*  Arguments:
63*  ==========
64*
65*> \param[in] JOBU1
66*> \verbatim
67*>          JOBU1 is CHARACTER
68*>          = 'Y':      U1 is computed;
69*>          otherwise:  U1 is not computed.
70*> \endverbatim
71*>
72*> \param[in] JOBU2
73*> \verbatim
74*>          JOBU2 is CHARACTER
75*>          = 'Y':      U2 is computed;
76*>          otherwise:  U2 is not computed.
77*> \endverbatim
78*>
79*> \param[in] JOBV1T
80*> \verbatim
81*>          JOBV1T is CHARACTER
82*>          = 'Y':      V1T is computed;
83*>          otherwise:  V1T is not computed.
84*> \endverbatim
85*>
86*> \param[in] M
87*> \verbatim
88*>          M is INTEGER
89*>          The number of rows in X.
90*> \endverbatim
91*>
92*> \param[in] P
93*> \verbatim
94*>          P is INTEGER
95*>          The number of rows in X11. 0 <= P <= M.
96*> \endverbatim
97*>
98*> \param[in] Q
99*> \verbatim
100*>          Q is INTEGER
101*>          The number of columns in X11 and X21. 0 <= Q <= M.
102*> \endverbatim
103*>
104*> \param[in,out] X11
105*> \verbatim
106*>          X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
107*>          On entry, part of the orthogonal matrix whose CSD is desired.
108*> \endverbatim
109*>
110*> \param[in] LDX11
111*> \verbatim
112*>          LDX11 is INTEGER
113*>          The leading dimension of X11. LDX11 >= MAX(1,P).
114*> \endverbatim
115*>
116*> \param[in,out] X21
117*> \verbatim
118*>          X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
119*>          On entry, part of the orthogonal matrix whose CSD is desired.
120*> \endverbatim
121*>
122*> \param[in] LDX21
123*> \verbatim
124*>          LDX21 is INTEGER
125*>          The leading dimension of X21. LDX21 >= MAX(1,M-P).
126*> \endverbatim
127*>
128*> \param[out] THETA
129*> \verbatim
130*>          THETA is DOUBLE PRECISION array, dimension (R), in which R =
131*>          MIN(P,M-P,Q,M-Q).
132*>          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
133*>          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
134*> \endverbatim
135*>
136*> \param[out] U1
137*> \verbatim
138*>          U1 is DOUBLE PRECISION array, dimension (P)
139*>          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
140*> \endverbatim
141*>
142*> \param[in] LDU1
143*> \verbatim
144*>          LDU1 is INTEGER
145*>          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
146*>          MAX(1,P).
147*> \endverbatim
148*>
149*> \param[out] U2
150*> \verbatim
151*>          U2 is DOUBLE PRECISION array, dimension (M-P)
152*>          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
153*>          matrix U2.
154*> \endverbatim
155*>
156*> \param[in] LDU2
157*> \verbatim
158*>          LDU2 is INTEGER
159*>          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
160*>          MAX(1,M-P).
161*> \endverbatim
162*>
163*> \param[out] V1T
164*> \verbatim
165*>          V1T is DOUBLE PRECISION array, dimension (Q)
166*>          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
167*>          matrix V1**T.
168*> \endverbatim
169*>
170*> \param[in] LDV1T
171*> \verbatim
172*>          LDV1T is INTEGER
173*>          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
174*>          MAX(1,Q).
175*> \endverbatim
176*>
177*> \param[out] WORK
178*> \verbatim
179*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
180*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
181*>          If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
182*>          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
183*>          define the matrix in intermediate bidiagonal-block form
184*>          remaining after nonconvergence. INFO specifies the number
185*>          of nonzero PHI's.
186*> \endverbatim
187*>
188*> \param[in] LWORK
189*> \verbatim
190*>          LWORK is INTEGER
191*>          The dimension of the array WORK.
192*>
193*>          If LWORK = -1, then a workspace query is assumed; the routine
194*>          only calculates the optimal size of the WORK array, returns
195*>          this value as the first entry of the work array, and no error
196*>          message related to LWORK is issued by XERBLA.
197*> \endverbatim
198*>
199*> \param[out] IWORK
200*> \verbatim
201*>          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
202*> \endverbatim
203*>
204*> \param[out] INFO
205*> \verbatim
206*>          INFO is INTEGER
207*>          = 0:  successful exit.
208*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
209*>          > 0:  DBBCSD did not converge. See the description of WORK
210*>                above for details.
211*> \endverbatim
212*
213*> \par References:
214*  ================
215*>
216*>  [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
217*>      Algorithms, 50(1):33-65, 2009.
218*
219*  Authors:
220*  ========
221*
222*> \author Univ. of Tennessee
223*> \author Univ. of California Berkeley
224*> \author Univ. of Colorado Denver
225*> \author NAG Ltd.
226*
227*> \ingroup doubleOTHERcomputational
228*
229*  =====================================================================
230      SUBROUTINE DORCSD2BY1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
231     $                       X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
232     $                       LDV1T, WORK, LWORK, IWORK, INFO )
233*
234*  -- LAPACK computational routine (3.5.0) --
235*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
236*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
237*
238*     .. Scalar Arguments ..
239      CHARACTER          JOBU1, JOBU2, JOBV1T
240      INTEGER            INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
241     $                   M, P, Q
242*     ..
243*     .. Array Arguments ..
244      DOUBLE PRECISION   THETA(*)
245      DOUBLE PRECISION   U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
246     $                   X11(LDX11,*), X21(LDX21,*)
247      INTEGER            IWORK(*)
248*     ..
249*
250*  =====================================================================
251*
252*     .. Parameters ..
253      DOUBLE PRECISION   ONE, ZERO
254      PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
255*     ..
256*     .. Local Scalars ..
257      INTEGER            CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
258     $                   IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
259     $                   IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
260     $                   J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN,
261     $                   LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT,
262     $                   LWORKMIN, LWORKOPT, R
263      LOGICAL            LQUERY, WANTU1, WANTU2, WANTV1T
264*     ..
265*     .. Local Arrays ..
266      DOUBLE PRECISION   DUM1(1), DUM2(1,1)
267*     ..
268*     .. External Subroutines ..
269      EXTERNAL           DBBCSD, DCOPY, DLACPY, DLAPMR, DLAPMT, DORBDB1,
270     $                   DORBDB2, DORBDB3, DORBDB4, DORGLQ, DORGQR,
271     $                   XERBLA
272*     ..
273*     .. External Functions ..
274      LOGICAL            LSAME
275      EXTERNAL           LSAME
276*     ..
277*     .. Intrinsic Function ..
278      INTRINSIC          INT, MAX, MIN
279*     ..
280*     .. Executable Statements ..
281*
282*     Test input arguments
283*
284      INFO = 0
285      WANTU1 = LSAME( JOBU1, 'Y' )
286      WANTU2 = LSAME( JOBU2, 'Y' )
287      WANTV1T = LSAME( JOBV1T, 'Y' )
288      LQUERY = LWORK .EQ. -1
289*
290      IF( M .LT. 0 ) THEN
291         INFO = -4
292      ELSE IF( P .LT. 0 .OR. P .GT. M ) THEN
293         INFO = -5
294      ELSE IF( Q .LT. 0 .OR. Q .GT. M ) THEN
295         INFO = -6
296      ELSE IF( LDX11 .LT. MAX( 1, P ) ) THEN
297         INFO = -8
298      ELSE IF( LDX21 .LT. MAX( 1, M-P ) ) THEN
299         INFO = -10
300      ELSE IF( WANTU1 .AND. LDU1 .LT. MAX( 1, P ) ) THEN
301         INFO = -13
302      ELSE IF( WANTU2 .AND. LDU2 .LT. MAX( 1, M - P ) ) THEN
303         INFO = -15
304      ELSE IF( WANTV1T .AND. LDV1T .LT. MAX( 1, Q ) ) THEN
305         INFO = -17
306      END IF
307*
308      R = MIN( P, M-P, Q, M-Q )
309*
310*     Compute workspace
311*
312*       WORK layout:
313*     |-------------------------------------------------------|
314*     | LWORKOPT (1)                                          |
315*     |-------------------------------------------------------|
316*     | PHI (MAX(1,R-1))                                      |
317*     |-------------------------------------------------------|
318*     | TAUP1 (MAX(1,P))                        | B11D (R)    |
319*     | TAUP2 (MAX(1,M-P))                      | B11E (R-1)  |
320*     | TAUQ1 (MAX(1,Q))                        | B12D (R)    |
321*     |-----------------------------------------| B12E (R-1)  |
322*     | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R)    |
323*     |             |             |             | B21E (R-1)  |
324*     |             |             |             | B22D (R)    |
325*     |             |             |             | B22E (R-1)  |
326*     |             |             |             | DBBCSD WORK |
327*     |-------------------------------------------------------|
328*
329      IF( INFO .EQ. 0 ) THEN
330         IPHI = 2
331         IB11D = IPHI + MAX( 1, R-1 )
332         IB11E = IB11D + MAX( 1, R )
333         IB12D = IB11E + MAX( 1, R - 1 )
334         IB12E = IB12D + MAX( 1, R )
335         IB21D = IB12E + MAX( 1, R - 1 )
336         IB21E = IB21D + MAX( 1, R )
337         IB22D = IB21E + MAX( 1, R - 1 )
338         IB22E = IB22D + MAX( 1, R )
339         IBBCSD = IB22E + MAX( 1, R - 1 )
340         ITAUP1 = IPHI + MAX( 1, R-1 )
341         ITAUP2 = ITAUP1 + MAX( 1, P )
342         ITAUQ1 = ITAUP2 + MAX( 1, M-P )
343         IORBDB = ITAUQ1 + MAX( 1, Q )
344         IORGQR = ITAUQ1 + MAX( 1, Q )
345         IORGLQ = ITAUQ1 + MAX( 1, Q )
346         LORGQRMIN = 1
347         LORGQROPT = 1
348         LORGLQMIN = 1
349         LORGLQOPT = 1
350         IF( R .EQ. Q ) THEN
351            CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
352     $                    DUM1, DUM1, DUM1, DUM1, WORK,
353     $                    -1, CHILDINFO )
354            LORBDB = INT( WORK(1) )
355            IF( WANTU1 .AND. P .GT. 0 ) THEN
356               CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1,
357     $                      CHILDINFO )
358               LORGQRMIN = MAX( LORGQRMIN, P )
359               LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
360            ENDIF
361            IF( WANTU2 .AND. M-P .GT. 0 ) THEN
362               CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1),
363     $                      -1, CHILDINFO )
364               LORGQRMIN = MAX( LORGQRMIN, M-P )
365               LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
366            END IF
367            IF( WANTV1T .AND. Q .GT. 0 ) THEN
368               CALL DORGLQ( Q-1, Q-1, Q-1, V1T, LDV1T,
369     $                      DUM1, WORK(1), -1, CHILDINFO )
370               LORGLQMIN = MAX( LORGLQMIN, Q-1 )
371               LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
372            END IF
373            CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
374     $                   DUM1, U1, LDU1, U2, LDU2, V1T, LDV1T,
375     $                   DUM2, 1, DUM1, DUM1, DUM1,
376     $                   DUM1, DUM1, DUM1, DUM1,
377     $                   DUM1, WORK(1), -1, CHILDINFO )
378            LBBCSD = INT( WORK(1) )
379         ELSE IF( R .EQ. P ) THEN
380            CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
381     $                    DUM1, DUM1, DUM1, DUM1,
382     $                    WORK(1), -1, CHILDINFO )
383            LORBDB = INT( WORK(1) )
384            IF( WANTU1 .AND. P .GT. 0 ) THEN
385               CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, DUM1,
386     $                      WORK(1), -1, CHILDINFO )
387               LORGQRMIN = MAX( LORGQRMIN, P-1 )
388               LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
389            END IF
390            IF( WANTU2 .AND. M-P .GT. 0 ) THEN
391               CALL DORGQR( M-P, M-P, Q, U2, LDU2, DUM1, WORK(1),
392     $                      -1, CHILDINFO )
393               LORGQRMIN = MAX( LORGQRMIN, M-P )
394               LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
395            END IF
396            IF( WANTV1T .AND. Q .GT. 0 ) THEN
397               CALL DORGLQ( Q, Q, R, V1T, LDV1T, DUM1, WORK(1), -1,
398     $                      CHILDINFO )
399               LORGLQMIN = MAX( LORGLQMIN, Q )
400               LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
401            END IF
402            CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
403     $                   DUM1, V1T, LDV1T, DUM2, 1, U1, LDU1,
404     $                   U2, LDU2, DUM1, DUM1, DUM1,
405     $                   DUM1, DUM1, DUM1, DUM1,
406     $                   DUM1, WORK(1), -1, CHILDINFO )
407            LBBCSD = INT( WORK(1) )
408         ELSE IF( R .EQ. M-P ) THEN
409            CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
410     $                    DUM1, DUM1, DUM1, DUM1,
411     $                    WORK(1), -1, CHILDINFO )
412            LORBDB = INT( WORK(1) )
413            IF( WANTU1 .AND. P .GT. 0 ) THEN
414               CALL DORGQR( P, P, Q, U1, LDU1, DUM1, WORK(1), -1,
415     $                      CHILDINFO )
416               LORGQRMIN = MAX( LORGQRMIN, P )
417               LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
418            END IF
419            IF( WANTU2 .AND. M-P .GT. 0 ) THEN
420               CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
421     $                      DUM1, WORK(1), -1, CHILDINFO )
422               LORGQRMIN = MAX( LORGQRMIN, M-P-1 )
423               LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
424            END IF
425            IF( WANTV1T .AND. Q .GT. 0 ) THEN
426               CALL DORGLQ( Q, Q, R, V1T, LDV1T, DUM1, WORK(1), -1,
427     $                      CHILDINFO )
428               LORGLQMIN = MAX( LORGLQMIN, Q )
429               LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
430            END IF
431            CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
432     $                   THETA, DUM1, DUM2, 1, V1T, LDV1T, U2,
433     $                   LDU2, U1, LDU1, DUM1, DUM1, DUM1,
434     $                   DUM1, DUM1, DUM1, DUM1,
435     $                   DUM1, WORK(1), -1, CHILDINFO )
436            LBBCSD = INT( WORK(1) )
437         ELSE
438            CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
439     $                    DUM1, DUM1, DUM1, DUM1,
440     $                    DUM1, WORK(1), -1, CHILDINFO )
441            LORBDB = M + INT( WORK(1) )
442            IF( WANTU1 .AND. P .GT. 0 ) THEN
443               CALL DORGQR( P, P, M-Q, U1, LDU1, DUM1, WORK(1), -1,
444     $                      CHILDINFO )
445               LORGQRMIN = MAX( LORGQRMIN, P )
446               LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
447            END IF
448            IF( WANTU2 .AND. M-P .GT. 0 ) THEN
449               CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, DUM1, WORK(1),
450     $                      -1, CHILDINFO )
451               LORGQRMIN = MAX( LORGQRMIN, M-P )
452               LORGQROPT = MAX( LORGQROPT, INT( WORK(1) ) )
453            END IF
454            IF( WANTV1T .AND. Q .GT. 0 ) THEN
455               CALL DORGLQ( Q, Q, Q, V1T, LDV1T, DUM1, WORK(1), -1,
456     $                      CHILDINFO )
457               LORGLQMIN = MAX( LORGLQMIN, Q )
458               LORGLQOPT = MAX( LORGLQOPT, INT( WORK(1) ) )
459            END IF
460            CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
461     $                   THETA, DUM1, U2, LDU2, U1, LDU1, DUM2,
462     $                   1, V1T, LDV1T, DUM1, DUM1, DUM1,
463     $                   DUM1, DUM1, DUM1, DUM1,
464     $                   DUM1, WORK(1), -1, CHILDINFO )
465            LBBCSD = INT( WORK(1) )
466         END IF
467         LWORKMIN = MAX( IORBDB+LORBDB-1,
468     $                   IORGQR+LORGQRMIN-1,
469     $                   IORGLQ+LORGLQMIN-1,
470     $                   IBBCSD+LBBCSD-1 )
471         LWORKOPT = MAX( IORBDB+LORBDB-1,
472     $                   IORGQR+LORGQROPT-1,
473     $                   IORGLQ+LORGLQOPT-1,
474     $                   IBBCSD+LBBCSD-1 )
475         WORK(1) = LWORKOPT
476         IF( LWORK .LT. LWORKMIN .AND. .NOT.LQUERY ) THEN
477            INFO = -19
478         END IF
479      END IF
480      IF( INFO .NE. 0 ) THEN
481         CALL XERBLA( 'DORCSD2BY1', -INFO )
482         RETURN
483      ELSE IF( LQUERY ) THEN
484         RETURN
485      END IF
486      LORGQR = LWORK-IORGQR+1
487      LORGLQ = LWORK-IORGLQ+1
488*
489*     Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
490*     in which R = MIN(P,M-P,Q,M-Q)
491*
492      IF( R .EQ. Q ) THEN
493*
494*        Case 1: R = Q
495*
496*        Simultaneously bidiagonalize X11 and X21
497*
498         CALL DORBDB1( M, P, Q, X11, LDX11, X21, LDX21, THETA,
499     $                 WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
500     $                 WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
501*
502*        Accumulate Householder reflectors
503*
504         IF( WANTU1 .AND. P .GT. 0 ) THEN
505            CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
506            CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
507     $                   LORGQR, CHILDINFO )
508         END IF
509         IF( WANTU2 .AND. M-P .GT. 0 ) THEN
510            CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
511            CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
512     $                   WORK(IORGQR), LORGQR, CHILDINFO )
513         END IF
514         IF( WANTV1T .AND. Q .GT. 0 ) THEN
515            V1T(1,1) = ONE
516            DO J = 2, Q
517               V1T(1,J) = ZERO
518               V1T(J,1) = ZERO
519            END DO
520            CALL DLACPY( 'U', Q-1, Q-1, X21(1,2), LDX21, V1T(2,2),
521     $                   LDV1T )
522            CALL DORGLQ( Q-1, Q-1, Q-1, V1T(2,2), LDV1T, WORK(ITAUQ1),
523     $                   WORK(IORGLQ), LORGLQ, CHILDINFO )
524         END IF
525*
526*        Simultaneously diagonalize X11 and X21.
527*
528         CALL DBBCSD( JOBU1, JOBU2, JOBV1T, 'N', 'N', M, P, Q, THETA,
529     $                WORK(IPHI), U1, LDU1, U2, LDU2, V1T, LDV1T,
530     $                DUM2, 1, WORK(IB11D), WORK(IB11E),
531     $                WORK(IB12D), WORK(IB12E), WORK(IB21D),
532     $                WORK(IB21E), WORK(IB22D), WORK(IB22E),
533     $                WORK(IBBCSD), LBBCSD, CHILDINFO )
534*
535*        Permute rows and columns to place zero submatrices in
536*        preferred positions
537*
538         IF( Q .GT. 0 .AND. WANTU2 ) THEN
539            DO I = 1, Q
540               IWORK(I) = M - P - Q + I
541            END DO
542            DO I = Q + 1, M - P
543               IWORK(I) = I - Q
544            END DO
545            CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
546         END IF
547      ELSE IF( R .EQ. P ) THEN
548*
549*        Case 2: R = P
550*
551*        Simultaneously bidiagonalize X11 and X21
552*
553         CALL DORBDB2( M, P, Q, X11, LDX11, X21, LDX21, THETA,
554     $                 WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
555     $                 WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
556*
557*        Accumulate Householder reflectors
558*
559         IF( WANTU1 .AND. P .GT. 0 ) THEN
560            U1(1,1) = ONE
561            DO J = 2, P
562               U1(1,J) = ZERO
563               U1(J,1) = ZERO
564            END DO
565            CALL DLACPY( 'L', P-1, P-1, X11(2,1), LDX11, U1(2,2), LDU1 )
566            CALL DORGQR( P-1, P-1, P-1, U1(2,2), LDU1, WORK(ITAUP1),
567     $                   WORK(IORGQR), LORGQR, CHILDINFO )
568         END IF
569         IF( WANTU2 .AND. M-P .GT. 0 ) THEN
570            CALL DLACPY( 'L', M-P, Q, X21, LDX21, U2, LDU2 )
571            CALL DORGQR( M-P, M-P, Q, U2, LDU2, WORK(ITAUP2),
572     $                   WORK(IORGQR), LORGQR, CHILDINFO )
573         END IF
574         IF( WANTV1T .AND. Q .GT. 0 ) THEN
575            CALL DLACPY( 'U', P, Q, X11, LDX11, V1T, LDV1T )
576            CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
577     $                   WORK(IORGLQ), LORGLQ, CHILDINFO )
578         END IF
579*
580*        Simultaneously diagonalize X11 and X21.
581*
582         CALL DBBCSD( JOBV1T, 'N', JOBU1, JOBU2, 'T', M, Q, P, THETA,
583     $                WORK(IPHI), V1T, LDV1T, DUM2, 1, U1, LDU1, U2,
584     $                LDU2, WORK(IB11D), WORK(IB11E), WORK(IB12D),
585     $                WORK(IB12E), WORK(IB21D), WORK(IB21E),
586     $                WORK(IB22D), WORK(IB22E), WORK(IBBCSD), LBBCSD,
587     $                CHILDINFO )
588*
589*        Permute rows and columns to place identity submatrices in
590*        preferred positions
591*
592         IF( Q .GT. 0 .AND. WANTU2 ) THEN
593            DO I = 1, Q
594               IWORK(I) = M - P - Q + I
595            END DO
596            DO I = Q + 1, M - P
597               IWORK(I) = I - Q
598            END DO
599            CALL DLAPMT( .FALSE., M-P, M-P, U2, LDU2, IWORK )
600         END IF
601      ELSE IF( R .EQ. M-P ) THEN
602*
603*        Case 3: R = M-P
604*
605*        Simultaneously bidiagonalize X11 and X21
606*
607         CALL DORBDB3( M, P, Q, X11, LDX11, X21, LDX21, THETA,
608     $                 WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
609     $                 WORK(ITAUQ1), WORK(IORBDB), LORBDB, CHILDINFO )
610*
611*        Accumulate Householder reflectors
612*
613         IF( WANTU1 .AND. P .GT. 0 ) THEN
614            CALL DLACPY( 'L', P, Q, X11, LDX11, U1, LDU1 )
615            CALL DORGQR( P, P, Q, U1, LDU1, WORK(ITAUP1), WORK(IORGQR),
616     $                   LORGQR, CHILDINFO )
617         END IF
618         IF( WANTU2 .AND. M-P .GT. 0 ) THEN
619            U2(1,1) = ONE
620            DO J = 2, M-P
621               U2(1,J) = ZERO
622               U2(J,1) = ZERO
623            END DO
624            CALL DLACPY( 'L', M-P-1, M-P-1, X21(2,1), LDX21, U2(2,2),
625     $                   LDU2 )
626            CALL DORGQR( M-P-1, M-P-1, M-P-1, U2(2,2), LDU2,
627     $                   WORK(ITAUP2), WORK(IORGQR), LORGQR, CHILDINFO )
628         END IF
629         IF( WANTV1T .AND. Q .GT. 0 ) THEN
630            CALL DLACPY( 'U', M-P, Q, X21, LDX21, V1T, LDV1T )
631            CALL DORGLQ( Q, Q, R, V1T, LDV1T, WORK(ITAUQ1),
632     $                   WORK(IORGLQ), LORGLQ, CHILDINFO )
633         END IF
634*
635*        Simultaneously diagonalize X11 and X21.
636*
637         CALL DBBCSD( 'N', JOBV1T, JOBU2, JOBU1, 'T', M, M-Q, M-P,
638     $                THETA, WORK(IPHI), DUM2, 1, V1T, LDV1T, U2,
639     $                LDU2, U1, LDU1, WORK(IB11D), WORK(IB11E),
640     $                WORK(IB12D), WORK(IB12E), WORK(IB21D),
641     $                WORK(IB21E), WORK(IB22D), WORK(IB22E),
642     $                WORK(IBBCSD), LBBCSD, CHILDINFO )
643*
644*        Permute rows and columns to place identity submatrices in
645*        preferred positions
646*
647         IF( Q .GT. R ) THEN
648            DO I = 1, R
649               IWORK(I) = Q - R + I
650            END DO
651            DO I = R + 1, Q
652               IWORK(I) = I - R
653            END DO
654            IF( WANTU1 ) THEN
655               CALL DLAPMT( .FALSE., P, Q, U1, LDU1, IWORK )
656            END IF
657            IF( WANTV1T ) THEN
658               CALL DLAPMR( .FALSE., Q, Q, V1T, LDV1T, IWORK )
659            END IF
660         END IF
661      ELSE
662*
663*        Case 4: R = M-Q
664*
665*        Simultaneously bidiagonalize X11 and X21
666*
667         CALL DORBDB4( M, P, Q, X11, LDX11, X21, LDX21, THETA,
668     $                 WORK(IPHI), WORK(ITAUP1), WORK(ITAUP2),
669     $                 WORK(ITAUQ1), WORK(IORBDB), WORK(IORBDB+M),
670     $                 LORBDB-M, CHILDINFO )
671*
672*        Accumulate Householder reflectors
673*
674         IF( WANTU2 .AND. M-P .GT. 0 ) THEN
675            CALL DCOPY( M-P, WORK(IORBDB+P), 1, U2, 1 )
676         END IF
677         IF( WANTU1 .AND. P .GT. 0 ) THEN
678            CALL DCOPY( P, WORK(IORBDB), 1, U1, 1 )
679            DO J = 2, P
680               U1(1,J) = ZERO
681            END DO
682            CALL DLACPY( 'L', P-1, M-Q-1, X11(2,1), LDX11, U1(2,2),
683     $                   LDU1 )
684            CALL DORGQR( P, P, M-Q, U1, LDU1, WORK(ITAUP1),
685     $                   WORK(IORGQR), LORGQR, CHILDINFO )
686         END IF
687         IF( WANTU2 .AND. M-P .GT. 0 ) THEN
688            DO J = 2, M-P
689               U2(1,J) = ZERO
690            END DO
691            CALL DLACPY( 'L', M-P-1, M-Q-1, X21(2,1), LDX21, U2(2,2),
692     $                   LDU2 )
693            CALL DORGQR( M-P, M-P, M-Q, U2, LDU2, WORK(ITAUP2),
694     $                   WORK(IORGQR), LORGQR, CHILDINFO )
695         END IF
696         IF( WANTV1T .AND. Q .GT. 0 ) THEN
697            CALL DLACPY( 'U', M-Q, Q, X21, LDX21, V1T, LDV1T )
698            CALL DLACPY( 'U', P-(M-Q), Q-(M-Q), X11(M-Q+1,M-Q+1), LDX11,
699     $                   V1T(M-Q+1,M-Q+1), LDV1T )
700            CALL DLACPY( 'U', -P+Q, Q-P, X21(M-Q+1,P+1), LDX21,
701     $                   V1T(P+1,P+1), LDV1T )
702            CALL DORGLQ( Q, Q, Q, V1T, LDV1T, WORK(ITAUQ1),
703     $                   WORK(IORGLQ), LORGLQ, CHILDINFO )
704         END IF
705*
706*        Simultaneously diagonalize X11 and X21.
707*
708         CALL DBBCSD( JOBU2, JOBU1, 'N', JOBV1T, 'N', M, M-P, M-Q,
709     $                THETA, WORK(IPHI), U2, LDU2, U1, LDU1, DUM2,
710     $                1, V1T, LDV1T, WORK(IB11D), WORK(IB11E),
711     $                WORK(IB12D), WORK(IB12E), WORK(IB21D),
712     $                WORK(IB21E), WORK(IB22D), WORK(IB22E),
713     $                WORK(IBBCSD), LBBCSD, CHILDINFO )
714*
715*        Permute rows and columns to place identity submatrices in
716*        preferred positions
717*
718         IF( P .GT. R ) THEN
719            DO I = 1, R
720               IWORK(I) = P - R + I
721            END DO
722            DO I = R + 1, P
723               IWORK(I) = I - R
724            END DO
725            IF( WANTU1 ) THEN
726               CALL DLAPMT( .FALSE., P, P, U1, LDU1, IWORK )
727            END IF
728            IF( WANTV1T ) THEN
729               CALL DLAPMR( .FALSE., P, Q, V1T, LDV1T, IWORK )
730            END IF
731         END IF
732      END IF
733*
734      RETURN
735*
736*     End of DORCSD2BY1
737*
738      END
739
740