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 cos, int, max, min, sin
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, 0, 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, 0, 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 claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
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
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 clapmr(FORWRD, M, N, X, LDX, K)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
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 clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
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.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine clapmt(FORWRD, M, N, X, LDX, K)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine cunglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGLQ
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR