LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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, and V1 are P-by-P,
58 *> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
59 *> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
60 *> R = MIN(P,M-P,Q,M-Q).
61 *> \endverbatim
62 *
63 * Arguments:
64 * ==========
65 *
66 *> \param[in] JOBU1
67 *> \verbatim
68 *> JOBU1 is CHARACTER
69 *> = 'Y': U1 is computed;
70 *> otherwise: U1 is not computed.
71 *> \endverbatim
72 *>
73 *> \param[in] JOBU2
74 *> \verbatim
75 *> JOBU2 is CHARACTER
76 *> = 'Y': U2 is computed;
77 *> otherwise: U2 is not computed.
78 *> \endverbatim
79 *>
80 *> \param[in] JOBV1T
81 *> \verbatim
82 *> JOBV1T is CHARACTER
83 *> = 'Y': V1T is computed;
84 *> otherwise: V1T is not computed.
85 *> \endverbatim
86 *>
87 *> \param[in] M
88 *> \verbatim
89 *> M is INTEGER
90 *> The number of rows in X.
91 *> \endverbatim
92 *>
93 *> \param[in] P
94 *> \verbatim
95 *> P is INTEGER
96 *> The number of rows in X11. 0 <= P <= M.
97 *> \endverbatim
98 *>
99 *> \param[in] Q
100 *> \verbatim
101 *> Q is INTEGER
102 *> The number of columns in X11 and X21. 0 <= Q <= M.
103 *> \endverbatim
104 *>
105 *> \param[in,out] X11
106 *> \verbatim
107 *> X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
108 *> On entry, part of the orthogonal matrix whose CSD is 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 DOUBLE PRECISION array, dimension (LDX21,Q)
120 *> On entry, part of the orthogonal matrix whose CSD is desired.
121 *> \endverbatim
122 *>
123 *> \param[in] LDX21
124 *> \verbatim
125 *> LDX21 is INTEGER
126 *> The leading dimension of X21. LDX21 >= MAX(1,M-P).
127 *> \endverbatim
128 *>
129 *> \param[out] THETA
130 *> \verbatim
131 *> THETA is DOUBLE PRECISION array, dimension (R), in which R =
132 *> MIN(P,M-P,Q,M-Q).
133 *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
134 *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
135 *> \endverbatim
136 *>
137 *> \param[out] U1
138 *> \verbatim
139 *> U1 is DOUBLE PRECISION array, dimension (P)
140 *> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
141 *> \endverbatim
142 *>
143 *> \param[in] LDU1
144 *> \verbatim
145 *> LDU1 is INTEGER
146 *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
147 *> MAX(1,P).
148 *> \endverbatim
149 *>
150 *> \param[out] U2
151 *> \verbatim
152 *> U2 is DOUBLE PRECISION array, dimension (M-P)
153 *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
154 *> matrix U2.
155 *> \endverbatim
156 *>
157 *> \param[in] LDU2
158 *> \verbatim
159 *> LDU2 is INTEGER
160 *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
161 *> MAX(1,M-P).
162 *> \endverbatim
163 *>
164 *> \param[out] V1T
165 *> \verbatim
166 *> V1T is DOUBLE PRECISION array, dimension (Q)
167 *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
168 *> matrix V1**T.
169 *> \endverbatim
170 *>
171 *> \param[in] LDV1T
172 *> \verbatim
173 *> LDV1T is INTEGER
174 *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
175 *> MAX(1,Q).
176 *> \endverbatim
177 *>
178 *> \param[out] WORK
179 *> \verbatim
180 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
181 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
182 *> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
183 *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
184 *> define the matrix in intermediate bidiagonal-block form
185 *> remaining after nonconvergence. INFO specifies the number
186 *> of nonzero PHI's.
187 *> \endverbatim
188 *>
189 *> \param[in] LWORK
190 *> \verbatim
191 *> LWORK is INTEGER
192 *> The dimension of the array WORK.
193 *>
194 *> If LWORK = -1, then a workspace query is assumed; the routine
195 *> only calculates the optimal size of the WORK array, returns
196 *> this value as the first entry of the work array, and no error
197 *> message related to LWORK is issued by XERBLA.
198 *> \endverbatim
199 *>
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: DBBCSD did not converge. See the description of WORK
211 *> above for details.
212 *> \endverbatim
213 *
214 *> \par References:
215 * ================
216 *>
217 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
218 *> Algorithms, 50(1):33-65, 2009.
219 *
220 * Authors:
221 * ========
222 *
223 *> \author Univ. of Tennessee
224 *> \author Univ. of California Berkeley
225 *> \author Univ. of Colorado Denver
226 *> \author NAG Ltd.
227 *
228 *> \date July 2012
229 *
230 *> \ingroup doubleOTHERcomputational
231 *
232 * =====================================================================
233  SUBROUTINE dorcsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
234  $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
235  $ ldv1t, work, lwork, iwork, info )
236 *
237 * -- LAPACK computational routine (3.5.0) --
238 * -- LAPACK is a software package provided by Univ. of Tennessee, --
239 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
240 * July 2012
241 *
242 * .. Scalar Arguments ..
243  CHARACTER JOBU1, JOBU2, JOBV1T
244  INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
245  $ m, p, q
246 * ..
247 * .. Array Arguments ..
248  DOUBLE PRECISION THETA(*)
249  DOUBLE PRECISION U1(ldu1,*), U2(ldu2,*), V1T(ldv1t,*), WORK(*),
250  $ x11(ldx11,*), x21(ldx21,*)
251  INTEGER IWORK(*)
252 * ..
253 *
254 * =====================================================================
255 *
256 * .. Parameters ..
257  DOUBLE PRECISION ONE, ZERO
258  parameter ( one = 1.0d0, zero = 0.0d0 )
259 * ..
260 * .. Local Scalars ..
261  INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
262  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
263  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
264  $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
265  $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
266  $ lworkmin, lworkopt, r
267  LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
268 * ..
269 * .. Local Arrays ..
270  DOUBLE PRECISION DUM1(1), DUM2(1,1)
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. max( 1, p ) ) THEN
305  info = -13
306  ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) ) THEN
307  info = -15
308  ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, 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  lorgqrmin = 1
351  lorgqropt = 1
352  lorglqmin = 1
353  lorglqopt = 1
354  IF( r .EQ. q ) THEN
355  CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
356  $ dum1, dum1, dum1, dum1, work,
357  $ -1, childinfo )
358  lorbdb = int( work(1) )
359  IF( wantu1 .AND. p .GT. 0 ) THEN
360  CALL dorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
361  $ childinfo )
362  lorgqrmin = max( lorgqrmin, p )
363  lorgqropt = max( lorgqropt, int( work(1) ) )
364  ENDIF
365  IF( wantu2 .AND. m-p .GT. 0 ) THEN
366  CALL dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
367  $ -1, childinfo )
368  lorgqrmin = max( lorgqrmin, m-p )
369  lorgqropt = max( lorgqropt, int( work(1) ) )
370  END IF
371  IF( wantv1t .AND. q .GT. 0 ) THEN
372  CALL dorglq( q-1, q-1, q-1, v1t, ldv1t,
373  $ dum1, work(1), -1, childinfo )
374  lorglqmin = max( lorglqmin, q-1 )
375  lorglqopt = max( lorglqopt, int( work(1) ) )
376  END IF
377  CALL dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
378  $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t,
379  $ dum2, 1, dum1, dum1, dum1,
380  $ dum1, dum1, dum1, dum1,
381  $ dum1, work(1), -1, childinfo )
382  lbbcsd = int( work(1) )
383  ELSE IF( r .EQ. p ) THEN
384  CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
385  $ dum1, dum1, dum1, dum1,
386  $ work(1), -1, childinfo )
387  lorbdb = int( work(1) )
388  IF( wantu1 .AND. p .GT. 0 ) THEN
389  CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, dum1,
390  $ work(1), -1, childinfo )
391  lorgqrmin = max( lorgqrmin, p-1 )
392  lorgqropt = max( lorgqropt, int( work(1) ) )
393  END IF
394  IF( wantu2 .AND. m-p .GT. 0 ) THEN
395  CALL dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
396  $ -1, childinfo )
397  lorgqrmin = max( lorgqrmin, m-p )
398  lorgqropt = max( lorgqropt, int( work(1) ) )
399  END IF
400  IF( wantv1t .AND. q .GT. 0 ) THEN
401  CALL dorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
402  $ childinfo )
403  lorglqmin = max( lorglqmin, q )
404  lorglqopt = max( lorglqopt, int( work(1) ) )
405  END IF
406  CALL dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
407  $ dum1, v1t, ldv1t, dum2, 1, u1, ldu1,
408  $ u2, ldu2, dum1, dum1, dum1,
409  $ dum1, dum1, dum1, dum1,
410  $ dum1, work(1), -1, childinfo )
411  lbbcsd = int( work(1) )
412  ELSE IF( r .EQ. m-p ) THEN
413  CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
414  $ dum1, dum1, dum1, dum1,
415  $ work(1), -1, childinfo )
416  lorbdb = int( work(1) )
417  IF( wantu1 .AND. p .GT. 0 ) THEN
418  CALL dorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
419  $ childinfo )
420  lorgqrmin = max( lorgqrmin, p )
421  lorgqropt = max( lorgqropt, int( work(1) ) )
422  END IF
423  IF( wantu2 .AND. m-p .GT. 0 ) THEN
424  CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
425  $ dum1, work(1), -1, childinfo )
426  lorgqrmin = max( lorgqrmin, m-p-1 )
427  lorgqropt = max( lorgqropt, int( work(1) ) )
428  END IF
429  IF( wantv1t .AND. q .GT. 0 ) THEN
430  CALL dorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
431  $ childinfo )
432  lorglqmin = max( lorglqmin, q )
433  lorglqopt = max( lorglqopt, int( work(1) ) )
434  END IF
435  CALL dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
436  $ theta, dum1, dum2, 1, v1t, ldv1t, u2,
437  $ ldu2, u1, ldu1, dum1, dum1, dum1,
438  $ dum1, dum1, dum1, dum1,
439  $ dum1, work(1), -1, childinfo )
440  lbbcsd = int( work(1) )
441  ELSE
442  CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
443  $ dum1, dum1, dum1, dum1,
444  $ dum1, work(1), -1, childinfo )
445  lorbdb = m + int( work(1) )
446  IF( wantu1 .AND. p .GT. 0 ) THEN
447  CALL dorgqr( p, p, m-q, u1, ldu1, dum1, work(1), -1,
448  $ childinfo )
449  lorgqrmin = max( lorgqrmin, p )
450  lorgqropt = max( lorgqropt, int( work(1) ) )
451  END IF
452  IF( wantu2 .AND. m-p .GT. 0 ) THEN
453  CALL dorgqr( m-p, m-p, m-q, u2, ldu2, dum1, work(1),
454  $ -1, childinfo )
455  lorgqrmin = max( lorgqrmin, m-p )
456  lorgqropt = max( lorgqropt, int( work(1) ) )
457  END IF
458  IF( wantv1t .AND. q .GT. 0 ) THEN
459  CALL dorglq( q, q, q, v1t, ldv1t, dum1, work(1), -1,
460  $ childinfo )
461  lorglqmin = max( lorglqmin, q )
462  lorglqopt = max( lorglqopt, int( work(1) ) )
463  END IF
464  CALL dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
465  $ theta, dum1, u2, ldu2, u1, ldu1, dum2,
466  $ 1, v1t, ldv1t, dum1, dum1, dum1,
467  $ dum1, dum1, dum1, dum1,
468  $ dum1, work(1), -1, childinfo )
469  lbbcsd = int( work(1) )
470  END IF
471  lworkmin = max( iorbdb+lorbdb-1,
472  $ iorgqr+lorgqrmin-1,
473  $ iorglq+lorglqmin-1,
474  $ ibbcsd+lbbcsd-1 )
475  lworkopt = max( iorbdb+lorbdb-1,
476  $ iorgqr+lorgqropt-1,
477  $ iorglq+lorglqopt-1,
478  $ ibbcsd+lbbcsd-1 )
479  work(1) = lworkopt
480  IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
481  info = -19
482  END IF
483  END IF
484  IF( info .NE. 0 ) THEN
485  CALL xerbla( 'DORCSD2BY1', -info )
486  RETURN
487  ELSE IF( lquery ) THEN
488  RETURN
489  END IF
490  lorgqr = lwork-iorgqr+1
491  lorglq = lwork-iorglq+1
492 *
493 * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
494 * in which R = MIN(P,M-P,Q,M-Q)
495 *
496  IF( r .EQ. q ) THEN
497 *
498 * Case 1: R = Q
499 *
500 * Simultaneously bidiagonalize X11 and X21
501 *
502  CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
503  $ work(iphi), work(itaup1), work(itaup2),
504  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
505 *
506 * Accumulate Householder reflectors
507 *
508  IF( wantu1 .AND. p .GT. 0 ) THEN
509  CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
510  CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
511  $ lorgqr, childinfo )
512  END IF
513  IF( wantu2 .AND. m-p .GT. 0 ) THEN
514  CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
515  CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
516  $ work(iorgqr), lorgqr, childinfo )
517  END IF
518  IF( wantv1t .AND. q .GT. 0 ) THEN
519  v1t(1,1) = one
520  DO j = 2, q
521  v1t(1,j) = zero
522  v1t(j,1) = zero
523  END DO
524  CALL dlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
525  $ ldv1t )
526  CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
527  $ work(iorglq), lorglq, childinfo )
528  END IF
529 *
530 * Simultaneously diagonalize X11 and X21.
531 *
532  CALL dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
533  $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t,
534  $ dum2, 1, work(ib11d), work(ib11e),
535  $ work(ib12d), work(ib12e), work(ib21d),
536  $ work(ib21e), work(ib22d), work(ib22e),
537  $ work(ibbcsd), lbbcsd, childinfo )
538 *
539 * Permute rows and columns to place zero submatrices in
540 * preferred positions
541 *
542  IF( q .GT. 0 .AND. wantu2 ) THEN
543  DO i = 1, q
544  iwork(i) = m - p - q + i
545  END DO
546  DO i = q + 1, m - p
547  iwork(i) = i - q
548  END DO
549  CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
550  END IF
551  ELSE IF( r .EQ. p ) THEN
552 *
553 * Case 2: R = P
554 *
555 * Simultaneously bidiagonalize X11 and X21
556 *
557  CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
558  $ work(iphi), work(itaup1), work(itaup2),
559  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
560 *
561 * Accumulate Householder reflectors
562 *
563  IF( wantu1 .AND. p .GT. 0 ) THEN
564  u1(1,1) = one
565  DO j = 2, p
566  u1(1,j) = zero
567  u1(j,1) = zero
568  END DO
569  CALL dlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
570  CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
571  $ work(iorgqr), lorgqr, childinfo )
572  END IF
573  IF( wantu2 .AND. m-p .GT. 0 ) THEN
574  CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
575  CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
576  $ work(iorgqr), lorgqr, childinfo )
577  END IF
578  IF( wantv1t .AND. q .GT. 0 ) THEN
579  CALL dlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
580  CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
581  $ work(iorglq), lorglq, childinfo )
582  END IF
583 *
584 * Simultaneously diagonalize X11 and X21.
585 *
586  CALL dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
587  $ work(iphi), v1t, ldv1t, dum2, 1, u1, ldu1, u2,
588  $ ldu2, work(ib11d), work(ib11e), work(ib12d),
589  $ work(ib12e), work(ib21d), work(ib21e),
590  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
591  $ childinfo )
592 *
593 * Permute rows and columns to place identity submatrices in
594 * preferred positions
595 *
596  IF( q .GT. 0 .AND. wantu2 ) THEN
597  DO i = 1, q
598  iwork(i) = m - p - q + i
599  END DO
600  DO i = q + 1, m - p
601  iwork(i) = i - q
602  END DO
603  CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
604  END IF
605  ELSE IF( r .EQ. m-p ) THEN
606 *
607 * Case 3: R = M-P
608 *
609 * Simultaneously bidiagonalize X11 and X21
610 *
611  CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
612  $ work(iphi), work(itaup1), work(itaup2),
613  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
614 *
615 * Accumulate Householder reflectors
616 *
617  IF( wantu1 .AND. p .GT. 0 ) THEN
618  CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
619  CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
620  $ lorgqr, childinfo )
621  END IF
622  IF( wantu2 .AND. m-p .GT. 0 ) THEN
623  u2(1,1) = one
624  DO j = 2, m-p
625  u2(1,j) = zero
626  u2(j,1) = zero
627  END DO
628  CALL dlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
629  $ ldu2 )
630  CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
631  $ work(itaup2), work(iorgqr), lorgqr, childinfo )
632  END IF
633  IF( wantv1t .AND. q .GT. 0 ) THEN
634  CALL dlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
635  CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
636  $ work(iorglq), lorglq, childinfo )
637  END IF
638 *
639 * Simultaneously diagonalize X11 and X21.
640 *
641  CALL dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
642  $ theta, work(iphi), dum2, 1, v1t, ldv1t, u2,
643  $ ldu2, u1, ldu1, work(ib11d), work(ib11e),
644  $ work(ib12d), work(ib12e), work(ib21d),
645  $ work(ib21e), work(ib22d), work(ib22e),
646  $ work(ibbcsd), lbbcsd, childinfo )
647 *
648 * Permute rows and columns to place identity submatrices in
649 * preferred positions
650 *
651  IF( q .GT. r ) THEN
652  DO i = 1, r
653  iwork(i) = q - r + i
654  END DO
655  DO i = r + 1, q
656  iwork(i) = i - r
657  END DO
658  IF( wantu1 ) THEN
659  CALL dlapmt( .false., p, q, u1, ldu1, iwork )
660  END IF
661  IF( wantv1t ) THEN
662  CALL dlapmr( .false., q, q, v1t, ldv1t, iwork )
663  END IF
664  END IF
665  ELSE
666 *
667 * Case 4: R = M-Q
668 *
669 * Simultaneously bidiagonalize X11 and X21
670 *
671  CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
672  $ work(iphi), work(itaup1), work(itaup2),
673  $ work(itauq1), work(iorbdb), work(iorbdb+m),
674  $ lorbdb-m, childinfo )
675 *
676 * Accumulate Householder reflectors
677 *
678  IF( wantu1 .AND. p .GT. 0 ) THEN
679  CALL dcopy( p, work(iorbdb), 1, u1, 1 )
680  DO j = 2, p
681  u1(1,j) = zero
682  END DO
683  CALL dlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
684  $ ldu1 )
685  CALL dorgqr( p, p, m-q, u1, ldu1, work(itaup1),
686  $ work(iorgqr), lorgqr, childinfo )
687  END IF
688  IF( wantu2 .AND. m-p .GT. 0 ) THEN
689  CALL dcopy( m-p, work(iorbdb+p), 1, u2, 1 )
690  DO j = 2, m-p
691  u2(1,j) = zero
692  END DO
693  CALL dlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
694  $ ldu2 )
695  CALL dorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
696  $ work(iorgqr), lorgqr, childinfo )
697  END IF
698  IF( wantv1t .AND. q .GT. 0 ) THEN
699  CALL dlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
700  CALL dlacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
701  $ v1t(m-q+1,m-q+1), ldv1t )
702  CALL dlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
703  $ v1t(p+1,p+1), ldv1t )
704  CALL dorglq( q, q, q, v1t, ldv1t, work(itauq1),
705  $ work(iorglq), lorglq, childinfo )
706  END IF
707 *
708 * Simultaneously diagonalize X11 and X21.
709 *
710  CALL dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
711  $ theta, work(iphi), u2, ldu2, u1, ldu1, dum2,
712  $ 1, v1t, ldv1t, work(ib11d), work(ib11e),
713  $ work(ib12d), work(ib12e), work(ib21d),
714  $ work(ib21e), work(ib22d), work(ib22e),
715  $ work(ibbcsd), lbbcsd, childinfo )
716 *
717 * Permute rows and columns to place identity submatrices in
718 * preferred positions
719 *
720  IF( p .GT. r ) THEN
721  DO i = 1, r
722  iwork(i) = p - r + i
723  END DO
724  DO i = r + 1, p
725  iwork(i) = i - r
726  END DO
727  IF( wantu1 ) THEN
728  CALL dlapmt( .false., p, p, u1, ldu1, iwork )
729  END IF
730  IF( wantv1t ) THEN
731  CALL dlapmr( .false., p, q, v1t, ldv1t, iwork )
732  END IF
733  END IF
734  END IF
735 *
736  RETURN
737 *
738 * End of DORCSD2BY1
739 *
740  END
741 
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
Definition: dorglq.f:129
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlapmt(FORWRD, M, N, X, LDX, K)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: dlapmt.f:106
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB1
Definition: dorbdb1.f:205
subroutine dorbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB3
Definition: dorbdb3.f:203
subroutine dorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB2
Definition: dorbdb2.f:204
subroutine dorbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
DORBDB4
Definition: dorbdb4.f:215
subroutine dlapmr(FORWRD, M, N, X, LDX, K)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: dlapmr.f:106
subroutine dbbcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, B22D, B22E, WORK, LWORK, INFO)
DBBCSD
Definition: dbbcsd.f:334
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
Definition: dorgqr.f:130
subroutine dorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
DORCSD2BY1
Definition: dorcsd2by1.f:236