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