316 RECURSIVE SUBROUTINE zuncsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
317 $ signs, m, p, q, x11, ldx11, x12,
318 $ ldx12, x21, ldx21, x22, ldx22, theta,
319 $ u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
320 $ ldv2t, work, lwork, rwork, lrwork,
329 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
330 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
331 $ ldx21, ldx22, lrwork, lwork, m, p, q
335 DOUBLE PRECISION THETA( * )
336 DOUBLE PRECISION RWORK( * )
337 COMPLEX*16 U1( ldu1, * ), U2( ldu2, * ), V1T( ldv1t, * ),
338 $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
339 $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
347 parameter( one = (1.0d0,0.0d0),
348 $ zero = (0.0d0,0.0d0) )
351 CHARACTER TRANST, SIGNST
352 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
353 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
354 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
355 $ itauq2, j, lbbcsdwork, lbbcsdworkmin,
356 $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
357 $ lorbdbworkopt, lorglqwork, lorglqworkmin,
358 $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
359 $ lorgqrworkopt, lworkmin, lworkopt, p1, q1
360 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
362 INTEGER LRWORKMIN, LRWORKOPT
374 INTRINSIC int, max, min
381 wantu1 = lsame( jobu1,
'Y' )
382 wantu2 = lsame( jobu2,
'Y' )
383 wantv1t = lsame( jobv1t,
'Y' )
384 wantv2t = lsame( jobv2t,
'Y' )
385 colmajor = .NOT. lsame( trans,
'T' )
386 defaultsigns = .NOT. lsame( signs,
'O' )
387 lquery = lwork .EQ. -1
388 lrquery = lrwork .EQ. -1
391 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
393 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
395 ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) )
THEN
397 ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) )
THEN
399 ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) )
THEN
401 ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) )
THEN
403 ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) )
THEN
405 ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) )
THEN
407 ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) )
THEN
409 ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) )
THEN
411 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
413 ELSE IF( wantu2 .AND. ldu2 .LT. m-p )
THEN
415 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
417 ELSE IF( wantv2t .AND. ldv2t .LT. m-q )
THEN
423 IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) )
THEN
429 IF( defaultsigns )
THEN
434 CALL
zuncsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
435 $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
436 $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
437 $ u2, ldu2, work, lwork, rwork, lrwork, iwork,
445 IF( info .EQ. 0 .AND. m-q .LT. q )
THEN
446 IF( defaultsigns )
THEN
451 CALL
zuncsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
452 $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
453 $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
454 $ ldv1t, work, lwork, rwork, lrwork, iwork, info )
460 IF( info .EQ. 0 )
THEN
465 ib11d = iphi + max( 1, q - 1 )
466 ib11e = ib11d + max( 1, q )
467 ib12d = ib11e + max( 1, q - 1 )
468 ib12e = ib12d + max( 1, q )
469 ib21d = ib12e + max( 1, q - 1 )
470 ib21e = ib21d + max( 1, q )
471 ib22d = ib21e + max( 1, q - 1 )
472 ib22e = ib22d + max( 1, q )
473 ibbcsd = ib22e + max( 1, q - 1 )
474 CALL
zbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
475 $ theta, theta, u1, ldu1, u2, ldu2, v1t, ldv1t,
476 $ v2t, ldv2t, theta, theta, theta, theta, theta,
477 $ theta, theta, theta, rwork, -1, childinfo )
478 lbbcsdworkopt = int( rwork(1) )
479 lbbcsdworkmin = lbbcsdworkopt
480 lrworkopt = ibbcsd + lbbcsdworkopt - 1
481 lrworkmin = ibbcsd + lbbcsdworkmin - 1
487 itaup2 = itaup1 + max( 1, p )
488 itauq1 = itaup2 + max( 1, m - p )
489 itauq2 = itauq1 + max( 1, q )
490 iorgqr = itauq2 + max( 1, m - q )
491 CALL
zungqr( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
493 lorgqrworkopt = int( work(1) )
494 lorgqrworkmin = max( 1, m - q )
495 iorglq = itauq2 + max( 1, m - q )
496 CALL
zunglq( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
498 lorglqworkopt = int( work(1) )
499 lorglqworkmin = max( 1, m - q )
500 iorbdb = itauq2 + max( 1, m - q )
501 CALL
zunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
502 $ x21, ldx21, x22, ldx22, theta, theta, u1, u2,
503 $ v1t, v2t, work, -1, childinfo )
504 lorbdbworkopt = int( work(1) )
505 lorbdbworkmin = lorbdbworkopt
506 lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
507 $ iorbdb + lorbdbworkopt ) - 1
508 lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
509 $ iorbdb + lorbdbworkmin ) - 1
510 work(1) = max(lworkopt,lworkmin)
512 IF( lwork .LT. lworkmin
513 $ .AND. .NOT. ( lquery .OR. lrquery ) )
THEN
515 ELSE IF( lrwork .LT. lrworkmin
516 $ .AND. .NOT. ( lquery .OR. lrquery ) )
THEN
519 lorgqrwork = lwork - iorgqr + 1
520 lorglqwork = lwork - iorglq + 1
521 lorbdbwork = lwork - iorbdb + 1
522 lbbcsdwork = lrwork - ibbcsd + 1
528 IF( info .NE. 0 )
THEN
529 CALL
xerbla(
'ZUNCSD', -info )
531 ELSE IF( lquery .OR. lrquery )
THEN
537 CALL
zunbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
538 $ ldx21, x22, ldx22, theta, rwork(iphi), work(itaup1),
539 $ work(itaup2), work(itauq1), work(itauq2),
540 $ work(iorbdb), lorbdbwork, childinfo )
545 IF( wantu1 .AND. p .GT. 0 )
THEN
546 CALL
zlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
547 CALL
zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
550 IF( wantu2 .AND. m-p .GT. 0 )
THEN
551 CALL
zlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
552 CALL
zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
553 $ work(iorgqr), lorgqrwork, info )
555 IF( wantv1t .AND. q .GT. 0 )
THEN
556 CALL
zlacpy(
'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
563 CALL
zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
564 $ work(iorglq), lorglqwork, info )
566 IF( wantv2t .AND. m-q .GT. 0 )
THEN
567 CALL
zlacpy(
'U', p, m-q, x12, ldx12, v2t, ldv2t )
569 CALL
zlacpy(
'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
570 $ v2t(p+1,p+1), ldv2t )
573 CALL
zunglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
574 $ work(iorglq), lorglqwork, info )
578 IF( wantu1 .AND. p .GT. 0 )
THEN
579 CALL
zlacpy(
'U', q, p, x11, ldx11, u1, ldu1 )
580 CALL
zunglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
583 IF( wantu2 .AND. m-p .GT. 0 )
THEN
584 CALL
zlacpy(
'U', q, m-p, x21, ldx21, u2, ldu2 )
585 CALL
zunglq( m-p, m-p, q, u2, ldu2, work(itaup2),
586 $ work(iorglq), lorglqwork, info )
588 IF( wantv1t .AND. q .GT. 0 )
THEN
589 CALL
zlacpy(
'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
596 CALL
zungqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
597 $ work(iorgqr), lorgqrwork, info )
599 IF( wantv2t .AND. m-q .GT. 0 )
THEN
602 CALL
zlacpy(
'L', m-q, p, x12, ldx12, v2t, ldv2t )
603 IF( m .GT. p+q )
THEN
604 CALL
zlacpy(
'L', m-p-q, m-p-q, x22(p1,q1), ldx22,
605 $ v2t(p+1,p+1), ldv2t )
607 CALL
zungqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
608 $ work(iorgqr), lorgqrwork, info )
614 CALL
zbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
615 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
616 $ ldv2t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
617 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
618 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
626 IF( q .GT. 0 .AND. wantu2 )
THEN
628 iwork(i) = m - p - q + i
634 CALL
zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
636 CALL
zlapmr( .false., m-p, m-p, u2, ldu2, iwork )
639 IF( m .GT. 0 .AND. wantv2t )
THEN
641 iwork(i) = m - p - q + i
646 IF( .NOT. colmajor )
THEN
647 CALL
zlapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
649 CALL
zlapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
subroutine zlapmt(FORWRD, M, N, X, LDX, K)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zunbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
ZUNBDB
recursive subroutine zuncsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, WORK, LWORK, RWORK, LRWORK, IWORK, INFO)
ZUNCSD
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ
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
subroutine zlapmr(FORWRD, M, N, X, LDX, K)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.