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 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 zunbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
ZUNBDB
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlapmt(FORWRD, M, N, X, LDX, K)
ZLAPMT performs a forward or backward permutation of the columns of a matrix.
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 zlapmr(FORWRD, M, N, X, LDX, K)
ZLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine zunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGLQ