LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
sorcsd.f
Go to the documentation of this file.
1 *> \brief \b SORCSD
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SORCSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sorcsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sorcsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sorcsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * RECURSIVE SUBROUTINE SORCSD( 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 * REAL THETA( * )
35 * REAL 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 *> SORCSD 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (LDU1,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 REAL array, dimension (LDU2,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 REAL array, dimension (LDV1T,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 REAL array, dimension (LDV2T,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 REAL 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: SBBCSD 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 *> \ingroup realOTHERcomputational
293 *
294 * =====================================================================
295  RECURSIVE SUBROUTINE sorcsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
296  $ SIGNS, M, P, Q, X11, LDX11, X12,
297  $ LDX12, X21, LDX21, X22, LDX22, THETA,
298  $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
299  $ LDV2T, WORK, LWORK, IWORK, INFO )
300 *
301 * -- LAPACK computational routine --
302 * -- LAPACK is a software package provided by Univ. of Tennessee, --
303 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
304 *
305 * .. Scalar Arguments ..
306  CHARACTER jobu1, jobu2, jobv1t, jobv2t, signs, trans
307  INTEGER info, ldu1, ldu2, ldv1t, ldv2t, ldx11, ldx12,
308  $ ldx21, ldx22, lwork, m, p, q
309 * ..
310 * .. Array Arguments ..
311  INTEGER iwork( * )
312  REAL theta( * )
313  REAL u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
314  $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
315  $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
316  $ * )
317 * ..
318 *
319 * ===================================================================
320 *
321 * .. Parameters ..
322  REAL one, zero
323  PARAMETER ( one = 1.0e+0,
324  $ zero = 0.0e+0 )
325 * ..
326 * .. Local Arrays ..
327  REAL dummy(1)
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 sbbcsd, slacpy, slapmr, slapmt,
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 sorcsd( 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 sorcsd( 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 sorgqr( m-q, m-q, m-q, dummy, max(1,m-q), dummy, work, -1,
444  $ childinfo )
445  lorgqrworkopt = int( work(1) )
446  lorgqrworkmin = max( 1, m - q )
447  iorglq = itauq2 + max( 1, m - q )
448  CALL sorglq( m-q, m-q, m-q, dummy, max(1,m-q), dummy, work, -1,
449  $ childinfo )
450  lorglqworkopt = int( work(1) )
451  lorglqworkmin = max( 1, m - q )
452  iorbdb = itauq2 + max( 1, m - q )
453  CALL sorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
454  $ x21, ldx21, x22, ldx22, dummy, dummy, dummy, dummy, dummy,
455  $ dummy,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 sbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
468  $ dummy, dummy, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
469  $ ldv2t, dummy, dummy, dummy, dummy, dummy, dummy,
470  $ dummy, dummy, work, -1, 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( 'SORCSD', -info )
493  RETURN
494  ELSE IF( lquery ) THEN
495  RETURN
496  END IF
497 *
498 * Transform to bidiagonal block form
499 *
500  CALL sorbdb( 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 slacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
510  CALL sorgqr( 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 slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
515  CALL sorgqr( 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 slacpy( '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 sorglq( 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 slacpy( 'U', p, m-q, x12, ldx12, v2t, ldv2t )
531  CALL slacpy( 'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
532  $ v2t(p+1,p+1), ldv2t )
533  CALL sorglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
534  $ work(iorglq), lorglqwork, info )
535  END IF
536  ELSE
537  IF( wantu1 .AND. p .GT. 0 ) THEN
538  CALL slacpy( 'U', q, p, x11, ldx11, u1, ldu1 )
539  CALL sorglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
540  $ lorglqwork, info)
541  END IF
542  IF( wantu2 .AND. m-p .GT. 0 ) THEN
543  CALL slacpy( 'U', q, m-p, x21, ldx21, u2, ldu2 )
544  CALL sorglq( m-p, m-p, q, u2, ldu2, work(itaup2),
545  $ work(iorglq), lorglqwork, info )
546  END IF
547  IF( wantv1t .AND. q .GT. 0 ) THEN
548  CALL slacpy( 'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
549  $ ldv1t )
550  v1t(1, 1) = one
551  DO j = 2, q
552  v1t(1,j) = zero
553  v1t(j,1) = zero
554  END DO
555  CALL sorgqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
556  $ work(iorgqr), lorgqrwork, info )
557  END IF
558  IF( wantv2t .AND. m-q .GT. 0 ) THEN
559  CALL slacpy( 'L', m-q, p, x12, ldx12, v2t, ldv2t )
560  CALL slacpy( 'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
561  $ v2t(p+1,p+1), ldv2t )
562  CALL sorgqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
563  $ work(iorgqr), lorgqrwork, info )
564  END IF
565  END IF
566 *
567 * Compute the CSD of the matrix in bidiagonal-block form
568 *
569  CALL sbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
570  $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
571  $ ldv2t, work(ib11d), work(ib11e), work(ib12d),
572  $ work(ib12e), work(ib21d), work(ib21e), work(ib22d),
573  $ work(ib22e), work(ibbcsd), lbbcsdwork, info )
574 *
575 * Permute rows and columns to place identity submatrices in top-
576 * left corner of (1,1)-block and/or bottom-right corner of (1,2)-
577 * block and/or bottom-right corner of (2,1)-block and/or top-left
578 * corner of (2,2)-block
579 *
580  IF( q .GT. 0 .AND. wantu2 ) THEN
581  DO i = 1, q
582  iwork(i) = m - p - q + i
583  END DO
584  DO i = q + 1, m - p
585  iwork(i) = i - q
586  END DO
587  IF( colmajor ) THEN
588  CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
589  ELSE
590  CALL slapmr( .false., m-p, m-p, u2, ldu2, iwork )
591  END IF
592  END IF
593  IF( m .GT. 0 .AND. wantv2t ) THEN
594  DO i = 1, p
595  iwork(i) = m - p - q + i
596  END DO
597  DO i = p + 1, m - q
598  iwork(i) = i - p
599  END DO
600  IF( .NOT. colmajor ) THEN
601  CALL slapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
602  ELSE
603  CALL slapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
604  END IF
605  END IF
606 *
607  RETURN
608 *
609 * End SORCSD
610 *
611  END
612 
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:103
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine slapmr(FORWRD, M, N, X, LDX, K)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: slapmr.f:104
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: slapmt.f:104
recursive subroutine sorcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, WORK, LWORK, IWORK, INFO)
SORCSD
Definition: sorcsd.f:300
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ
Definition: sorglq.f:127
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:128
subroutine sbbcsd(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)
SBBCSD
Definition: sbbcsd.f:332
subroutine sorbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
SORBDB
Definition: sorbdb.f:287