LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ zuncsd2by1()

subroutine zuncsd2by1 ( character jobu1,
character jobu2,
character jobv1t,
integer m,
integer p,
integer q,
complex*16, dimension(ldx11,*) x11,
integer ldx11,
complex*16, dimension(ldx21,*) x21,
integer ldx21,
double precision, dimension(*) theta,
complex*16, dimension(ldu1,*) u1,
integer ldu1,
complex*16, dimension(ldu2,*) u2,
integer ldu2,
complex*16, dimension(ldv1t,*) v1t,
integer ldv1t,
complex*16, dimension(*) work,
integer lwork,
double precision, dimension(*) rwork,
integer lrwork,
integer, dimension(*) iwork,
integer info )

ZUNCSD2BY1

Download ZUNCSD2BY1 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> ZUNCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
!> orthonormal columns that has been partitioned into a 2-by-1 block
!> structure:
!>
!>                                [  I1 0  0 ]
!>                                [  0  C  0 ]
!>          [ X11 ]   [ U1 |    ] [  0  0  0 ]
!>      X = [-----] = [---------] [----------] V1**T .
!>          [ X21 ]   [    | U2 ] [  0  0  0 ]
!>                                [  0  S  0 ]
!>                                [  0  0  I2]
!>
!> X11 is P-by-Q. The unitary matrices U1, U2, and V1 are P-by-P,
!> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
!> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
!> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a
!> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
!>
Parameters
[in]JOBU1
!>          JOBU1 is CHARACTER
!>          = 'Y':      U1 is computed;
!>          otherwise:  U1 is not computed.
!>
[in]JOBU2
!>          JOBU2 is CHARACTER
!>          = 'Y':      U2 is computed;
!>          otherwise:  U2 is not computed.
!>
[in]JOBV1T
!>          JOBV1T is CHARACTER
!>          = 'Y':      V1T is computed;
!>          otherwise:  V1T is not computed.
!>
[in]M
!>          M is INTEGER
!>          The number of rows in X.
!>
[in]P
!>          P is INTEGER
!>          The number of rows in X11. 0 <= P <= M.
!>
[in]Q
!>          Q is INTEGER
!>          The number of columns in X11 and X21. 0 <= Q <= M.
!>
[in,out]X11
!>          X11 is COMPLEX*16 array, dimension (LDX11,Q)
!>          On entry, part of the unitary matrix whose CSD is desired.
!>
[in]LDX11
!>          LDX11 is INTEGER
!>          The leading dimension of X11. LDX11 >= MAX(1,P).
!>
[in,out]X21
!>          X21 is COMPLEX*16 array, dimension (LDX21,Q)
!>          On entry, part of the unitary matrix whose CSD is desired.
!>
[in]LDX21
!>          LDX21 is INTEGER
!>          The leading dimension of X21. LDX21 >= MAX(1,M-P).
!>
[out]THETA
!>          THETA is DOUBLE PRECISION array, dimension (R), in which R =
!>          MIN(P,M-P,Q,M-Q).
!>          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
!>          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
!>
[out]U1
!>          U1 is COMPLEX*16 array, dimension (P)
!>          If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
!>
[in]LDU1
!>          LDU1 is INTEGER
!>          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
!>          MAX(1,P).
!>
[out]U2
!>          U2 is COMPLEX*16 array, dimension (M-P)
!>          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
!>          matrix U2.
!>
[in]LDU2
!>          LDU2 is INTEGER
!>          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
!>          MAX(1,M-P).
!>
[out]V1T
!>          V1T is COMPLEX*16 array, dimension (Q)
!>          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
!>          matrix V1**T.
!>
[in]LDV1T
!>          LDV1T is INTEGER
!>          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
!>          MAX(1,Q).
!>
[out]WORK
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK and RWORK
!>          arrays, returns this value as the first entry of the WORK
!>          and RWORK array, respectively, and no error message related
!>          to LWORK or LRWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is DOUBLE PRECISION array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>          If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
!>          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
!>          define the matrix in intermediate bidiagonal-block form
!>          remaining after nonconvergence. INFO specifies the number
!>          of nonzero PHI's.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of the array RWORK.
!>
!>          If LRWORK=-1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK and RWORK
!>          arrays, returns this value as the first entry of the WORK
!>          and RWORK array, respectively, and no error message related
!>          to LWORK or LRWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  ZBBCSD did not converge. See the description of WORK
!>                above for details.
!>
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 250 of file zuncsd2by1.f.

255*
256* -- LAPACK computational routine --
257* -- LAPACK is a software package provided by Univ. of Tennessee, --
258* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
259*
260* .. Scalar Arguments ..
261 CHARACTER JOBU1, JOBU2, JOBV1T
262 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
263 $ M, P, Q
264 INTEGER LRWORK, LRWORKMIN, LRWORKOPT
265* ..
266* .. Array Arguments ..
267 DOUBLE PRECISION RWORK(*)
268 DOUBLE PRECISION THETA(*)
269 COMPLEX*16 U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
270 $ X11(LDX11,*), X21(LDX21,*)
271 INTEGER IWORK(*)
272* ..
273*
274* =====================================================================
275*
276* .. Parameters ..
277 COMPLEX*16 ONE, ZERO
278 parameter( one = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
279* ..
280* .. Local Scalars ..
281 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
282 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
283 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
284 $ J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN,
285 $ LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT,
286 $ LWORKMIN, LWORKOPT, R
287 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
288* ..
289* .. Local Arrays ..
290 DOUBLE PRECISION DUM( 1 )
291 COMPLEX*16 CDUM( 1, 1 )
292* ..
293* .. External Subroutines ..
294 EXTERNAL zbbcsd, zcopy, zlacpy, zlapmr, zlapmt,
295 $ 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,
392 $ dum,
393 $ cdum, cdum, cdum, work, -1, childinfo )
394 lorbdb = int( work(1) )
395 IF( wantu1 .AND. p .GT. 0 ) THEN
396 CALL zungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
397 $ childinfo )
398 lorgqrmin = max( lorgqrmin, p )
399 lorgqropt = max( lorgqropt, int( work(1) ) )
400 ENDIF
401 IF( wantu2 .AND. m-p .GT. 0 ) THEN
402 CALL zungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
403 $ childinfo )
404 lorgqrmin = max( lorgqrmin, m-p )
405 lorgqropt = max( lorgqropt, int( work(1) ) )
406 END IF
407 IF( wantv1t .AND. q .GT. 0 ) THEN
408 CALL zunglq( q-1, q-1, q-1, v1t, ldv1t,
409 $ cdum, work(1), -1, childinfo )
410 lorglqmin = max( lorglqmin, q-1 )
411 lorglqopt = max( lorglqopt, int( work(1) ) )
412 END IF
413 CALL zbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q,
414 $ theta,
415 $ dum, u1, ldu1, u2, ldu2, v1t, ldv1t, cdum, 1,
416 $ dum, dum, dum, dum, dum, dum, dum, dum,
417 $ rwork(1), -1, childinfo )
418 lbbcsd = int( rwork(1) )
419 ELSE IF( r .EQ. p ) THEN
420 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
421 $ dum,
422 $ cdum, cdum, cdum, work(1), -1, childinfo )
423 lorbdb = int( work(1) )
424 IF( wantu1 .AND. p .GT. 0 ) THEN
425 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, cdum,
426 $ work(1),
427 $ -1, childinfo )
428 lorgqrmin = max( lorgqrmin, p-1 )
429 lorgqropt = max( lorgqropt, int( work(1) ) )
430 END IF
431 IF( wantu2 .AND. m-p .GT. 0 ) THEN
432 CALL zungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
433 $ childinfo )
434 lorgqrmin = max( lorgqrmin, m-p )
435 lorgqropt = max( lorgqropt, int( work(1) ) )
436 END IF
437 IF( wantv1t .AND. q .GT. 0 ) THEN
438 CALL zunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
439 $ childinfo )
440 lorglqmin = max( lorglqmin, q )
441 lorglqopt = max( lorglqopt, int( work(1) ) )
442 END IF
443 CALL zbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p,
444 $ theta,
445 $ dum, v1t, ldv1t, cdum, 1, u1, ldu1, u2, ldu2,
446 $ dum, dum, dum, dum, dum, dum, dum, dum,
447 $ rwork(1), -1, childinfo )
448 lbbcsd = int( rwork(1) )
449 ELSE IF( r .EQ. m-p ) THEN
450 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
451 $ dum,
452 $ cdum, cdum, cdum, work(1), -1, childinfo )
453 lorbdb = int( work(1) )
454 IF( wantu1 .AND. p .GT. 0 ) THEN
455 CALL zungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
456 $ childinfo )
457 lorgqrmin = max( lorgqrmin, p )
458 lorgqropt = max( lorgqropt, int( work(1) ) )
459 END IF
460 IF( wantu2 .AND. m-p .GT. 0 ) THEN
461 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, cdum,
462 $ work(1), -1, childinfo )
463 lorgqrmin = max( lorgqrmin, m-p-1 )
464 lorgqropt = max( lorgqropt, int( work(1) ) )
465 END IF
466 IF( wantv1t .AND. q .GT. 0 ) THEN
467 CALL zunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
468 $ childinfo )
469 lorglqmin = max( lorglqmin, q )
470 lorglqopt = max( lorglqopt, int( work(1) ) )
471 END IF
472 CALL zbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
473 $ theta, dum, cdum, 1, v1t, ldv1t, u2, ldu2, u1,
474 $ ldu1, dum, dum, dum, dum, dum, dum, dum, dum,
475 $ rwork(1), -1, childinfo )
476 lbbcsd = int( rwork(1) )
477 ELSE
478 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
479 $ dum,
480 $ cdum, cdum, cdum, cdum, work(1), -1, childinfo
481 $ )
482 lorbdb = m + int( work(1) )
483 IF( wantu1 .AND. p .GT. 0 ) THEN
484 CALL zungqr( p, p, m-q, u1, ldu1, cdum, work(1), -1,
485 $ childinfo )
486 lorgqrmin = max( lorgqrmin, p )
487 lorgqropt = max( lorgqropt, int( work(1) ) )
488 END IF
489 IF( wantu2 .AND. m-p .GT. 0 ) THEN
490 CALL zungqr( m-p, m-p, m-q, u2, ldu2, cdum, work(1),
491 $ -1,
492 $ childinfo )
493 lorgqrmin = max( lorgqrmin, m-p )
494 lorgqropt = max( lorgqropt, int( work(1) ) )
495 END IF
496 IF( wantv1t .AND. q .GT. 0 ) THEN
497 CALL zunglq( q, q, q, v1t, ldv1t, cdum, work(1), -1,
498 $ childinfo )
499 lorglqmin = max( lorglqmin, q )
500 lorglqopt = max( lorglqopt, int( work(1) ) )
501 END IF
502 CALL zbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
503 $ theta, dum, u2, ldu2, u1, ldu1, cdum, 1, v1t,
504 $ ldv1t, dum, dum, dum, dum, dum, dum, dum, dum,
505 $ rwork(1), -1, childinfo )
506 lbbcsd = int( rwork(1) )
507 END IF
508 lrworkmin = ibbcsd+lbbcsd-1
509 lrworkopt = lrworkmin
510 rwork(1) = lrworkopt
511 lworkmin = max( iorbdb+lorbdb-1,
512 $ iorgqr+lorgqrmin-1,
513 $ iorglq+lorglqmin-1 )
514 lworkopt = max( iorbdb+lorbdb-1,
515 $ iorgqr+lorgqropt-1,
516 $ iorglq+lorglqopt-1 )
517 work(1) = lworkopt
518 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
519 info = -19
520 END IF
521 IF( lrwork .LT. lrworkmin .AND. .NOT.lquery ) THEN
522 info = -21
523 END IF
524 END IF
525 IF( info .NE. 0 ) THEN
526 CALL xerbla( 'ZUNCSD2BY1', -info )
527 RETURN
528 ELSE IF( lquery ) THEN
529 RETURN
530 END IF
531 lorgqr = lwork-iorgqr+1
532 lorglq = lwork-iorglq+1
533*
534* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
535* in which R = MIN(P,M-P,Q,M-Q)
536*
537 IF( r .EQ. q ) THEN
538*
539* Case 1: R = Q
540*
541* Simultaneously bidiagonalize X11 and X21
542*
543 CALL zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
544 $ rwork(iphi), work(itaup1), work(itaup2),
545 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
546*
547* Accumulate Householder reflectors
548*
549 IF( wantu1 .AND. p .GT. 0 ) THEN
550 CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
551 CALL zungqr( p, p, q, u1, ldu1, work(itaup1),
552 $ work(iorgqr),
553 $ lorgqr, childinfo )
554 END IF
555 IF( wantu2 .AND. m-p .GT. 0 ) THEN
556 CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
557 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
558 $ work(iorgqr), lorgqr, childinfo )
559 END IF
560 IF( wantv1t .AND. q .GT. 0 ) THEN
561 v1t(1,1) = one
562 DO j = 2, q
563 v1t(1,j) = zero
564 v1t(j,1) = zero
565 END DO
566 CALL zlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
567 $ ldv1t )
568 CALL zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t,
569 $ work(itauq1),
570 $ work(iorglq), lorglq, childinfo )
571 END IF
572*
573* Simultaneously diagonalize X11 and X21.
574*
575 CALL zbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
576 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
577 $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
578 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
579 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
580 $ lrwork-ibbcsd+1, childinfo )
581*
582* Permute rows and columns to place zero submatrices in
583* preferred positions
584*
585 IF( q .GT. 0 .AND. wantu2 ) THEN
586 DO i = 1, q
587 iwork(i) = m - p - q + i
588 END DO
589 DO i = q + 1, m - p
590 iwork(i) = i - q
591 END DO
592 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
593 END IF
594 ELSE IF( r .EQ. p ) THEN
595*
596* Case 2: R = P
597*
598* Simultaneously bidiagonalize X11 and X21
599*
600 CALL zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
601 $ rwork(iphi), work(itaup1), work(itaup2),
602 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
603*
604* Accumulate Householder reflectors
605*
606 IF( wantu1 .AND. p .GT. 0 ) THEN
607 u1(1,1) = one
608 DO j = 2, p
609 u1(1,j) = zero
610 u1(j,1) = zero
611 END DO
612 CALL zlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2),
613 $ ldu1 )
614 CALL zungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
615 $ work(iorgqr), lorgqr, childinfo )
616 END IF
617 IF( wantu2 .AND. m-p .GT. 0 ) THEN
618 CALL zlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
619 CALL zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
620 $ work(iorgqr), lorgqr, childinfo )
621 END IF
622 IF( wantv1t .AND. q .GT. 0 ) THEN
623 CALL zlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
624 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
625 $ work(iorglq), lorglq, childinfo )
626 END IF
627*
628* Simultaneously diagonalize X11 and X21.
629*
630 CALL zbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
631 $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
632 $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
633 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
634 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
635 $ childinfo )
636*
637* Permute rows and columns to place identity submatrices in
638* preferred positions
639*
640 IF( q .GT. 0 .AND. wantu2 ) THEN
641 DO i = 1, q
642 iwork(i) = m - p - q + i
643 END DO
644 DO i = q + 1, m - p
645 iwork(i) = i - q
646 END DO
647 CALL zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
648 END IF
649 ELSE IF( r .EQ. m-p ) THEN
650*
651* Case 3: R = M-P
652*
653* Simultaneously bidiagonalize X11 and X21
654*
655 CALL zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
656 $ rwork(iphi), work(itaup1), work(itaup2),
657 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
658*
659* Accumulate Householder reflectors
660*
661 IF( wantu1 .AND. p .GT. 0 ) THEN
662 CALL zlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
663 CALL zungqr( p, p, q, u1, ldu1, work(itaup1),
664 $ work(iorgqr),
665 $ lorgqr, childinfo )
666 END IF
667 IF( wantu2 .AND. m-p .GT. 0 ) THEN
668 u2(1,1) = one
669 DO j = 2, m-p
670 u2(1,j) = zero
671 u2(j,1) = zero
672 END DO
673 CALL zlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
674 $ ldu2 )
675 CALL zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
676 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
677 END IF
678 IF( wantv1t .AND. q .GT. 0 ) THEN
679 CALL zlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
680 CALL zunglq( q, q, r, v1t, ldv1t, work(itauq1),
681 $ work(iorglq), lorglq, childinfo )
682 END IF
683*
684* Simultaneously diagonalize X11 and X21.
685*
686 CALL zbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
687 $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
688 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
689 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
690 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
691 $ rwork(ibbcsd), lbbcsd, childinfo )
692*
693* Permute rows and columns to place identity submatrices in
694* preferred positions
695*
696 IF( q .GT. r ) THEN
697 DO i = 1, r
698 iwork(i) = q - r + i
699 END DO
700 DO i = r + 1, q
701 iwork(i) = i - r
702 END DO
703 IF( wantu1 ) THEN
704 CALL zlapmt( .false., p, q, u1, ldu1, iwork )
705 END IF
706 IF( wantv1t ) THEN
707 CALL zlapmr( .false., q, q, v1t, ldv1t, iwork )
708 END IF
709 END IF
710 ELSE
711*
712* Case 4: R = M-Q
713*
714* Simultaneously bidiagonalize X11 and X21
715*
716 CALL zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
717 $ rwork(iphi), work(itaup1), work(itaup2),
718 $ work(itauq1), work(iorbdb), work(iorbdb+m),
719 $ lorbdb-m, childinfo )
720*
721* Accumulate Householder reflectors
722*
723 IF( wantu2 .AND. m-p .GT. 0 ) THEN
724 CALL zcopy( m-p, work(iorbdb+p), 1, u2, 1 )
725 END IF
726 IF( wantu1 .AND. p .GT. 0 ) THEN
727 CALL zcopy( p, work(iorbdb), 1, u1, 1 )
728 DO j = 2, p
729 u1(1,j) = zero
730 END DO
731 CALL zlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
732 $ ldu1 )
733 CALL zungqr( p, p, m-q, u1, ldu1, work(itaup1),
734 $ work(iorgqr), lorgqr, childinfo )
735 END IF
736 IF( wantu2 .AND. m-p .GT. 0 ) THEN
737 DO j = 2, m-p
738 u2(1,j) = zero
739 END DO
740 CALL zlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
741 $ ldu2 )
742 CALL zungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
743 $ work(iorgqr), lorgqr, childinfo )
744 END IF
745 IF( wantv1t .AND. q .GT. 0 ) THEN
746 CALL zlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
747 CALL zlacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1),
748 $ ldx11,
749 $ v1t(m-q+1,m-q+1), ldv1t )
750 CALL zlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
751 $ v1t(p+1,p+1), ldv1t )
752 CALL zunglq( q, q, q, v1t, ldv1t, work(itauq1),
753 $ work(iorglq), lorglq, childinfo )
754 END IF
755*
756* Simultaneously diagonalize X11 and X21.
757*
758 CALL zbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
759 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
760 $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
761 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
762 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
763 $ rwork(ibbcsd), lbbcsd, childinfo )
764*
765* Permute rows and columns to place identity submatrices in
766* preferred positions
767*
768 IF( p .GT. r ) THEN
769 DO i = 1, r
770 iwork(i) = p - r + i
771 END DO
772 DO i = r + 1, p
773 iwork(i) = i - r
774 END DO
775 IF( wantu1 ) THEN
776 CALL zlapmt( .false., p, p, u1, ldu1, iwork )
777 END IF
778 IF( wantv1t ) THEN
779 CALL zlapmr( .false., p, q, v1t, ldv1t, iwork )
780 END IF
781 END IF
782 END IF
783*
784 RETURN
785*
786* End of ZUNCSD2BY1
787*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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:331
zcopy
subroutine zcopy(n, zx, incx, zy, incy)
ZCOPY
Definition zcopy.f:81
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:101
zlapmr
subroutine zlapmr(forwrd, m, n, x, ldx, k)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition zlapmr.f:102
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:102
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
zunbdb1
subroutine zunbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB1
Definition zunbdb1.f:202
zunbdb2
subroutine zunbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB2
Definition zunbdb2.f:200
zunbdb3
subroutine zunbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
ZUNBDB3
Definition zunbdb3.f:200
zunbdb4
subroutine zunbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
ZUNBDB4
Definition zunbdb4.f:212
zunglq
subroutine zunglq(m, n, k, a, lda, tau, work, lwork, info)
ZUNGLQ
Definition zunglq.f:125
zungqr
subroutine zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
Definition zungqr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: