LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
sorcsd2by1.f
Go to the documentation of this file.
1 *> \brief \b SORCSD2BY1
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SORCSD2BY1 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorcsd2by1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorcsd2by1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorcsd2by1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SORCSD2BY1( 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 * REAL THETA(*)
32 * REAL 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 *> SORCSD2BY1 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 *> [ I 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 I ]
54 *>
55 *> X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P,
56 *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
57 *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
58 *> which R = MIN(P,M-P,Q,M-Q).
59 *>
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 and columns in X.
90 *> \endverbatim
91 *>
92 *> \param[in] P
93 *> \verbatim
94 *> P is INTEGER
95 *> The number of rows in X11 and X12. 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 REAL array, dimension (LDX11,Q)
107 *> On entry, part of the orthogonal matrix whose CSD is
108 *> desired.
109 *> \endverbatim
110 *>
111 *> \param[in] LDX11
112 *> \verbatim
113 *> LDX11 is INTEGER
114 *> The leading dimension of X11. LDX11 >= MAX(1,P).
115 *> \endverbatim
116 *>
117 *> \param[in,out] X21
118 *> \verbatim
119 *> X21 is REAL array, dimension (LDX21,Q)
120 *> On entry, part of the orthogonal matrix whose CSD is
121 *> desired.
122 *> \endverbatim
123 *>
124 *> \param[in] LDX21
125 *> \verbatim
126 *> LDX21 is INTEGER
127 *> The leading dimension of X21. LDX21 >= MAX(1,M-P).
128 *> \endverbatim
129 *>
130 *> \param[out] THETA
131 *> \verbatim
132 *> THETA is REAL array, dimension (R), in which R =
133 *> MIN(P,M-P,Q,M-Q).
134 *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
135 *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
136 *> \endverbatim
137 *>
138 *> \param[out] U1
139 *> \verbatim
140 *> U1 is REAL array, dimension (P)
141 *> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
142 *> \endverbatim
143 *>
144 *> \param[in] LDU1
145 *> \verbatim
146 *> LDU1 is INTEGER
147 *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
148 *> MAX(1,P).
149 *> \endverbatim
150 *>
151 *> \param[out] U2
152 *> \verbatim
153 *> U2 is REAL array, dimension (M-P)
154 *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
155 *> matrix U2.
156 *> \endverbatim
157 *>
158 *> \param[in] LDU2
159 *> \verbatim
160 *> LDU2 is INTEGER
161 *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
162 *> MAX(1,M-P).
163 *> \endverbatim
164 *>
165 *> \param[out] V1T
166 *> \verbatim
167 *> V1T is REAL array, dimension (Q)
168 *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
169 *> matrix V1**T.
170 *> \endverbatim
171 *>
172 *> \param[in] LDV1T
173 *> \verbatim
174 *> LDV1T is INTEGER
175 *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
176 *> MAX(1,Q).
177 *> \endverbatim
178 *>
179 *> \param[out] WORK
180 *> \verbatim
181 *> WORK is REAL array, dimension (MAX(1,LWORK))
182 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
183 *> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
184 *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
185 *> define the matrix in intermediate bidiagonal-block form
186 *> remaining after nonconvergence. INFO specifies the number
187 *> of nonzero PHI's.
188 *> \endverbatim
189 *>
190 *> \param[in] LWORK
191 *> \verbatim
192 *> LWORK is INTEGER
193 *> The dimension of the array WORK.
194 *> \endverbatim
195 *>
196 *> If LWORK = -1, then a workspace query is assumed; the routine
197 *> only calculates the optimal size of the WORK array, returns
198 *> this value as the first entry of the work array, and no error
199 *> message related to LWORK is issued by XERBLA.
200 *> \param[out] IWORK
201 *> \verbatim
202 *> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
203 *> \endverbatim
204 *>
205 *> \param[out] INFO
206 *> \verbatim
207 *> INFO is INTEGER
208 *> = 0: successful exit.
209 *> < 0: if INFO = -i, the i-th argument had an illegal value.
210 *> > 0: SBBCSD did not converge. See the description of WORK
211 *> above for details.
212 *> \endverbatim
213 *>
214 *> \par Reference:
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 *> \date July 2012
228 *
229 *> \ingroup realOTHERcomputational
230 *
231 * =====================================================================
232  SUBROUTINE sorcsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
233  $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
234  $ ldv1t, work, lwork, iwork, info )
235 *
236 * -- LAPACK computational routine (version 3.5.0) --
237 * -- LAPACK is a software package provided by Univ. of Tennessee, --
238 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
239 * July 2012
240 *
241 * .. Scalar Arguments ..
242  CHARACTER jobu1, jobu2, jobv1t
243  INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
244  $ m, p, q
245 * ..
246 * .. Array Arguments ..
247  REAL theta(*)
248  REAL u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
249  $ x11(ldx11,*), x21(ldx21,*)
250  INTEGER iwork(*)
251 * ..
252 *
253 * =====================================================================
254 *
255 * .. Parameters ..
256  REAL one, zero
257  parameter( one = 1.0e0, zero = 0.0e0 )
258 * ..
259 * .. Local Scalars ..
260  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
261  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
262  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
263  $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
264  $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
265  $ lworkmin, lworkopt, r
266  LOGICAL lquery, wantu1, wantu2, wantv1t
267 * ..
268 * .. External Subroutines ..
269  EXTERNAL sbbcsd, scopy, slacpy, slapmr, slapmt, sorbdb1,
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. p ) THEN
301  info = -13
302  ELSE IF( wantu2 .AND. ldu2 .LT. m - p ) THEN
303  info = -15
304  ELSE IF( wantv1t .AND. ldv1t .LT. 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 * | SORBDB WORK | SORGQR WORK | SORGLQ WORK | B21D (R) |
323 * | | | | B21E (R-1) |
324 * | | | | B22D (R) |
325 * | | | | B22E (R-1) |
326 * | | | | SBBCSD 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  IF( r .EQ. q ) THEN
347  CALL sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
348  $ 0, 0, work, -1, childinfo )
349  lorbdb = int( work(1) )
350  IF( p .GE. m-p ) THEN
351  CALL sorgqr( p, p, q, u1, ldu1, 0, work(1), -1,
352  $ childinfo )
353  lorgqrmin = max( 1, p )
354  lorgqropt = int( work(1) )
355  ELSE
356  CALL sorgqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
357  $ childinfo )
358  lorgqrmin = max( 1, m-p )
359  lorgqropt = int( work(1) )
360  END IF
361  CALL sorglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
362  $ 0, work(1), -1, childinfo )
363  lorglqmin = max( 1, q-1 )
364  lorglqopt = int( work(1) )
365  CALL sbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
366  $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
367  $ 0, 0, 0, 0, 0, 0, work(1), -1, childinfo )
368  lbbcsd = int( work(1) )
369  ELSE IF( r .EQ. p ) THEN
370  CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
371  $ 0, 0, work(1), -1, childinfo )
372  lorbdb = int( work(1) )
373  IF( p-1 .GE. m-p ) THEN
374  CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
375  $ -1, childinfo )
376  lorgqrmin = max( 1, p-1 )
377  lorgqropt = int( work(1) )
378  ELSE
379  CALL sorgqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
380  $ childinfo )
381  lorgqrmin = max( 1, m-p )
382  lorgqropt = int( work(1) )
383  END IF
384  CALL sorglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
385  $ childinfo )
386  lorglqmin = max( 1, q )
387  lorglqopt = int( work(1) )
388  CALL sbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
389  $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
390  $ 0, 0, 0, 0, 0, 0, work(1), -1, childinfo )
391  lbbcsd = int( work(1) )
392  ELSE IF( r .EQ. m-p ) THEN
393  CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
394  $ 0, 0, work(1), -1, childinfo )
395  lorbdb = int( work(1) )
396  IF( p .GE. m-p-1 ) THEN
397  CALL sorgqr( p, p, q, u1, ldu1, 0, work(1), -1,
398  $ childinfo )
399  lorgqrmin = max( 1, p )
400  lorgqropt = int( work(1) )
401  ELSE
402  CALL sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
403  $ work(1), -1, childinfo )
404  lorgqrmin = max( 1, m-p-1 )
405  lorgqropt = int( work(1) )
406  END IF
407  CALL sorglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
408  $ childinfo )
409  lorglqmin = max( 1, q )
410  lorglqopt = int( work(1) )
411  CALL sbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
412  $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
413  $ 0, 0, 0, 0, 0, 0, 0, 0, work(1), -1,
414  $ childinfo )
415  lbbcsd = int( work(1) )
416  ELSE
417  CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
418  $ 0, 0, 0, work(1), -1, childinfo )
419  lorbdb = m + int( work(1) )
420  IF( p .GE. m-p ) THEN
421  CALL sorgqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
422  $ childinfo )
423  lorgqrmin = max( 1, p )
424  lorgqropt = int( work(1) )
425  ELSE
426  CALL sorgqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
427  $ childinfo )
428  lorgqrmin = max( 1, m-p )
429  lorgqropt = int( work(1) )
430  END IF
431  CALL sorglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
432  $ childinfo )
433  lorglqmin = max( 1, q )
434  lorglqopt = int( work(1) )
435  CALL sbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
436  $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
437  $ 0, 0, 0, 0, 0, 0, 0, 0, work(1), -1,
438  $ childinfo )
439  lbbcsd = int( work(1) )
440  END IF
441  lworkmin = max( iorbdb+lorbdb-1,
442  $ iorgqr+lorgqrmin-1,
443  $ iorglq+lorglqmin-1,
444  $ ibbcsd+lbbcsd-1 )
445  lworkopt = max( iorbdb+lorbdb-1,
446  $ iorgqr+lorgqropt-1,
447  $ iorglq+lorglqopt-1,
448  $ ibbcsd+lbbcsd-1 )
449  work(1) = lworkopt
450  IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
451  info = -19
452  END IF
453  END IF
454  IF( info .NE. 0 ) THEN
455  CALL xerbla( 'SORCSD2BY1', -info )
456  RETURN
457  ELSE IF( lquery ) THEN
458  RETURN
459  END IF
460  lorgqr = lwork-iorgqr+1
461  lorglq = lwork-iorglq+1
462 *
463 * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
464 * in which R = MIN(P,M-P,Q,M-Q)
465 *
466  IF( r .EQ. q ) THEN
467 *
468 * Case 1: R = Q
469 *
470 * Simultaneously bidiagonalize X11 and X21
471 *
472  CALL sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
473  $ work(iphi), work(itaup1), work(itaup2),
474  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
475 *
476 * Accumulate Householder reflectors
477 *
478  IF( wantu1 .AND. p .GT. 0 ) THEN
479  CALL slacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
480  CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
481  $ lorgqr, childinfo )
482  END IF
483  IF( wantu2 .AND. m-p .GT. 0 ) THEN
484  CALL slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
485  CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
486  $ work(iorgqr), lorgqr, childinfo )
487  END IF
488  IF( wantv1t .AND. q .GT. 0 ) THEN
489  v1t(1,1) = one
490  DO j = 2, q
491  v1t(1,j) = zero
492  v1t(j,1) = zero
493  END DO
494  CALL slacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
495  $ ldv1t )
496  CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
497  $ work(iorglq), lorglq, childinfo )
498  END IF
499 *
500 * Simultaneously diagonalize X11 and X21.
501 *
502  CALL sbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
503  $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
504  $ work(ib11d), work(ib11e), work(ib12d),
505  $ work(ib12e), work(ib21d), work(ib21e),
506  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
507  $ childinfo )
508 *
509 * Permute rows and columns to place zero submatrices in
510 * preferred positions
511 *
512  IF( q .GT. 0 .AND. wantu2 ) THEN
513  DO i = 1, q
514  iwork(i) = m - p - q + i
515  END DO
516  DO i = q + 1, m - p
517  iwork(i) = i - q
518  END DO
519  CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
520  END IF
521  ELSE IF( r .EQ. p ) THEN
522 *
523 * Case 2: R = P
524 *
525 * Simultaneously bidiagonalize X11 and X21
526 *
527  CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
528  $ work(iphi), work(itaup1), work(itaup2),
529  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
530 *
531 * Accumulate Householder reflectors
532 *
533  IF( wantu1 .AND. p .GT. 0 ) THEN
534  u1(1,1) = one
535  DO j = 2, p
536  u1(1,j) = zero
537  u1(j,1) = zero
538  END DO
539  CALL slacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
540  CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
541  $ work(iorgqr), lorgqr, childinfo )
542  END IF
543  IF( wantu2 .AND. m-p .GT. 0 ) THEN
544  CALL slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
545  CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
546  $ work(iorgqr), lorgqr, childinfo )
547  END IF
548  IF( wantv1t .AND. q .GT. 0 ) THEN
549  CALL slacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
550  CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
551  $ work(iorglq), lorglq, childinfo )
552  END IF
553 *
554 * Simultaneously diagonalize X11 and X21.
555 *
556  CALL sbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
557  $ work(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
558  $ work(ib11d), work(ib11e), work(ib12d),
559  $ work(ib12e), work(ib21d), work(ib21e),
560  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
561  $ childinfo )
562 *
563 * Permute rows and columns to place identity submatrices in
564 * preferred positions
565 *
566  IF( q .GT. 0 .AND. wantu2 ) THEN
567  DO i = 1, q
568  iwork(i) = m - p - q + i
569  END DO
570  DO i = q + 1, m - p
571  iwork(i) = i - q
572  END DO
573  CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
574  END IF
575  ELSE IF( r .EQ. m-p ) THEN
576 *
577 * Case 3: R = M-P
578 *
579 * Simultaneously bidiagonalize X11 and X21
580 *
581  CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
582  $ work(iphi), work(itaup1), work(itaup2),
583  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
584 *
585 * Accumulate Householder reflectors
586 *
587  IF( wantu1 .AND. p .GT. 0 ) THEN
588  CALL slacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
589  CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
590  $ lorgqr, childinfo )
591  END IF
592  IF( wantu2 .AND. m-p .GT. 0 ) THEN
593  u2(1,1) = one
594  DO j = 2, m-p
595  u2(1,j) = zero
596  u2(j,1) = zero
597  END DO
598  CALL slacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
599  $ ldu2 )
600  CALL sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
601  $ work(itaup2), work(iorgqr), lorgqr, childinfo )
602  END IF
603  IF( wantv1t .AND. q .GT. 0 ) THEN
604  CALL slacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
605  CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
606  $ work(iorglq), lorglq, childinfo )
607  END IF
608 *
609 * Simultaneously diagonalize X11 and X21.
610 *
611  CALL sbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
612  $ theta, work(iphi), 0, 1, v1t, ldv1t, u2, ldu2, u1,
613  $ ldu1, work(ib11d), work(ib11e), work(ib12d),
614  $ work(ib12e), work(ib21d), work(ib21e),
615  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
616  $ childinfo )
617 *
618 * Permute rows and columns to place identity submatrices in
619 * preferred positions
620 *
621  IF( q .GT. r ) THEN
622  DO i = 1, r
623  iwork(i) = q - r + i
624  END DO
625  DO i = r + 1, q
626  iwork(i) = i - r
627  END DO
628  IF( wantu1 ) THEN
629  CALL slapmt( .false., p, q, u1, ldu1, iwork )
630  END IF
631  IF( wantv1t ) THEN
632  CALL slapmr( .false., q, q, v1t, ldv1t, iwork )
633  END IF
634  END IF
635  ELSE
636 *
637 * Case 4: R = M-Q
638 *
639 * Simultaneously bidiagonalize X11 and X21
640 *
641  CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
642  $ work(iphi), work(itaup1), work(itaup2),
643  $ work(itauq1), work(iorbdb), work(iorbdb+m),
644  $ lorbdb-m, childinfo )
645 *
646 * Accumulate Householder reflectors
647 *
648  IF( wantu1 .AND. p .GT. 0 ) THEN
649  CALL scopy( p, work(iorbdb), 1, u1, 1 )
650  DO j = 2, p
651  u1(1,j) = zero
652  END DO
653  CALL slacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
654  $ ldu1 )
655  CALL sorgqr( p, p, m-q, u1, ldu1, work(itaup1),
656  $ work(iorgqr), lorgqr, childinfo )
657  END IF
658  IF( wantu2 .AND. m-p .GT. 0 ) THEN
659  CALL scopy( m-p, work(iorbdb+p), 1, u2, 1 )
660  DO j = 2, m-p
661  u2(1,j) = zero
662  END DO
663  CALL slacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
664  $ ldu2 )
665  CALL sorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
666  $ work(iorgqr), lorgqr, childinfo )
667  END IF
668  IF( wantv1t .AND. q .GT. 0 ) THEN
669  CALL slacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
670  CALL slacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
671  $ v1t(m-q+1,m-q+1), ldv1t )
672  CALL slacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
673  $ v1t(p+1,p+1), ldv1t )
674  CALL sorglq( q, q, q, v1t, ldv1t, work(itauq1),
675  $ work(iorglq), lorglq, childinfo )
676  END IF
677 *
678 * Simultaneously diagonalize X11 and X21.
679 *
680  CALL sbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
681  $ theta, work(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
682  $ ldv1t, work(ib11d), work(ib11e), work(ib12d),
683  $ work(ib12e), work(ib21d), work(ib21e),
684  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
685  $ childinfo )
686 *
687 * Permute rows and columns to place identity submatrices in
688 * preferred positions
689 *
690  IF( p .GT. r ) THEN
691  DO i = 1, r
692  iwork(i) = p - r + i
693  END DO
694  DO i = r + 1, p
695  iwork(i) = i - r
696  END DO
697  IF( wantu1 ) THEN
698  CALL slapmt( .false., p, p, u1, ldu1, iwork )
699  END IF
700  IF( wantv1t ) THEN
701  CALL slapmr( .false., p, q, v1t, ldv1t, iwork )
702  END IF
703  END IF
704  END IF
705 *
706  RETURN
707 *
708 * End of SORCSD2BY1
709 *
710  END
711