LAPACK  3.4.2 LAPACK: Linear Algebra PACKage
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
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 2011
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.4.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 2011
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) ) .OR.
372  \$ ( .NOT.colmajor .AND. ldx11 .LT. max(1,q) ) ) THEN
373  info = -11
374  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
375  info = -20
376  ELSE IF( wantu2 .AND. ldu2 .LT. m-p ) THEN
377  info = -22
378  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
379  info = -24
380  ELSE IF( wantv2t .AND. ldv2t .LT. m-q ) THEN
381  info = -26
382  END IF
383 *
384 * Work with transpose if convenient
385 *
386  IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) ) THEN
387  IF( colmajor ) THEN
388  transt = 'T'
389  ELSE
390  transt = 'N'
391  END IF
392  IF( defaultsigns ) THEN
393  signst = 'O'
394  ELSE
395  signst = 'D'
396  END IF
397  CALL dorcsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
398  \$ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
399  \$ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
400  \$ u2, ldu2, work, lwork, iwork, info )
401  return
402  END IF
403 *
404 * Work with permutation [ 0 I; I 0 ] * X * [ 0 I; I 0 ] if
405 * convenient
406 *
407  IF( info .EQ. 0 .AND. m-q .LT. q ) THEN
408  IF( defaultsigns ) THEN
409  signst = 'O'
410  ELSE
411  signst = 'D'
412  END IF
413  CALL dorcsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
414  \$ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
415  \$ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
416  \$ ldv1t, work, lwork, iwork, info )
417  return
418  END IF
419 *
420 * Compute workspace
421 *
422  IF( info .EQ. 0 ) THEN
423 *
424  iphi = 2
425  itaup1 = iphi + max( 1, q - 1 )
426  itaup2 = itaup1 + max( 1, p )
427  itauq1 = itaup2 + max( 1, m - p )
428  itauq2 = itauq1 + max( 1, q )
429  iorgqr = itauq2 + max( 1, m - q )
430  CALL dorgqr( m-q, m-q, m-q, 0, max(1,m-q), 0, work, -1,
431  \$ childinfo )
432  lorgqrworkopt = int( work(1) )
433  lorgqrworkmin = max( 1, m - q )
434  iorglq = itauq2 + max( 1, m - q )
435  CALL dorglq( m-q, m-q, m-q, 0, max(1,m-q), 0, work, -1,
436  \$ childinfo )
437  lorglqworkopt = int( work(1) )
438  lorglqworkmin = max( 1, m - q )
439  iorbdb = itauq2 + max( 1, m - q )
440  CALL dorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
441  \$ x21, ldx21, x22, ldx22, 0, 0, 0, 0, 0, 0, work,
442  \$ -1, childinfo )
443  lorbdbworkopt = int( work(1) )
444  lorbdbworkmin = lorbdbworkopt
445  ib11d = itauq2 + max( 1, m - q )
446  ib11e = ib11d + max( 1, q )
447  ib12d = ib11e + max( 1, q - 1 )
448  ib12e = ib12d + max( 1, q )
449  ib21d = ib12e + max( 1, q - 1 )
450  ib21e = ib21d + max( 1, q )
451  ib22d = ib21e + max( 1, q - 1 )
452  ib22e = ib22d + max( 1, q )
453  ibbcsd = ib22e + max( 1, q - 1 )
454  CALL dbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, 0,
455  \$ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t, ldv2t, 0,
456  \$ 0, 0, 0, 0, 0, 0, 0, work, -1, childinfo )
457  lbbcsdworkopt = int( work(1) )
458  lbbcsdworkmin = lbbcsdworkopt
459  lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
460  \$ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkopt ) - 1
461  lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
462  \$ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkmin ) - 1
463  work(1) = max(lworkopt,lworkmin)
464 *
465  IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
466  info = -22
467  ELSE
468  lorgqrwork = lwork - iorgqr + 1
469  lorglqwork = lwork - iorglq + 1
470  lorbdbwork = lwork - iorbdb + 1
471  lbbcsdwork = lwork - ibbcsd + 1
472  END IF
473  END IF
474 *
475 * Abort if any illegal arguments
476 *
477  IF( info .NE. 0 ) THEN
478  CALL xerbla( 'DORCSD', -info )
479  return
480  ELSE IF( lquery ) THEN
481  return
482  END IF
483 *
484 * Transform to bidiagonal block form
485 *
486  CALL dorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
487  \$ ldx21, x22, ldx22, theta, work(iphi), work(itaup1),
488  \$ work(itaup2), work(itauq1), work(itauq2),
489  \$ work(iorbdb), lorbdbwork, childinfo )
490 *
491 * Accumulate Householder reflectors
492 *
493  IF( colmajor ) THEN
494  IF( wantu1 .AND. p .GT. 0 ) THEN
495  CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
496  CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
497  \$ lorgqrwork, info)
498  END IF
499  IF( wantu2 .AND. m-p .GT. 0 ) THEN
500  CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
501  CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
502  \$ work(iorgqr), lorgqrwork, info )
503  END IF
504  IF( wantv1t .AND. q .GT. 0 ) THEN
505  CALL dlacpy( 'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
506  \$ ldv1t )
507  v1t(1, 1) = one
508  DO j = 2, q
509  v1t(1,j) = zero
510  v1t(j,1) = zero
511  END DO
512  CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
513  \$ work(iorglq), lorglqwork, info )
514  END IF
515  IF( wantv2t .AND. m-q .GT. 0 ) THEN
516  CALL dlacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
517  CALL dlacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
518  \$ v2t(p+1,p+1), ldv2t )
519  CALL dorglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
520  \$ work(iorglq), lorglqwork, info )
521  END IF
522  ELSE
523  IF( wantu1 .AND. p .GT. 0 ) THEN
524  CALL dlacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
525  CALL dorglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
526  \$ lorglqwork, info)
527  END IF
528  IF( wantu2 .AND. m-p .GT. 0 ) THEN
529  CALL dlacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
530  CALL dorglq( m-p, m-p, q, u2, ldu2, work(itaup2),
531  \$ work(iorglq), lorglqwork, info )
532  END IF
533  IF( wantv1t .AND. q .GT. 0 ) THEN
534  CALL dlacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
535  \$ ldv1t )
536  v1t(1, 1) = one
537  DO j = 2, q
538  v1t(1,j) = zero
539  v1t(j,1) = zero
540  END DO
541  CALL dorgqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
542  \$ work(iorgqr), lorgqrwork, info )
543  END IF
544  IF( wantv2t .AND. m-q .GT. 0 ) THEN
545  CALL dlacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
546  CALL dlacpy( 'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
547  \$ v2t(p+1,p+1), ldv2t )
548  CALL dorgqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
549  \$ work(iorgqr), lorgqrwork, info )
550  END IF
551  END IF
552 *
553 * Compute the CSD of the matrix in bidiagonal-block form
554 *
555  CALL dbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
556  \$ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
557  \$ ldv2t, work(ib11d), work(ib11e), work(ib12d),
558  \$ work(ib12e), work(ib21d), work(ib21e), work(ib22d),
559  \$ work(ib22e), work(ibbcsd), lbbcsdwork, info )
560 *
561 * Permute rows and columns to place identity submatrices in top-
562 * left corner of (1,1)-block and/or bottom-right corner of (1,2)-
563 * block and/or bottom-right corner of (2,1)-block and/or top-left
564 * corner of (2,2)-block
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  IF( colmajor ) THEN
574  CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
575  ELSE
576  CALL dlapmr( .false., m-p, m-p, u2, ldu2, iwork )
577  END IF
578  END IF
579  IF( m .GT. 0 .AND. wantv2t ) THEN
580  DO i = 1, p
581  iwork(i) = m - p - q + i
582  END DO
583  DO i = p + 1, m - q
584  iwork(i) = i - p
585  END DO
586  IF( .NOT. colmajor ) THEN
587  CALL dlapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
588  ELSE
589  CALL dlapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
590  END IF
591  END IF
592 *
593  return
594 *
595 * End DORCSD
596 *
597  END
598