LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dorcsd.f
Go to the documentation of this file.
1 *> \brief \b DORCSD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DORCSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * RECURSIVE SUBROUTINE DORCSD( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
22 * SIGNS, M, P, Q, X11, LDX11, X12,
23 * LDX12, X21, LDX21, X22, LDX22, THETA,
24 * U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
25 * LDV2T, WORK, LWORK, IWORK, INFO )
26 *
27 * .. Scalar Arguments ..
28 * CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
29 * INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
30 * $ LDX21, LDX22, LWORK, M, P, Q
31 * ..
32 * .. Array Arguments ..
33 * INTEGER IWORK( * )
34 * DOUBLE PRECISION THETA( * )
35 * DOUBLE PRECISION U1( LDU1, * ), U2( LDU2, * ), V1T( LDV1T, * ),
36 * $ V2T( LDV2T, * ), WORK( * ), X11( LDX11, * ),
37 * $ X12( LDX12, * ), X21( LDX21, * ), X22( LDX22,
38 * $ * )
39 * ..
40 *
41 *
42 *> \par Purpose:
43 * =============
44 *>
45 *> \verbatim
46 *>
47 *> DORCSD computes the CS decomposition of an M-by-M partitioned
48 *> orthogonal matrix X:
49 *>
50 *> [ I 0 0 | 0 0 0 ]
51 *> [ 0 C 0 | 0 -S 0 ]
52 *> [ X11 | X12 ] [ U1 | ] [ 0 0 0 | 0 0 -I ] [ V1 | ]**T
53 *> X = [-----------] = [---------] [---------------------] [---------] .
54 *> [ X21 | X22 ] [ | U2 ] [ 0 0 0 | I 0 0 ] [ | V2 ]
55 *> [ 0 S 0 | 0 C 0 ]
56 *> [ 0 0 I | 0 0 0 ]
57 *>
58 *> X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P,
59 *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
60 *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
61 *> which R = MIN(P,M-P,Q,M-Q).
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] JOBV2T
89 *> \verbatim
90 *> JOBV2T is CHARACTER
91 *> = 'Y': V2T is computed;
92 *> otherwise: V2T is not computed.
93 *> \endverbatim
94 *>
95 *> \param[in] TRANS
96 *> \verbatim
97 *> TRANS is CHARACTER
98 *> = 'T': X, U1, U2, V1T, and V2T are stored in row-major
99 *> order;
100 *> otherwise: X, U1, U2, V1T, and V2T are stored in column-
101 *> major order.
102 *> \endverbatim
103 *>
104 *> \param[in] SIGNS
105 *> \verbatim
106 *> SIGNS is CHARACTER
107 *> = 'O': The lower-left block is made nonpositive (the
108 *> "other" convention);
109 *> otherwise: The upper-right block is made nonpositive (the
110 *> "default" convention).
111 *> \endverbatim
112 *>
113 *> \param[in] M
114 *> \verbatim
115 *> M is INTEGER
116 *> The number of rows and columns in X.
117 *> \endverbatim
118 *>
119 *> \param[in] P
120 *> \verbatim
121 *> P is INTEGER
122 *> The number of rows in X11 and X12. 0 <= P <= M.
123 *> \endverbatim
124 *>
125 *> \param[in] Q
126 *> \verbatim
127 *> Q is INTEGER
128 *> The number of columns in X11 and X21. 0 <= Q <= M.
129 *> \endverbatim
130 *>
131 *> \param[in,out] X11
132 *> \verbatim
133 *> X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
134 *> On entry, part of the orthogonal matrix whose CSD is desired.
135 *> \endverbatim
136 *>
137 *> \param[in] LDX11
138 *> \verbatim
139 *> LDX11 is INTEGER
140 *> The leading dimension of X11. LDX11 >= MAX(1,P).
141 *> \endverbatim
142 *>
143 *> \param[in,out] X12
144 *> \verbatim
145 *> X12 is DOUBLE PRECISION array, dimension (LDX12,M-Q)
146 *> On entry, part of the orthogonal matrix whose CSD is desired.
147 *> \endverbatim
148 *>
149 *> \param[in] LDX12
150 *> \verbatim
151 *> LDX12 is INTEGER
152 *> The leading dimension of X12. LDX12 >= MAX(1,P).
153 *> \endverbatim
154 *>
155 *> \param[in,out] X21
156 *> \verbatim
157 *> X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
158 *> On entry, part of the orthogonal matrix whose CSD is desired.
159 *> \endverbatim
160 *>
161 *> \param[in] LDX21
162 *> \verbatim
163 *> LDX21 is INTEGER
164 *> The leading dimension of X11. LDX21 >= MAX(1,M-P).
165 *> \endverbatim
166 *>
167 *> \param[in,out] X22
168 *> \verbatim
169 *> X22 is DOUBLE PRECISION array, dimension (LDX22,M-Q)
170 *> On entry, part of the orthogonal matrix whose CSD is desired.
171 *> \endverbatim
172 *>
173 *> \param[in] LDX22
174 *> \verbatim
175 *> LDX22 is INTEGER
176 *> The leading dimension of X11. LDX22 >= MAX(1,M-P).
177 *> \endverbatim
178 *>
179 *> \param[out] THETA
180 *> \verbatim
181 *> THETA is DOUBLE PRECISION array, dimension (R), in which R =
182 *> MIN(P,M-P,Q,M-Q).
183 *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
184 *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
185 *> \endverbatim
186 *>
187 *> \param[out] U1
188 *> \verbatim
189 *> U1 is DOUBLE PRECISION array, dimension (P)
190 *> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
191 *> \endverbatim
192 *>
193 *> \param[in] LDU1
194 *> \verbatim
195 *> LDU1 is INTEGER
196 *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
197 *> MAX(1,P).
198 *> \endverbatim
199 *>
200 *> \param[out] U2
201 *> \verbatim
202 *> U2 is DOUBLE PRECISION array, dimension (M-P)
203 *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
204 *> matrix U2.
205 *> \endverbatim
206 *>
207 *> \param[in] LDU2
208 *> \verbatim
209 *> LDU2 is INTEGER
210 *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
211 *> MAX(1,M-P).
212 *> \endverbatim
213 *>
214 *> \param[out] V1T
215 *> \verbatim
216 *> V1T is DOUBLE PRECISION array, dimension (Q)
217 *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
218 *> matrix V1**T.
219 *> \endverbatim
220 *>
221 *> \param[in] LDV1T
222 *> \verbatim
223 *> LDV1T is INTEGER
224 *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
225 *> MAX(1,Q).
226 *> \endverbatim
227 *>
228 *> \param[out] V2T
229 *> \verbatim
230 *> V2T is DOUBLE PRECISION array, dimension (M-Q)
231 *> If JOBV2T = 'Y', V2T contains the (M-Q)-by-(M-Q) orthogonal
232 *> matrix V2**T.
233 *> \endverbatim
234 *>
235 *> \param[in] LDV2T
236 *> \verbatim
237 *> LDV2T is INTEGER
238 *> The leading dimension of V2T. If JOBV2T = 'Y', LDV2T >=
239 *> MAX(1,M-Q).
240 *> \endverbatim
241 *>
242 *> \param[out] WORK
243 *> \verbatim
244 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
245 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
246 *> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
247 *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
248 *> define the matrix in intermediate bidiagonal-block form
249 *> remaining after nonconvergence. INFO specifies the number
250 *> of nonzero PHI's.
251 *> \endverbatim
252 *>
253 *> \param[in] LWORK
254 *> \verbatim
255 *> LWORK is INTEGER
256 *> The dimension of the array WORK.
257 *>
258 *> If LWORK = -1, then a workspace query is assumed; the routine
259 *> only calculates the optimal size of the WORK array, returns
260 *> this value as the first entry of the work array, and no error
261 *> message related to LWORK is issued by XERBLA.
262 *> \endverbatim
263 *>
264 *> \param[out] IWORK
265 *> \verbatim
266 *> IWORK is INTEGER array, dimension (M-MIN(P, M-P, Q, M-Q))
267 *> \endverbatim
268 *>
269 *> \param[out] INFO
270 *> \verbatim
271 *> INFO is INTEGER
272 *> = 0: successful exit.
273 *> < 0: if INFO = -i, the i-th argument had an illegal value.
274 *> > 0: DBBCSD did not converge. See the description of WORK
275 *> above for details.
276 *> \endverbatim
277 *
278 *> \par References:
279 * ================
280 *>
281 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
282 *> Algorithms, 50(1):33-65, 2009.
283 *
284 * Authors:
285 * ========
286 *
287 *> \author Univ. of Tennessee
288 *> \author Univ. of California Berkeley
289 *> \author Univ. of Colorado Denver
290 *> \author NAG Ltd.
291 *
292 *> \date November 2013
293 *
294 *> \ingroup doubleOTHERcomputational
295 *
296 * =====================================================================
297  RECURSIVE SUBROUTINE dorcsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
298  $ signs, m, p, q, x11, ldx11, x12,
299  $ ldx12, x21, ldx21, x22, ldx22, theta,
300  $ u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
301  $ ldv2t, work, lwork, iwork, info )
302 *
303 * -- LAPACK computational routine (version 3.5.0) --
304 * -- LAPACK is a software package provided by Univ. of Tennessee, --
305 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
306 * November 2013
307 *
308 * .. Scalar Arguments ..
309  CHARACTER jobu1, jobu2, jobv1t, jobv2t, signs, trans
310  INTEGER info, ldu1, ldu2, ldv1t, ldv2t, ldx11, ldx12,
311  $ ldx21, ldx22, lwork, m, p, q
312 * ..
313 * .. Array Arguments ..
314  INTEGER iwork( * )
315  DOUBLE PRECISION theta( * )
316  DOUBLE PRECISION u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
317  $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
318  $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
319  $ * )
320 * ..
321 *
322 * ===================================================================
323 *
324 * .. Parameters ..
325  DOUBLE PRECISION one, zero
326  parameter( one = 1.0d0,
327  $ zero = 0.0d0 )
328 * ..
329 * .. Local Scalars ..
330  CHARACTER transt, signst
331  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
332  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
333  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
334  $ itauq2, j, lbbcsdwork, lbbcsdworkmin,
335  $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
336  $ lorbdbworkopt, lorglqwork, lorglqworkmin,
337  $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
338  $ lorgqrworkopt, lworkmin, lworkopt
339  LOGICAL colmajor, defaultsigns, lquery, wantu1, wantu2,
340  $ wantv1t, wantv2t
341 * ..
342 * .. External Subroutines ..
343  EXTERNAL dbbcsd, dlacpy, dlapmr, dlapmt, dlascl, dlaset,
345 * ..
346 * .. External Functions ..
347  LOGICAL lsame
348  EXTERNAL lsame
349 * ..
350 * .. Intrinsic Functions
351  INTRINSIC int, max, min
352 * ..
353 * .. Executable Statements ..
354 *
355 * Test input arguments
356 *
357  info = 0
358  wantu1 = lsame( jobu1, 'Y' )
359  wantu2 = lsame( jobu2, 'Y' )
360  wantv1t = lsame( jobv1t, 'Y' )
361  wantv2t = lsame( jobv2t, 'Y' )
362  colmajor = .NOT. lsame( trans, 'T' )
363  defaultsigns = .NOT. lsame( signs, 'O' )
364  lquery = lwork .EQ. -1
365  IF( m .LT. 0 ) THEN
366  info = -7
367  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
368  info = -8
369  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
370  info = -9
371  ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
372  info = -11
373  ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
374  info = -11
375  ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
376  info = -13
377  ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
378  info = -13
379  ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
380  info = -15
381  ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
382  info = -15
383  ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
384  info = -17
385  ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
386  info = -17
387  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
388  info = -20
389  ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
390  info = -22
391  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
392  info = -24
393  ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
394  info = -26
395  END IF
396 *
397 * Work with transpose if convenient
398 *
399  IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
400  IF( colmajor ) THEN
401  transt = 'T'
402  ELSE
403  transt = 'N'
404  END IF
405  IF( defaultsigns ) THEN
406  signst = 'O'
407  ELSE
408  signst = 'D'
409  END IF
410  CALL dorcsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
411  $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
412  $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
413  $ u2, ldu2, work, lwork, iwork, info )
414  RETURN
415  END IF
416 *
417 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
418 * convenient
419 *
420  IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
421  IF( defaultsigns ) THEN
422  signst = 'O'
423  ELSE
424  signst = 'D'
425  END IF
426  CALL dorcsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
427  $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
428  $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
429  $ ldv1t, work, lwork, iwork, info )
430  RETURN
431  END IF
432 *
433 * Compute workspace
434 *
435  IF( info .EQ. 0 ) THEN
436 *
437  iphi = 2
438  itaup1 = iphi + max( 1, q - 1 )
439  itaup2 = itaup1 + max( 1, p )
440  itauq1 = itaup2 + max( 1, m - p )
441  itauq2 = itauq1 + max( 1, q )
442  iorgqr = itauq2 + max( 1, m - q )
443  CALL dorgqr( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
444  $ childinfo )
445  lorgqrworkopt = int( work(1) )
446  lorgqrworkmin = max( 1, m - q )
447  iorglq = itauq2 + max( 1, m - q )
448  CALL dorglq( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
449  $ childinfo )
450  lorglqworkopt = int( work(1) )
451  lorglqworkmin = max( 1, m - q )
452  iorbdb = itauq2 + max( 1, m - q )
453  CALL dorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
454  $ x21, ldx21, x22, ldx22, theta, v1t, u1, u2, v1t,
455  $ v2t, work, -1, childinfo )
456  lorbdbworkopt = int( work(1) )
457  lorbdbworkmin = lorbdbworkopt
458  ib11d = itauq2 + max( 1, m - q )
459  ib11e = ib11d + max( 1, q )
460  ib12d = ib11e + max( 1, q - 1 )
461  ib12e = ib12d + max( 1, q )
462  ib21d = ib12e + max( 1, q - 1 )
463  ib21e = ib21d + max( 1, q )
464  ib22d = ib21e + max( 1, q - 1 )
465  ib22e = ib22d + max( 1, q )
466  ibbcsd = ib22e + max( 1, q - 1 )
467  CALL dbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
468  $ theta, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
469  $ ldv2t, u1, u1, u1, u1, u1, u1, u1, u1, work, -1,
470  $ childinfo )
471  lbbcsdworkopt = int( work(1) )
472  lbbcsdworkmin = lbbcsdworkopt
473  lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
474  $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkopt ) - 1
475  lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
476  $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkmin ) - 1
477  work(1) = max(lworkopt,lworkmin)
478 *
479  IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
480  info = -22
481  ELSE
482  lorgqrwork = lwork - iorgqr + 1
483  lorglqwork = lwork - iorglq + 1
484  lorbdbwork = lwork - iorbdb + 1
485  lbbcsdwork = lwork - ibbcsd + 1
486  END IF
487  END IF
488 *
489 * Abort if any illegal arguments
490 *
491  IF( info .NE. 0 ) THEN
492  CALL xerbla( 'DORCSD', -info )
493  RETURN
494  ELSE IF( lquery ) THEN
495  RETURN
496  END IF
497 *
498 * Transform to bidiagonal block form
499 *
500  CALL dorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
501  $ ldx21, x22, ldx22, theta, work(iphi), work(itaup1),
502  $ work(itaup2), work(itauq1), work(itauq2),
503  $ work(iorbdb), lorbdbwork, childinfo )
504 *
505 * Accumulate Householder reflectors
506 *
507  IF( colmajor ) THEN
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  $ lorgqrwork, info)
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), lorgqrwork, info )
517  END IF
518  IF( wantv1t .AND. q .GT. 0 ) THEN
519  CALL dlacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
520  $ ldv1t )
521  v1t(1, 1) = one
522  DO j = 2, q
523  v1t(1,j) = zero
524  v1t(j,1) = zero
525  END DO
526  CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
527  $ work(iorglq), lorglqwork, info )
528  END IF
529  IF( wantv2t .AND. m-q .GT. 0 ) THEN
530  CALL dlacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
531  IF (m-p .GT. q) Then
532  CALL dlacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
533  $ v2t(p+1,p+1), ldv2t )
534  END IF
535  IF (m .GT. q) THEN
536  CALL dorglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
537  $ work(iorglq), lorglqwork, info )
538  END IF
539  END IF
540  ELSE
541  IF( wantu1 .AND. p .GT. 0 ) THEN
542  CALL dlacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
543  CALL dorglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
544  $ lorglqwork, info)
545  END IF
546  IF( wantu2 .AND. m-p .GT. 0 ) THEN
547  CALL dlacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
548  CALL dorglq( m-p, m-p, q, u2, ldu2, work(itaup2),
549  $ work(iorglq), lorglqwork, info )
550  END IF
551  IF( wantv1t .AND. q .GT. 0 ) THEN
552  CALL dlacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
553  $ ldv1t )
554  v1t(1, 1) = one
555  DO j = 2, q
556  v1t(1,j) = zero
557  v1t(j,1) = zero
558  END DO
559  CALL dorgqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
560  $ work(iorgqr), lorgqrwork, info )
561  END IF
562  IF( wantv2t .AND. m-q .GT. 0 ) THEN
563  CALL dlacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
564  CALL dlacpy( 'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
565  $ v2t(p+1,p+1), ldv2t )
566  CALL dorgqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
567  $ work(iorgqr), lorgqrwork, info )
568  END IF
569  END IF
570 *
571 * Compute the CSD of the matrix in bidiagonal-block form
572 *
573  CALL dbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
574  $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
575  $ ldv2t, work(ib11d), work(ib11e), work(ib12d),
576  $ work(ib12e), work(ib21d), work(ib21e), work(ib22d),
577  $ work(ib22e), work(ibbcsd), lbbcsdwork, info )
578 *
579 * Permute rows and columns to place identity submatrices in top-
580 * left corner of (1,1)-block and/or bottom-right corner of (1,2)-
581 * block and/or bottom-right corner of (2,1)-block and/or top-left
582 * corner of (2,2)-block
583 *
584  IF( q .GT. 0 .AND. wantu2 ) THEN
585  DO i = 1, q
586  iwork(i) = m - p - q + i
587  END DO
588  DO i = q + 1, m - p
589  iwork(i) = i - q
590  END DO
591  IF( colmajor ) THEN
592  CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
593  ELSE
594  CALL dlapmr( .false., m-p, m-p, u2, ldu2, iwork )
595  END IF
596  END IF
597  IF( m .GT. 0 .AND. wantv2t ) THEN
598  DO i = 1, p
599  iwork(i) = m - p - q + i
600  END DO
601  DO i = p + 1, m - q
602  iwork(i) = i - p
603  END DO
604  IF( .NOT. colmajor ) THEN
605  CALL dlapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
606  ELSE
607  CALL dlapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
608  END IF
609  END IF
610 *
611  RETURN
612 *
613 * End DORCSD
614 *
615  END
616