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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine slapmr(FORWRD, M, N, X, LDX, K)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine xerbla(SRNAME, INFO)
XERBLA
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
logical function lsame(CA, CB)
LSAME
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ
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
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