LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dorcsd2by1.f
Go to the documentation of this file.
1 *> \brief \b DORCSD2BY1
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DORCSD2BY1 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dorcsd2by1.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dorcsd2by1.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dorcsd2by1.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DORCSD2BY1( 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 * DOUBLE PRECISION THETA(*)
32 * DOUBLE PRECISION 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 *> Purpose:
43 *> ========
44 *>
45 *> DORCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
46 *> orthonormal columns that has been partitioned into a 2-by-1 block
47 *> structure:
48 *>
49 *> [ I 0 0 ]
50 *> [ 0 C 0 ]
51 *> [ X11 ] [ U1 | ] [ 0 0 0 ]
52 *> X = [-----] = [---------] [----------] V1**T .
53 *> [ X21 ] [ | U2 ] [ 0 0 0 ]
54 *> [ 0 S 0 ]
55 *> [ 0 0 I ]
56 *>
57 *> X11 is P-by-Q. The orthogonal matrices U1, U2, V1, and V2 are P-by-P,
58 *> (M-P)-by-(M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. C and S are
59 *> R-by-R nonnegative diagonal matrices satisfying C^2 + S^2 = I, in
60 *> which R = MIN(P,M-P,Q,M-Q).
61 *>
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] M
89 *> \verbatim
90 *> M is INTEGER
91 *> The number of rows and columns in X.
92 *> \endverbatim
93 *>
94 *> \param[in] P
95 *> \verbatim
96 *> P is INTEGER
97 *> The number of rows in X11 and X12. 0 <= P <= M.
98 *> \endverbatim
99 *>
100 *> \param[in] Q
101 *> \verbatim
102 *> Q is INTEGER
103 *> The number of columns in X11 and X21. 0 <= Q <= M.
104 *> \endverbatim
105 *>
106 *> \param[in,out] X11
107 *> \verbatim
108 *> X11 is DOUBLE PRECISION array, dimension (LDX11,Q)
109 *> On entry, part of the orthogonal matrix whose CSD is
110 *> desired.
111 *> \endverbatim
112 *>
113 *> \param[in] LDX11
114 *> \verbatim
115 *> LDX11 is INTEGER
116 *> The leading dimension of X11. LDX11 >= MAX(1,P).
117 *> \endverbatim
118 *>
119 *> \param[in,out] X21
120 *> \verbatim
121 *> X21 is DOUBLE PRECISION array, dimension (LDX21,Q)
122 *> On entry, part of the orthogonal matrix whose CSD is
123 *> desired.
124 *> \endverbatim
125 *>
126 *> \param[in] LDX21
127 *> \verbatim
128 *> LDX21 is INTEGER
129 *> The leading dimension of X21. LDX21 >= MAX(1,M-P).
130 *> \endverbatim
131 *>
132 *> \param[out] THETA
133 *> \verbatim
134 *> THETA is DOUBLE PRECISION array, dimension (R), in which R =
135 *> MIN(P,M-P,Q,M-Q).
136 *> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
137 *> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
138 *> \endverbatim
139 *>
140 *> \param[out] U1
141 *> \verbatim
142 *> U1 is DOUBLE PRECISION array, dimension (P)
143 *> If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
144 *> \endverbatim
145 *>
146 *> \param[in] LDU1
147 *> \verbatim
148 *> LDU1 is INTEGER
149 *> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
150 *> MAX(1,P).
151 *> \endverbatim
152 *>
153 *> \param[out] U2
154 *> \verbatim
155 *> U2 is DOUBLE PRECISION array, dimension (M-P)
156 *> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
157 *> matrix U2.
158 *> \endverbatim
159 *>
160 *> \param[in] LDU2
161 *> \verbatim
162 *> LDU2 is INTEGER
163 *> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
164 *> MAX(1,M-P).
165 *> \endverbatim
166 *>
167 *> \param[out] V1T
168 *> \verbatim
169 *> V1T is DOUBLE PRECISION array, dimension (Q)
170 *> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
171 *> matrix V1**T.
172 *> \endverbatim
173 *>
174 *> \param[in] LDV1T
175 *> \verbatim
176 *> LDV1T is INTEGER
177 *> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
178 *> MAX(1,Q).
179 *> \endverbatim
180 *>
181 *> \param[out] WORK
182 *> \verbatim
183 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
184 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
185 *> If INFO > 0 on exit, WORK(2:R) contains the values PHI(1),
186 *> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
187 *> define the matrix in intermediate bidiagonal-block form
188 *> remaining after nonconvergence. INFO specifies the number
189 *> of nonzero PHI's.
190 *> \endverbatim
191 *>
192 *> \param[in] LWORK
193 *> \verbatim
194 *> LWORK is INTEGER
195 *> The dimension of the array WORK.
196 *> \endverbatim
197 *>
198 *> If LWORK = -1, then a workspace query is assumed; the routine
199 *> only calculates the optimal size of the WORK array, returns
200 *> this value as the first entry of the work array, and no error
201 *> message related to LWORK is issued by XERBLA.
202 *> \param[out] IWORK
203 *> \verbatim
204 *> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
205 *> \endverbatim
206 *>
207 *> \param[out] INFO
208 *> \verbatim
209 *> INFO is INTEGER
210 *> = 0: successful exit.
211 *> < 0: if INFO = -i, the i-th argument had an illegal value.
212 *> > 0: DBBCSD did not converge. See the description of WORK
213 *> above for details.
214 *> \endverbatim
215 *
216 * Authors:
217 * ========
218 *
219 *> \author Univ. of Tennessee
220 *> \author Univ. of California Berkeley
221 *> \author Univ. of Colorado Denver
222 *> \author NAG Ltd.
223 *
224 *> \date July 2012
225 *
226 *> \ingroup doubleOTHERcomputational
227 *
228 *> \par References:
229 * ================
230 *>
231 *> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
232 *> Algorithms, 50(1):33-65, 2009.
233 *> \endverbatim
234 *>
235 * =====================================================================
236  SUBROUTINE dorcsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
237  $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
238  $ ldv1t, work, lwork, iwork, info )
239 *
240 * -- LAPACK computational routine (3.5.0) --
241 * -- LAPACK is a software package provided by Univ. of Tennessee, --
242 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
243 * July 2012
244 *
245 * .. Scalar Arguments ..
246  CHARACTER jobu1, jobu2, jobv1t
247  INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
248  $ m, p, q
249 * ..
250 * .. Array Arguments ..
251  DOUBLE PRECISION theta(*)
252  DOUBLE PRECISION u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
253  $ x11(ldx11,*), x21(ldx21,*)
254  INTEGER iwork(*)
255 * ..
256 *
257 * =====================================================================
258 *
259 * .. Parameters ..
260  DOUBLE PRECISION one, zero
261  parameter( one = 1.0d0, zero = 0.0d0 )
262 * ..
263 * .. Local Scalars ..
264  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
265  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
266  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
267  $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
268  $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
269  $ lworkmin, lworkopt, r
270  LOGICAL lquery, wantu1, wantu2, wantv1t
271 * ..
272 * .. External Subroutines ..
273  EXTERNAL dbbcsd, dcopy, dlacpy, dlapmr, dlapmt, dorbdb1,
275  $ xerbla
276 * ..
277 * .. External Functions ..
278  LOGICAL lsame
279  EXTERNAL lsame
280 * ..
281 * .. Intrinsic Function ..
282  INTRINSIC int, max, min
283 * ..
284 * .. Executable Statements ..
285 *
286 * Test input arguments
287 *
288  info = 0
289  wantu1 = lsame( jobu1, 'Y' )
290  wantu2 = lsame( jobu2, 'Y' )
291  wantv1t = lsame( jobv1t, 'Y' )
292  lquery = lwork .EQ. -1
293 *
294  IF( m .LT. 0 ) THEN
295  info = -4
296  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
297  info = -5
298  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
299  info = -6
300  ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
301  info = -8
302  ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
303  info = -10
304  ELSE IF( wantu1 .AND. ldu1 .LT. p ) THEN
305  info = -13
306  ELSE IF( wantu2 .AND. ldu2 .LT. m - p ) THEN
307  info = -15
308  ELSE IF( wantv1t .AND. ldv1t .LT. q ) THEN
309  info = -17
310  END IF
311 *
312  r = min( p, m-p, q, m-q )
313 *
314 * Compute workspace
315 *
316 * WORK layout:
317 * |-------------------------------------------------------|
318 * | LWORKOPT (1) |
319 * |-------------------------------------------------------|
320 * | PHI (MAX(1,R-1)) |
321 * |-------------------------------------------------------|
322 * | TAUP1 (MAX(1,P)) | B11D (R) |
323 * | TAUP2 (MAX(1,M-P)) | B11E (R-1) |
324 * | TAUQ1 (MAX(1,Q)) | B12D (R) |
325 * |-----------------------------------------| B12E (R-1) |
326 * | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R) |
327 * | | | | B21E (R-1) |
328 * | | | | B22D (R) |
329 * | | | | B22E (R-1) |
330 * | | | | DBBCSD WORK |
331 * |-------------------------------------------------------|
332 *
333  IF( info .EQ. 0 ) THEN
334  iphi = 2
335  ib11d = iphi + max( 1, r-1 )
336  ib11e = ib11d + max( 1, r )
337  ib12d = ib11e + max( 1, r - 1 )
338  ib12e = ib12d + max( 1, r )
339  ib21d = ib12e + max( 1, r - 1 )
340  ib21e = ib21d + max( 1, r )
341  ib22d = ib21e + max( 1, r - 1 )
342  ib22e = ib22d + max( 1, r )
343  ibbcsd = ib22e + max( 1, r - 1 )
344  itaup1 = iphi + max( 1, r-1 )
345  itaup2 = itaup1 + max( 1, p )
346  itauq1 = itaup2 + max( 1, m-p )
347  iorbdb = itauq1 + max( 1, q )
348  iorgqr = itauq1 + max( 1, q )
349  iorglq = itauq1 + max( 1, q )
350  IF( r .EQ. q ) THEN
351  CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
352  $ 0, 0, work, -1, childinfo )
353  lorbdb = int( work(1) )
354  IF( p .GE. m-p ) THEN
355  CALL dorgqr( p, p, q, u1, ldu1, 0, work(1), -1,
356  $ childinfo )
357  lorgqrmin = max( 1, p )
358  lorgqropt = int( work(1) )
359  ELSE
360  CALL dorgqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
361  $ childinfo )
362  lorgqrmin = max( 1, m-p )
363  lorgqropt = int( work(1) )
364  END IF
365  CALL dorglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
366  $ 0, work(1), -1, childinfo )
367  lorglqmin = max( 1, q-1 )
368  lorglqopt = int( work(1) )
369  CALL dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
370  $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
371  $ 0, 0, 0, 0, 0, 0, work(1), -1, childinfo )
372  lbbcsd = int( work(1) )
373  ELSE IF( r .EQ. p ) THEN
374  CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
375  $ 0, 0, work(1), -1, childinfo )
376  lorbdb = int( work(1) )
377  IF( p-1 .GE. m-p ) THEN
378  CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
379  $ -1, childinfo )
380  lorgqrmin = max( 1, p-1 )
381  lorgqropt = int( work(1) )
382  ELSE
383  CALL dorgqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
384  $ childinfo )
385  lorgqrmin = max( 1, m-p )
386  lorgqropt = int( work(1) )
387  END IF
388  CALL dorglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
389  $ childinfo )
390  lorglqmin = max( 1, q )
391  lorglqopt = int( work(1) )
392  CALL dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
393  $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
394  $ 0, 0, 0, 0, 0, 0, work(1), -1, childinfo )
395  lbbcsd = int( work(1) )
396  ELSE IF( r .EQ. m-p ) THEN
397  CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
398  $ 0, 0, work(1), -1, childinfo )
399  lorbdb = int( work(1) )
400  IF( p .GE. m-p-1 ) THEN
401  CALL dorgqr( p, p, q, u1, ldu1, 0, work(1), -1,
402  $ childinfo )
403  lorgqrmin = max( 1, p )
404  lorgqropt = int( work(1) )
405  ELSE
406  CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
407  $ work(1), -1, childinfo )
408  lorgqrmin = max( 1, m-p-1 )
409  lorgqropt = int( work(1) )
410  END IF
411  CALL dorglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
412  $ childinfo )
413  lorglqmin = max( 1, q )
414  lorglqopt = int( work(1) )
415  CALL dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
416  $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
417  $ 0, 0, 0, 0, 0, 0, 0, 0, work(1), -1,
418  $ childinfo )
419  lbbcsd = int( work(1) )
420  ELSE
421  CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
422  $ 0, 0, 0, work(1), -1, childinfo )
423  lorbdb = m + int( work(1) )
424  IF( p .GE. m-p ) THEN
425  CALL dorgqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
426  $ childinfo )
427  lorgqrmin = max( 1, p )
428  lorgqropt = int( work(1) )
429  ELSE
430  CALL dorgqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
431  $ childinfo )
432  lorgqrmin = max( 1, m-p )
433  lorgqropt = int( work(1) )
434  END IF
435  CALL dorglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
436  $ childinfo )
437  lorglqmin = max( 1, q )
438  lorglqopt = int( work(1) )
439  CALL dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
440  $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
441  $ 0, 0, 0, 0, 0, 0, 0, 0, work(1), -1,
442  $ childinfo )
443  lbbcsd = int( work(1) )
444  END IF
445  lworkmin = max( iorbdb+lorbdb-1,
446  $ iorgqr+lorgqrmin-1,
447  $ iorglq+lorglqmin-1,
448  $ ibbcsd+lbbcsd-1 )
449  lworkopt = max( iorbdb+lorbdb-1,
450  $ iorgqr+lorgqropt-1,
451  $ iorglq+lorglqopt-1,
452  $ ibbcsd+lbbcsd-1 )
453  work(1) = lworkopt
454  IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
455  info = -19
456  END IF
457  END IF
458  IF( info .NE. 0 ) THEN
459  CALL xerbla( 'DORCSD2BY1', -info )
460  RETURN
461  ELSE IF( lquery ) THEN
462  RETURN
463  END IF
464  lorgqr = lwork-iorgqr+1
465  lorglq = lwork-iorglq+1
466 *
467 * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
468 * in which R = MIN(P,M-P,Q,M-Q)
469 *
470  IF( r .EQ. q ) THEN
471 *
472 * Case 1: R = Q
473 *
474 * Simultaneously bidiagonalize X11 and X21
475 *
476  CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
477  $ work(iphi), work(itaup1), work(itaup2),
478  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
479 *
480 * Accumulate Householder reflectors
481 *
482  IF( wantu1 .AND. p .GT. 0 ) THEN
483  CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
484  CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
485  $ lorgqr, childinfo )
486  END IF
487  IF( wantu2 .AND. m-p .GT. 0 ) THEN
488  CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
489  CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
490  $ work(iorgqr), lorgqr, childinfo )
491  END IF
492  IF( wantv1t .AND. q .GT. 0 ) THEN
493  v1t(1,1) = one
494  DO j = 2, q
495  v1t(1,j) = zero
496  v1t(j,1) = zero
497  END DO
498  CALL dlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
499  $ ldv1t )
500  CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
501  $ work(iorglq), lorglq, childinfo )
502  END IF
503 *
504 * Simultaneously diagonalize X11 and X21.
505 *
506  CALL dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
507  $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
508  $ work(ib11d), work(ib11e), work(ib12d),
509  $ work(ib12e), work(ib21d), work(ib21e),
510  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
511  $ childinfo )
512 *
513 * Permute rows and columns to place zero submatrices in
514 * preferred positions
515 *
516  IF( q .GT. 0 .AND. wantu2 ) THEN
517  DO i = 1, q
518  iwork(i) = m - p - q + i
519  END DO
520  DO i = q + 1, m - p
521  iwork(i) = i - q
522  END DO
523  CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
524  END IF
525  ELSE IF( r .EQ. p ) THEN
526 *
527 * Case 2: R = P
528 *
529 * Simultaneously bidiagonalize X11 and X21
530 *
531  CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
532  $ work(iphi), work(itaup1), work(itaup2),
533  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
534 *
535 * Accumulate Householder reflectors
536 *
537  IF( wantu1 .AND. p .GT. 0 ) THEN
538  u1(1,1) = one
539  DO j = 2, p
540  u1(1,j) = zero
541  u1(j,1) = zero
542  END DO
543  CALL dlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
544  CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
545  $ work(iorgqr), lorgqr, childinfo )
546  END IF
547  IF( wantu2 .AND. m-p .GT. 0 ) THEN
548  CALL dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
549  CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
550  $ work(iorgqr), lorgqr, childinfo )
551  END IF
552  IF( wantv1t .AND. q .GT. 0 ) THEN
553  CALL dlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
554  CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
555  $ work(iorglq), lorglq, childinfo )
556  END IF
557 *
558 * Simultaneously diagonalize X11 and X21.
559 *
560  CALL dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
561  $ work(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
562  $ work(ib11d), work(ib11e), work(ib12d),
563  $ work(ib12e), work(ib21d), work(ib21e),
564  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
565  $ childinfo )
566 *
567 * Permute rows and columns to place identity submatrices in
568 * preferred positions
569 *
570  IF( q .GT. 0 .AND. wantu2 ) THEN
571  DO i = 1, q
572  iwork(i) = m - p - q + i
573  END DO
574  DO i = q + 1, m - p
575  iwork(i) = i - q
576  END DO
577  CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
578  END IF
579  ELSE IF( r .EQ. m-p ) THEN
580 *
581 * Case 3: R = M-P
582 *
583 * Simultaneously bidiagonalize X11 and X21
584 *
585  CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
586  $ work(iphi), work(itaup1), work(itaup2),
587  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
588 *
589 * Accumulate Householder reflectors
590 *
591  IF( wantu1 .AND. p .GT. 0 ) THEN
592  CALL dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
593  CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
594  $ lorgqr, childinfo )
595  END IF
596  IF( wantu2 .AND. m-p .GT. 0 ) THEN
597  u2(1,1) = one
598  DO j = 2, m-p
599  u2(1,j) = zero
600  u2(j,1) = zero
601  END DO
602  CALL dlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
603  $ ldu2 )
604  CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
605  $ work(itaup2), work(iorgqr), lorgqr, childinfo )
606  END IF
607  IF( wantv1t .AND. q .GT. 0 ) THEN
608  CALL dlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
609  CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
610  $ work(iorglq), lorglq, childinfo )
611  END IF
612 *
613 * Simultaneously diagonalize X11 and X21.
614 *
615  CALL dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
616  $ theta, work(iphi), 0, 1, v1t, ldv1t, u2, ldu2, u1,
617  $ ldu1, work(ib11d), work(ib11e), work(ib12d),
618  $ work(ib12e), work(ib21d), work(ib21e),
619  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
620  $ childinfo )
621 *
622 * Permute rows and columns to place identity submatrices in
623 * preferred positions
624 *
625  IF( q .GT. r ) THEN
626  DO i = 1, r
627  iwork(i) = q - r + i
628  END DO
629  DO i = r + 1, q
630  iwork(i) = i - r
631  END DO
632  IF( wantu1 ) THEN
633  CALL dlapmt( .false., p, q, u1, ldu1, iwork )
634  END IF
635  IF( wantv1t ) THEN
636  CALL dlapmr( .false., q, q, v1t, ldv1t, iwork )
637  END IF
638  END IF
639  ELSE
640 *
641 * Case 4: R = M-Q
642 *
643 * Simultaneously bidiagonalize X11 and X21
644 *
645  CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
646  $ work(iphi), work(itaup1), work(itaup2),
647  $ work(itauq1), work(iorbdb), work(iorbdb+m),
648  $ lorbdb-m, childinfo )
649 *
650 * Accumulate Householder reflectors
651 *
652  IF( wantu1 .AND. p .GT. 0 ) THEN
653  CALL dcopy( p, work(iorbdb), 1, u1, 1 )
654  DO j = 2, p
655  u1(1,j) = zero
656  END DO
657  CALL dlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
658  $ ldu1 )
659  CALL dorgqr( p, p, m-q, u1, ldu1, work(itaup1),
660  $ work(iorgqr), lorgqr, childinfo )
661  END IF
662  IF( wantu2 .AND. m-p .GT. 0 ) THEN
663  CALL dcopy( m-p, work(iorbdb+p), 1, u2, 1 )
664  DO j = 2, m-p
665  u2(1,j) = zero
666  END DO
667  CALL dlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
668  $ ldu2 )
669  CALL dorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
670  $ work(iorgqr), lorgqr, childinfo )
671  END IF
672  IF( wantv1t .AND. q .GT. 0 ) THEN
673  CALL dlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
674  CALL dlacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
675  $ v1t(m-q+1,m-q+1), ldv1t )
676  CALL dlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
677  $ v1t(p+1,p+1), ldv1t )
678  CALL dorglq( q, q, q, v1t, ldv1t, work(itauq1),
679  $ work(iorglq), lorglq, childinfo )
680  END IF
681 *
682 * Simultaneously diagonalize X11 and X21.
683 *
684  CALL dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
685  $ theta, work(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
686  $ ldv1t, work(ib11d), work(ib11e), work(ib12d),
687  $ work(ib12e), work(ib21d), work(ib21e),
688  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
689  $ childinfo )
690 *
691 * Permute rows and columns to place identity submatrices in
692 * preferred positions
693 *
694  IF( p .GT. r ) THEN
695  DO i = 1, r
696  iwork(i) = p - r + i
697  END DO
698  DO i = r + 1, p
699  iwork(i) = i - r
700  END DO
701  IF( wantu1 ) THEN
702  CALL dlapmt( .false., p, p, u1, ldu1, iwork )
703  END IF
704  IF( wantv1t ) THEN
705  CALL dlapmr( .false., p, q, v1t, ldv1t, iwork )
706  END IF
707  END IF
708  END IF
709 *
710  RETURN
711 *
712 * End of DORCSD2BY1
713 *
714  END
715 
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
Definition: dorglq.f:128
subroutine dorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB1
Definition: dorbdb1.f:203
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:52
subroutine dorbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
DORBDB4
Definition: dorbdb4.f:212
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:61
subroutine dorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
DORCSD2BY1
Definition: dorcsd2by1.f:236
subroutine dorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB2
Definition: dorbdb2.f:202
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:54
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
Definition: dorgqr.f:129
subroutine dorbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB3
Definition: dorbdb3.f:201
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
Definition: xerbla-fortran:9
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:104
subroutine dbbcsd(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)
DBBCSD
Definition: dbbcsd.f:330
subroutine dlapmt(FORWRD, M, N, X, LDX, K)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: dlapmt.f:105
subroutine dlapmr(FORWRD, M, N, X, LDX, K)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: dlapmr.f:105