LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
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 
subroutine slapmr(FORWRD, M, N, X, LDX, K)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: slapmr.f:106
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB2
Definition: sorbdb2.f:203
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:334
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53
subroutine sorbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
SORBDB4
Definition: sorbdb4.f:216
subroutine sorbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB3
Definition: sorbdb3.f:204
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: slapmt.f:106
subroutine sorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB1
Definition: sorbdb1.f:205
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:130
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ
Definition: sorglq.f:129
subroutine sorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
SORCSD2BY1
Definition: sorcsd2by1.f:235