OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
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*> [ I1 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 I2]
57*>
58*> X11 is P-by-Q. The unitary matrices U1, U2, and V1 are P-by-P,
59*> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
60*> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
61*> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a
62*> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
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 in X.
93*> \endverbatim
94*>
95*> \param[in] P
96*> \verbatim
97*> P is INTEGER
98*> The number of rows in X11. 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 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 COMPLEX*16 array, dimension (LDX21,Q)
122*> On entry, part of the unitary matrix whose CSD is desired.
123*> \endverbatim
124*>
125*> \param[in] LDX21
126*> \verbatim
127*> LDX21 is INTEGER
128*> The leading dimension of X21. LDX21 >= MAX(1,M-P).
129*> \endverbatim
130*>
131*> \param[out] THETA
132*> \verbatim
133*> THETA is DOUBLE PRECISION array, dimension (R), in which R =
134*> MIN(P,M-P,Q,M-Q).
135*> C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
136*> S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
137*> \endverbatim
138*>
139*> \param[out] U1
140*> \verbatim
141*> U1 is COMPLEX*16 array, dimension (P)
142*> If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
143*> \endverbatim
144*>
145*> \param[in] LDU1
146*> \verbatim
147*> LDU1 is INTEGER
148*> The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
149*> MAX(1,P).
150*> \endverbatim
151*>
152*> \param[out] U2
153*> \verbatim
154*> U2 is COMPLEX*16 array, dimension (M-P)
155*> If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
156*> matrix U2.
157*> \endverbatim
158*>
159*> \param[in] LDU2
160*> \verbatim
161*> LDU2 is INTEGER
162*> The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
163*> MAX(1,M-P).
164*> \endverbatim
165*>
166*> \param[out] V1T
167*> \verbatim
168*> V1T is COMPLEX*16 array, dimension (Q)
169*> If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
170*> matrix V1**T.
171*> \endverbatim
172*>
173*> \param[in] LDV1T
174*> \verbatim
175*> LDV1T is INTEGER
176*> The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
177*> MAX(1,Q).
178*> \endverbatim
179*>
180*> \param[out] WORK
181*> \verbatim
182*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
183*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
184*> \endverbatim
185*>
186*> \param[in] LWORK
187*> \verbatim
188*> LWORK is INTEGER
189*> The dimension of the array WORK.
190*>
191*> If LWORK = -1, then a workspace query is assumed; the routine
192*> only calculates the optimal size of the WORK and RWORK
193*> arrays, returns this value as the first entry of the WORK
194*> and RWORK array, respectively, and no error message related
195*> to LWORK or LRWORK is issued by XERBLA.
196*> \endverbatim
197*>
198*> \param[out] RWORK
199*> \verbatim
200*> RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
201*> On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
202*> If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
203*> ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
204*> define the matrix in intermediate bidiagonal-block form
205*> remaining after nonconvergence. INFO specifies the number
206*> of nonzero PHI's.
207*> \endverbatim
208*>
209*> \param[in] LRWORK
210*> \verbatim
211*> LRWORK is INTEGER
212*> The dimension of the array RWORK.
213*>
214*> If LRWORK=-1, then a workspace query is assumed; the routine
215*> only calculates the optimal size of the WORK and RWORK
216*> arrays, returns this value as the first entry of the WORK
217*> and RWORK array, respectively, and no error message related
218*> to LWORK or LRWORK is issued by XERBLA.
219*> \endverbatim
220*
221*> \param[out] IWORK
222*> \verbatim
223*> IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
224*> \endverbatim
225*>
226*> \param[out] INFO
227*> \verbatim
228*> INFO is INTEGER
229*> = 0: successful exit.
230*> < 0: if INFO = -i, the i-th argument had an illegal value.
231*> > 0: ZBBCSD did not converge. See the description of WORK
232*> above for details.
233*> \endverbatim
234*
235*> \par References:
236* ================
237*>
238*> [1] Brian D. Sutton. Computing the complete CS decomposition. Numer.
239*> Algorithms, 50(1):33-65, 2009.
240*
241* Authors:
242* ========
243*
244*> \author Univ. of Tennessee
245*> \author Univ. of California Berkeley
246*> \author Univ. of Colorado Denver
247*> \author NAG Ltd.
248*
249*> \ingroup complex16OTHERcomputational
250*
251* =====================================================================
252 SUBROUTINE zuncsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
253 $ X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T,
254 $ LDV1T, WORK, LWORK, RWORK, LRWORK, IWORK,
255 $ INFO )
256*
257* -- LAPACK computational routine --
258* -- LAPACK is a software package provided by Univ. of Tennessee, --
259* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
260*
261* .. Scalar Arguments ..
262 CHARACTER JOBU1, JOBU2, JOBV1T
263 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
264 $ M, P, Q
265 INTEGER LRWORK, LRWORKMIN, LRWORKOPT
266* ..
267* .. Array Arguments ..
268 DOUBLE PRECISION RWORK(*)
269 DOUBLE PRECISION THETA(*)
270 COMPLEX*16 U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
271 $ x11(ldx11,*), x21(ldx21,*)
272 INTEGER IWORK(*)
273* ..
274*
275* =====================================================================
276*
277* .. Parameters ..
278 COMPLEX*16 ONE, ZERO
279 PARAMETER ( ONE = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
280* ..
281* .. Local Scalars ..
282 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
283 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
284 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
285 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
286 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
287 $ lworkmin, lworkopt, r
288 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
289* ..
290* .. Local Arrays ..
291 DOUBLE PRECISION DUM( 1 )
292 COMPLEX*16 CDUM( 1, 1 )
293* ..
294* .. External Subroutines ..
295 EXTERNAL zbbcsd, zcopy, zlacpy, zlapmr, zlapmt, zunbdb1,
296 $ zunbdb2, zunbdb3, zunbdb4, zunglq, zungqr,
297 $ xerbla
298* ..
299* .. External Functions ..
300 LOGICAL LSAME
301 EXTERNAL LSAME
302* ..
303* .. Intrinsic Function ..
304 INTRINSIC int, max, min
305* ..
306* .. Executable Statements ..
307*
308* Test input arguments
309*
310 info = 0
311 wantu1 = lsame( jobu1, 'Y' )
312 wantu2 = lsame( jobu2, 'Y' )
313 wantv1t = lsame( jobv1t, 'Y' )
314 lquery = ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 )
315*
316 IF( m .LT. 0 ) THEN
317 info = -4
318 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
319 info = -5
320 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
321 info = -6
322 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
323 info = -8
324 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
325 info = -10
326 ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) ) THEN
327 info = -13
328 ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) ) THEN
329 info = -15
330 ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) ) THEN
331 info = -17
332 END IF
333*
334 r = min( p, m-p, q, m-q )
335*
336* Compute workspace
337*
338* WORK layout:
339* |-----------------------------------------|
340* | LWORKOPT (1) |
341* |-----------------------------------------|
342* | TAUP1 (MAX(1,P)) |
343* | TAUP2 (MAX(1,M-P)) |
344* | TAUQ1 (MAX(1,Q)) |
345* |-----------------------------------------|
346* | ZUNBDB WORK | ZUNGQR WORK | ZUNGLQ WORK |
347* | | | |
348* | | | |
349* | | | |
350* | | | |
351* |-----------------------------------------|
352* RWORK layout:
353* |------------------|
354* | LRWORKOPT (1) |
355* |------------------|
356* | PHI (MAX(1,R-1)) |
357* |------------------|
358* | B11D (R) |
359* | B11E (R-1) |
360* | B12D (R) |
361* | B12E (R-1) |
362* | B21D (R) |
363* | B21E (R-1) |
364* | B22D (R) |
365* | B22E (R-1) |
366* | ZBBCSD RWORK |
367* |------------------|
368*
369 IF( info .EQ. 0 ) THEN
370 iphi = 2
371 ib11d = iphi + max( 1, r-1 )
372 ib11e = ib11d + max( 1, r )
373 ib12d = ib11e + max( 1, r - 1 )
374 ib12e = ib12d + max( 1, r )
375 ib21d = ib12e + max( 1, r - 1 )
376 ib21e = ib21d + max( 1, r )
377 ib22d = ib21e + max( 1, r - 1 )
378 ib22e = ib22d + max( 1, r )
379 ibbcsd = ib22e + max( 1, r - 1 )
380 itaup1 = 2
381 itaup2 = itaup1 + max( 1, p )
382 itauq1 = itaup2 + max( 1, m-p )
383 iorbdb = itauq1 + max( 1, q )
384 iorgqr = itauq1 + max( 1, q )
385 iorglq = itauq1 + max( 1, q )
386 lorgqrmin = 1
387 lorgqropt = 1
388 lorglqmin = 1
389 lorglqopt = 1
390 IF( r .EQ. q ) THEN
391 CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
392 $ cdum, cdum, cdum, work, -1, childinfo )
393 lorbdb = int( work(1) )
394 IF( wantu1 .AND. p .GT. 0 ) THEN
395 CALL zungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
396 $ childinfo )
397 lorgqrmin = max( lorgqrmin, p )
398 lorgqropt = max( lorgqropt, int( work(1) ) )
399 ENDIF
400 IF( wantu2 .AND. m-p .GT. 0 ) THEN
401 CALL zungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
402 $ childinfo )
403 lorgqrmin = max( lorgqrmin, m-p )
404 lorgqropt = max( lorgqropt, int( work(1) ) )
405 END IF
406 IF( wantv1t .AND. q .GT. 0 ) THEN
407 CALL zunglq( q-1, q-1, q-1, v1t, ldv1t,
408 $ cdum, work(1), -1, childinfo )
409 lorglqmin = max( lorglqmin, q-1 )
410 lorglqopt = max( lorglqopt, int( work(1) ) )
411 END IF
412 CALL zbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
413 $ dum, u1, ldu1, u2, ldu2, v1t, ldv1t, cdum, 1,
414 $ dum, dum, dum, dum, dum, dum, dum, dum,
415 $ rwork(1), -1, childinfo )
416 lbbcsd = int( rwork(1) )
417 ELSE IF( r .EQ. p ) THEN
418 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
419 $ cdum, cdum, cdum, work(1), -1, childinfo )
420 lorbdb = int( work(1) )
421 IF( wantu1 .AND. p .GT. 0 ) THEN
422 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, cdum, work(1),
423 $ -1, childinfo )
424 lorgqrmin = max( lorgqrmin, p-1 )
425 lorgqropt = max( lorgqropt, int( work(1) ) )
426 END IF
427 IF( wantu2 .AND. m-p .GT. 0 ) THEN
428 CALL zungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
429 $ childinfo )
430 lorgqrmin = max( lorgqrmin, m-p )
431 lorgqropt = max( lorgqropt, int( work(1) ) )
432 END IF
433 IF( wantv1t .AND. q .GT. 0 ) THEN
434 CALL zunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
435 $ childinfo )
436 lorglqmin = max( lorglqmin, q )
437 lorglqopt = max( lorglqopt, int( work(1) ) )
438 END IF
439 CALL zbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
440 $ dum, v1t, ldv1t, cdum, 1, u1, ldu1, u2, ldu2,
441 $ dum, dum, dum, dum, dum, dum, dum, dum,
442 $ rwork(1), -1, childinfo )
443 lbbcsd = int( rwork(1) )
444 ELSE IF( r .EQ. m-p ) THEN
445 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
446 $ cdum, cdum, cdum, work(1), -1, childinfo )
447 lorbdb = int( work(1) )
448 IF( wantu1 .AND. p .GT. 0 ) THEN
449 CALL zungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
450 $ childinfo )
451 lorgqrmin = max( lorgqrmin, p )
452 lorgqropt = max( lorgqropt, int( work(1) ) )
453 END IF
454 IF( wantu2 .AND. m-p .GT. 0 ) THEN
455 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, cdum,
456 $ work(1), -1, childinfo )
457 lorgqrmin = max( lorgqrmin, m-p-1 )
458 lorgqropt = max( lorgqropt, int( work(1) ) )
459 END IF
460 IF( wantv1t .AND. q .GT. 0 ) THEN
461 CALL zunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
462 $ childinfo )
463 lorglqmin = max( lorglqmin, q )
464 lorglqopt = max( lorglqopt, int( work(1) ) )
465 END IF
466 CALL zbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
467 $ theta, dum, cdum, 1, v1t, ldv1t, u2, ldu2, u1,
468 $ ldu1, dum, dum, dum, dum, dum, dum, dum, dum,
469 $ rwork(1), -1, childinfo )
470 lbbcsd = int( rwork(1) )
471 ELSE
472 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
473 $ cdum, cdum, cdum, cdum, work(1), -1, childinfo
474 $ )
475 lorbdb = m + int( work(1) )
476 IF( wantu1 .AND. p .GT. 0 ) THEN
477 CALL zungqr( p, p, m-q, u1, ldu1, cdum, work(1), -1,
478 $ childinfo )
479 lorgqrmin = max( lorgqrmin, p )
480 lorgqropt = max( lorgqropt, int( work(1) ) )
481 END IF
482 IF( wantu2 .AND. m-p .GT. 0 ) THEN
483 CALL zungqr( m-p, m-p, m-q, u2, ldu2, cdum, work(1), -1,
484 $ childinfo )
485 lorgqrmin = max( lorgqrmin, m-p )
486 lorgqropt = max( lorgqropt, int( work(1) ) )
487 END IF
488 IF( wantv1t .AND. q .GT. 0 ) THEN
489 CALL zunglq( q, q, q, v1t, ldv1t, cdum, work(1), -1,
490 $ childinfo )
491 lorglqmin = max( lorglqmin, q )
492 lorglqopt = max( lorglqopt, int( work(1) ) )
493 END IF
494 CALL zbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
495 $ theta, dum, u2, ldu2, u1, ldu1, cdum, 1, v1t,
496 $ ldv1t, dum, dum, dum, dum, dum, dum, dum, dum,
497 $ rwork(1), -1, childinfo )
498 lbbcsd = int( rwork(1) )
499 END IF
500 lrworkmin = ibbcsd+lbbcsd-1
501 lrworkopt = lrworkmin
502 rwork(1) = lrworkopt
503 lworkmin = max( iorbdb+lorbdb-1,
504 $ iorgqr+lorgqrmin-1,
505 $ iorglq+lorglqmin-1 )
506 lworkopt = max( iorbdb+lorbdb-1,
507 $ iorgqr+lorgqropt-1,
508 $ iorglq+lorglqopt-1 )
509 work(1) = lworkopt
510 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
511 info = -19
512 END IF
513 IF( lrwork .LT. lrworkmin .AND. .NOT.lquery ) THEN
514 info = -21
515 END IF
516 END IF
517 IF( info .NE. 0 ) THEN
518 CALL xerbla( 'ZUNCSD2BY1', -info )
519 RETURN
520 ELSE IF( lquery ) THEN
521 RETURN
522 END IF
523 lorgqr = lwork-iorgqr+1
524 lorglq = lwork-iorglq+1
525*
526* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
527* in which R = MIN(P,M-P,Q,M-Q)
528*
529 IF( r .EQ. q ) THEN
530*
531* Case 1: R = Q
532*
533* Simultaneously bidiagonalize X11 and X21
534*
535 CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
536 $ rwork(iphi), work(itaup1), work(itaup2),
537 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
538*
539* Accumulate Householder reflectors
540*
541 IF( wantu1 .AND. p .GT. 0 ) THEN
542 CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
543 CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
544 $ lorgqr, childinfo )
545 END IF
546 IF( wantu2 .AND. m-p .GT. 0 ) THEN
547 CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
548 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
549 $ work(iorgqr), lorgqr, childinfo )
550 END IF
551 IF( wantv1t .AND. q .GT. 0 ) THEN
552 v1t(1,1) = one
553 DO j = 2, q
554 v1t(1,j) = zero
555 v1t(j,1) = zero
556 END DO
557 CALL zlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
558 $ ldv1t )
559 CALL zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
560 $ work(iorglq), lorglq, childinfo )
561 END IF
562*
563* Simultaneously diagonalize X11 and X21.
564*
565 CALL zbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
566 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
567 $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
568 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
569 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
570 $ lrwork-ibbcsd+1, childinfo )
571*
572* Permute rows and columns to place zero submatrices in
573* preferred positions
574*
575 IF( q .GT. 0 .AND. wantu2 ) THEN
576 DO i = 1, q
577 iwork(i) = m - p - q + i
578 END DO
579 DO i = q + 1, m - p
580 iwork(i) = i - q
581 END DO
582 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
583 END IF
584 ELSE IF( r .EQ. p ) THEN
585*
586* Case 2: R = P
587*
588* Simultaneously bidiagonalize X11 and X21
589*
590 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
591 $ rwork(iphi), work(itaup1), work(itaup2),
592 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
593*
594* Accumulate Householder reflectors
595*
596 IF( wantu1 .AND. p .GT. 0 ) THEN
597 u1(1,1) = one
598 DO j = 2, p
599 u1(1,j) = zero
600 u1(j,1) = zero
601 END DO
602 CALL zlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
603 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
604 $ work(iorgqr), lorgqr, childinfo )
605 END IF
606 IF( wantu2 .AND. m-p .GT. 0 ) THEN
607 CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
608 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
609 $ work(iorgqr), lorgqr, childinfo )
610 END IF
611 IF( wantv1t .AND. q .GT. 0 ) THEN
612 CALL zlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
613 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
614 $ work(iorglq), lorglq, childinfo )
615 END IF
616*
617* Simultaneously diagonalize X11 and X21.
618*
619 CALL zbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
620 $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
621 $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
622 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
623 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
624 $ childinfo )
625*
626* Permute rows and columns to place identity submatrices in
627* preferred positions
628*
629 IF( q .GT. 0 .AND. wantu2 ) THEN
630 DO i = 1, q
631 iwork(i) = m - p - q + i
632 END DO
633 DO i = q + 1, m - p
634 iwork(i) = i - q
635 END DO
636 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
637 END IF
638 ELSE IF( r .EQ. m-p ) THEN
639*
640* Case 3: R = M-P
641*
642* Simultaneously bidiagonalize X11 and X21
643*
644 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
645 $ rwork(iphi), work(itaup1), work(itaup2),
646 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
647*
648* Accumulate Householder reflectors
649*
650 IF( wantu1 .AND. p .GT. 0 ) THEN
651 CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
652 CALL zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
653 $ lorgqr, childinfo )
654 END IF
655 IF( wantu2 .AND. m-p .GT. 0 ) THEN
656 u2(1,1) = one
657 DO j = 2, m-p
658 u2(1,j) = zero
659 u2(j,1) = zero
660 END DO
661 CALL zlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
662 $ ldu2 )
663 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
664 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
665 END IF
666 IF( wantv1t .AND. q .GT. 0 ) THEN
667 CALL zlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
668 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
669 $ work(iorglq), lorglq, childinfo )
670 END IF
671*
672* Simultaneously diagonalize X11 and X21.
673*
674 CALL zbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
675 $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
676 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
677 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
678 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
679 $ rwork(ibbcsd), lbbcsd, childinfo )
680*
681* Permute rows and columns to place identity submatrices in
682* preferred positions
683*
684 IF( q .GT. r ) THEN
685 DO i = 1, r
686 iwork(i) = q - r + i
687 END DO
688 DO i = r + 1, q
689 iwork(i) = i - r
690 END DO
691 IF( wantu1 ) THEN
692 CALL zlapmt( .false., p, q, u1, ldu1, iwork )
693 END IF
694 IF( wantv1t ) THEN
695 CALL zlapmr( .false., q, q, v1t, ldv1t, iwork )
696 END IF
697 END IF
698 ELSE
699*
700* Case 4: R = M-Q
701*
702* Simultaneously bidiagonalize X11 and X21
703*
704 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
705 $ rwork(iphi), work(itaup1), work(itaup2),
706 $ work(itauq1), work(iorbdb), work(iorbdb+m),
707 $ lorbdb-m, childinfo )
708*
709* Accumulate Householder reflectors
710*
711 IF( wantu2 .AND. m-p .GT. 0 ) THEN
712 CALL zcopy( m-p, work(iorbdb+p), 1, u2, 1 )
713 END IF
714 IF( wantu1 .AND. p .GT. 0 ) THEN
715 CALL zcopy( p, work(iorbdb), 1, u1, 1 )
716 DO j = 2, p
717 u1(1,j) = zero
718 END DO
719 CALL zlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
720 $ ldu1 )
721 CALL zungqr( p, p, m-q, u1, ldu1, work(itaup1),
722 $ work(iorgqr), lorgqr, childinfo )
723 END IF
724 IF( wantu2 .AND. m-p .GT. 0 ) THEN
725 DO j = 2, m-p
726 u2(1,j) = zero
727 END DO
728 CALL zlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
729 $ ldu2 )
730 CALL zungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
731 $ work(iorgqr), lorgqr, childinfo )
732 END IF
733 IF( wantv1t .AND. q .GT. 0 ) THEN
734 CALL zlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
735 CALL zlacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
736 $ v1t(m-q+1,m-q+1), ldv1t )
737 CALL zlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
738 $ v1t(p+1,p+1), ldv1t )
739 CALL zunglq( q, q, q, v1t, ldv1t, work(itauq1),
740 $ work(iorglq), lorglq, childinfo )
741 END IF
742*
743* Simultaneously diagonalize X11 and X21.
744*
745 CALL zbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
746 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
747 $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
748 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
749 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
750 $ rwork(ibbcsd), lbbcsd, childinfo )
751*
752* Permute rows and columns to place identity submatrices in
753* preferred positions
754*
755 IF( p .GT. r ) THEN
756 DO i = 1, r
757 iwork(i) = p - r + i
758 END DO
759 DO i = r + 1, p
760 iwork(i) = i - r
761 END DO
762 IF( wantu1 ) THEN
763 CALL zlapmt( .false., p, p, u1, ldu1, iwork )
764 END IF
765 IF( wantv1t ) THEN
766 CALL zlapmr( .false., p, q, v1t, ldv1t, iwork )
767 END IF
768 END IF
769 END IF
770*
771 RETURN
772*
773* End of ZUNCSD2BY1
774*
775 END
776
xerbla
subroutine xerbla(srname, info)
XERBLA
Definition xerbla.f:60
zlacpy
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
zlapmt
subroutine zlapmt(forwrd, m, n, x, ldx, k)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition zlapmt.f:104
zlapmr
subroutine zlapmr(forwrd, m, n, x, ldx, k)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition zlapmr.f:104
zuncsd2by1
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:256
zunbdb3
subroutine zunbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB3
Definition zunbdb3.f:201
zunbdb4
subroutine zunbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
ZUNBDB4
Definition zunbdb4.f:213
zunglq
subroutine zunglq(m, n, k, a, lda, tau, work, lwork, info)
ZUNGLQ
Definition zunglq.f:127
zungqr
subroutine zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
Definition zungqr.f:128
zunbdb1
subroutine zunbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB1
Definition zunbdb1.f:203
zbbcsd
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:332
zunbdb2
subroutine zunbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB2
Definition zunbdb2.f:201
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21