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