297 RECURSIVE SUBROUTINE sorcsd( JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS,
298 $ signs, m, p, q, x11, ldx11, x12,
299 $ ldx12, x21, ldx21, x22, ldx22, theta,
300 $ u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
301 $ ldv2t, work, lwork, iwork, info )
309 CHARACTER JOBU1, JOBU2, JOBV1T, JOBV2T, SIGNS, TRANS
310 INTEGER INFO, LDU1, LDU2, LDV1T, LDV2T, LDX11, LDX12,
311 $ ldx21, ldx22, lwork, m, p, q
316 REAL U1( ldu1, * ), U2( ldu2, * ), V1T( ldv1t, * ),
317 $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
318 $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
326 parameter ( one = 1.0e+0,
333 CHARACTER TRANST, SIGNST
334 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
335 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
336 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
337 $ itauq2, j, lbbcsdwork, lbbcsdworkmin,
338 $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
339 $ lorbdbworkopt, lorglqwork, lorglqworkmin,
340 $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
341 $ lorgqrworkopt, lworkmin, lworkopt
342 LOGICAL COLMAJOR, DEFAULTSIGNS, LQUERY, WANTU1, WANTU2,
354 INTRINSIC int, max, min
361 wantu1 = lsame( jobu1,
'Y' )
362 wantu2 = lsame( jobu2,
'Y' )
363 wantv1t = lsame( jobv1t,
'Y' )
364 wantv2t = lsame( jobv2t,
'Y' )
365 colmajor = .NOT. lsame( trans,
'T' )
366 defaultsigns = .NOT. lsame( signs,
'O' )
367 lquery = lwork .EQ. -1
370 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
372 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
374 ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) )
THEN
376 ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) )
THEN
378 ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) )
THEN
380 ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) )
THEN
382 ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) )
THEN
384 ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) )
THEN
386 ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) )
THEN
388 ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) )
THEN
390 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
392 ELSE IF( wantu2 .AND. ldu2 .LT. m-p )
THEN
394 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
396 ELSE IF( wantv2t .AND. ldv2t .LT. m-q )
THEN
402 IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) )
THEN
408 IF( defaultsigns )
THEN
413 CALL sorcsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
414 $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
415 $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
416 $ u2, ldu2, work, lwork, iwork, info )
423 IF( info .EQ. 0 .AND. m-q .LT. q )
THEN
424 IF( defaultsigns )
THEN
429 CALL sorcsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
430 $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
431 $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
432 $ ldv1t, work, lwork, iwork, info )
438 IF( info .EQ. 0 )
THEN
441 itaup1 = iphi + max( 1, q - 1 )
442 itaup2 = itaup1 + max( 1, p )
443 itauq1 = itaup2 + max( 1, m - p )
444 itauq2 = itauq1 + max( 1, q )
445 iorgqr = itauq2 + max( 1, m - q )
446 CALL sorgqr( m-q, m-q, m-q, dummy, max(1,m-q), dummy, work, -1,
448 lorgqrworkopt = int( work(1) )
449 lorgqrworkmin = max( 1, m - q )
450 iorglq = itauq2 + max( 1, m - q )
451 CALL sorglq( m-q, m-q, m-q, dummy, max(1,m-q), dummy, work, -1,
453 lorglqworkopt = int( work(1) )
454 lorglqworkmin = max( 1, m - q )
455 iorbdb = itauq2 + max( 1, m - q )
456 CALL sorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
457 $ x21, ldx21, x22, ldx22, dummy, dummy, dummy, dummy, dummy,
458 $ dummy,work,-1,childinfo )
459 lorbdbworkopt = int( work(1) )
460 lorbdbworkmin = lorbdbworkopt
461 ib11d = itauq2 + max( 1, m - q )
462 ib11e = ib11d + max( 1, q )
463 ib12d = ib11e + max( 1, q - 1 )
464 ib12e = ib12d + max( 1, q )
465 ib21d = ib12e + max( 1, q - 1 )
466 ib21e = ib21d + max( 1, q )
467 ib22d = ib21e + max( 1, q - 1 )
468 ib22e = ib22d + max( 1, q )
469 ibbcsd = ib22e + max( 1, q - 1 )
470 CALL sbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
471 $ dummy, dummy, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
472 $ ldv2t, dummy, dummy, dummy, dummy, dummy, dummy,
473 $ dummy, dummy, work, -1, childinfo )
474 lbbcsdworkopt = int( work(1) )
475 lbbcsdworkmin = lbbcsdworkopt
476 lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
477 $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkopt ) - 1
478 lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
479 $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkmin ) - 1
480 work(1) = max(lworkopt,lworkmin)
482 IF( lwork .LT. lworkmin .AND. .NOT. lquery )
THEN
485 lorgqrwork = lwork - iorgqr + 1
486 lorglqwork = lwork - iorglq + 1
487 lorbdbwork = lwork - iorbdb + 1
488 lbbcsdwork = lwork - ibbcsd + 1
494 IF( info .NE. 0 )
THEN
495 CALL xerbla(
'SORCSD', -info )
497 ELSE IF( lquery )
THEN
503 CALL sorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
504 $ ldx21, x22, ldx22, theta, work(iphi), work(itaup1),
505 $ work(itaup2), work(itauq1), work(itauq2),
506 $ work(iorbdb), lorbdbwork, childinfo )
511 IF( wantu1 .AND. p .GT. 0 )
THEN
512 CALL slacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
513 CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
516 IF( wantu2 .AND. m-p .GT. 0 )
THEN
517 CALL slacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
518 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
519 $ work(iorgqr), lorgqrwork, info )
521 IF( wantv1t .AND. q .GT. 0 )
THEN
522 CALL slacpy(
'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
529 CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
530 $ work(iorglq), lorglqwork, info )
532 IF( wantv2t .AND. m-q .GT. 0 )
THEN
533 CALL slacpy(
'U', p, m-q, x12, ldx12, v2t, ldv2t )
534 CALL slacpy(
'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
535 $ v2t(p+1,p+1), ldv2t )
536 CALL sorglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
537 $ work(iorglq), lorglqwork, info )
540 IF( wantu1 .AND. p .GT. 0 )
THEN
541 CALL slacpy(
'U', q, p, x11, ldx11, u1, ldu1 )
542 CALL sorglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
545 IF( wantu2 .AND. m-p .GT. 0 )
THEN
546 CALL slacpy(
'U', q, m-p, x21, ldx21, u2, ldu2 )
547 CALL sorglq( m-p, m-p, q, u2, ldu2, work(itaup2),
548 $ work(iorglq), lorglqwork, info )
550 IF( wantv1t .AND. q .GT. 0 )
THEN
551 CALL slacpy(
'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
558 CALL sorgqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
559 $ work(iorgqr), lorgqrwork, info )
561 IF( wantv2t .AND. m-q .GT. 0 )
THEN
562 CALL slacpy(
'L', m-q, p, x12, ldx12, v2t, ldv2t )
563 CALL slacpy(
'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
564 $ v2t(p+1,p+1), ldv2t )
565 CALL sorgqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
566 $ work(iorgqr), lorgqrwork, info )
572 CALL sbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
573 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
574 $ ldv2t, work(ib11d), work(ib11e), work(ib12d),
575 $ work(ib12e), work(ib21d), work(ib21e), work(ib22d),
576 $ work(ib22e), work(ibbcsd), lbbcsdwork, info )
583 IF( q .GT. 0 .AND. wantu2 )
THEN
585 iwork(i) = m - p - q + i
591 CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
593 CALL slapmr( .false., m-p, m-p, u2, ldu2, iwork )
596 IF( m .GT. 0 .AND. wantv2t )
THEN
598 iwork(i) = m - p - q + i
603 IF( .NOT. colmajor )
THEN
604 CALL slapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
606 CALL slapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
subroutine sbbcsd(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, WORK, LWORK, INFO)
SBBCSD
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ
subroutine slapmr(FORWRD, M, N, X, LDX, K)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
recursive subroutine sorcsd(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, IWORK, INFO)
SORCSD
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sorbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
SORBDB