316 RECURSIVE SUBROUTINE cuncsd( 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
337 COMPLEX U1( ldu1, * ), U2( ldu2, * ), V1T( ldv1t, * ),
338 $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
339 $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
347 parameter ( one = (1.0e0,0.0e0),
348 $ zero = (0.0e0,0.0e0) )
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 cuncsd( 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 cuncsd( 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 cbbcsd( 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 cungqr( 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 cunglq( 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 cunbdb( 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(
'CUNCSD', -info )
531 ELSE IF( lquery .OR. lrquery )
THEN
537 CALL cunbdb( 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 clacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
547 CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
550 IF( wantu2 .AND. m-p .GT. 0 )
THEN
551 CALL clacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
552 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
553 $ work(iorgqr), lorgqrwork, info )
555 IF( wantv1t .AND. q .GT. 0 )
THEN
556 CALL clacpy(
'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
563 CALL cunglq( 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 clacpy(
'U', p, m-q, x12, ldx12, v2t, ldv2t )
568 IF( m-p .GT. q )
THEN
569 CALL clacpy(
'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
570 $ v2t(p+1,p+1), ldv2t )
573 CALL cunglq( 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 clacpy(
'U', q, p, x11, ldx11, u1, ldu1 )
580 CALL cunglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
583 IF( wantu2 .AND. m-p .GT. 0 )
THEN
584 CALL clacpy(
'U', q, m-p, x21, ldx21, u2, ldu2 )
585 CALL cunglq( m-p, m-p, q, u2, ldu2, work(itaup2),
586 $ work(iorglq), lorglqwork, info )
588 IF( wantv1t .AND. q .GT. 0 )
THEN
589 CALL clacpy(
'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
596 CALL cungqr( 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 clacpy(
'L', m-q, p, x12, ldx12, v2t, ldv2t )
603 IF ( m .GT. p+q )
THEN
604 CALL clacpy(
'L', m-p-q, m-p-q, x22(p1,q1), ldx22,
605 $ v2t(p+1,p+1), ldv2t )
607 CALL cungqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
608 $ work(iorgqr), lorgqrwork, info )
614 CALL cbbcsd( 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 clapmt( .false., m-p, m-p, u2, ldu2, iwork )
636 CALL clapmr( .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 clapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
649 CALL clapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
subroutine cbbcsd(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)
CBBCSD
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine cunbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
CUNBDB
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
subroutine clapmr(FORWRD, M, N, X, LDX, K)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
recursive subroutine cuncsd(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)
CUNCSD