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