297 RECURSIVE SUBROUTINE dorcsd( 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
315 DOUBLE PRECISION theta( * )
316 DOUBLE PRECISION u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
317 $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
318 $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
325 DOUBLE PRECISION one, zero
326 parameter( one = 1.0d0,
330 CHARACTER transt, signst
331 INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
332 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
333 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
334 $ itauq2,
j, lbbcsdwork, lbbcsdworkmin,
335 $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
336 $ lorbdbworkopt, lorglqwork, lorglqworkmin,
337 $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
338 $ lorgqrworkopt, lworkmin, lworkopt
339 LOGICAL colmajor, defaultsigns, lquery, wantu1, wantu2,
351 INTRINSIC int, max, min
358 wantu1 =
lsame( jobu1,
'Y' )
359 wantu2 =
lsame( jobu2,
'Y' )
360 wantv1t =
lsame( jobv1t,
'Y' )
361 wantv2t =
lsame( jobv2t,
'Y' )
362 colmajor = .NOT.
lsame( trans,
'T' )
363 defaultsigns = .NOT.
lsame( signs,
'O' )
364 lquery = lwork .EQ. -1
367 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
369 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
371 ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) )
THEN
373 ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) )
THEN
375 ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) )
THEN
377 ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) )
THEN
379 ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) )
THEN
381 ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) )
THEN
383 ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) )
THEN
385 ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) )
THEN
387 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
389 ELSE IF( wantu2 .AND. ldu2 .LT. m-p )
THEN
391 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
393 ELSE IF( wantv2t .AND. ldv2t .LT. m-q )
THEN
399 IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) )
THEN
405 IF( defaultsigns )
THEN
410 CALL
dorcsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst, m,
411 $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
412 $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
413 $ u2, ldu2, work, lwork, iwork, info )
420 IF( info .EQ. 0 .AND. m-q .LT. q )
THEN
421 IF( defaultsigns )
THEN
426 CALL
dorcsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
427 $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
428 $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
429 $ ldv1t, work, lwork, iwork, info )
435 IF( info .EQ. 0 )
THEN
438 itaup1 = iphi + max( 1, q - 1 )
439 itaup2 = itaup1 + max( 1, p )
440 itauq1 = itaup2 + max( 1, m - p )
441 itauq2 = itauq1 + max( 1, q )
442 iorgqr = itauq2 + max( 1, m - q )
443 CALL
dorgqr( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
445 lorgqrworkopt = int( work(1) )
446 lorgqrworkmin = max( 1, m - q )
447 iorglq = itauq2 + max( 1, m - q )
448 CALL
dorglq( m-q, m-q, m-q, u1, max(1,m-q), u1, work, -1,
450 lorglqworkopt = int( work(1) )
451 lorglqworkmin = max( 1, m - q )
452 iorbdb = itauq2 + max( 1, m - q )
453 CALL
dorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
454 $ x21, ldx21, x22, ldx22, theta, v1t, u1, u2, v1t,
455 $ v2t, work, -1, childinfo )
456 lorbdbworkopt = int( work(1) )
457 lorbdbworkmin = lorbdbworkopt
458 ib11d = itauq2 + max( 1, m - q )
459 ib11e = ib11d + max( 1, q )
460 ib12d = ib11e + max( 1, q - 1 )
461 ib12e = ib12d + max( 1, q )
462 ib21d = ib12e + max( 1, q - 1 )
463 ib21e = ib21d + max( 1, q )
464 ib22d = ib21e + max( 1, q - 1 )
465 ib22e = ib22d + max( 1, q )
466 ibbcsd = ib22e + max( 1, q - 1 )
467 CALL
dbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
468 $ theta, theta, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
469 $ ldv2t, u1, u1, u1, u1, u1, u1, u1, u1, work, -1,
471 lbbcsdworkopt = int( work(1) )
472 lbbcsdworkmin = lbbcsdworkopt
473 lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
474 $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkopt ) - 1
475 lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
476 $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkmin ) - 1
477 work(1) = max(lworkopt,lworkmin)
479 IF( lwork .LT. lworkmin .AND. .NOT. lquery )
THEN
482 lorgqrwork = lwork - iorgqr + 1
483 lorglqwork = lwork - iorglq + 1
484 lorbdbwork = lwork - iorbdb + 1
485 lbbcsdwork = lwork - ibbcsd + 1
491 IF( info .NE. 0 )
THEN
492 CALL
xerbla(
'DORCSD', -info )
494 ELSE IF( lquery )
THEN
500 CALL
dorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12, x21,
501 $ ldx21, x22, ldx22, theta, work(iphi), work(itaup1),
502 $ work(itaup2), work(itauq1), work(itauq2),
503 $ work(iorbdb), lorbdbwork, childinfo )
508 IF( wantu1 .AND. p .GT. 0 )
THEN
509 CALL
dlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
510 CALL
dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
513 IF( wantu2 .AND. m-p .GT. 0 )
THEN
514 CALL
dlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
515 CALL
dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
516 $ work(iorgqr), lorgqrwork, info )
518 IF( wantv1t .AND. q .GT. 0 )
THEN
519 CALL
dlacpy(
'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
526 CALL
dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
527 $ work(iorglq), lorglqwork, info )
529 IF( wantv2t .AND. m-q .GT. 0 )
THEN
530 CALL
dlacpy(
'U', p, m-q, x12, ldx12, v2t, ldv2t )
532 CALL
dlacpy(
'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
533 $ v2t(p+1,p+1), ldv2t )
536 CALL
dorglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
537 $ work(iorglq), lorglqwork, info )
541 IF( wantu1 .AND. p .GT. 0 )
THEN
542 CALL
dlacpy(
'U', q, p, x11, ldx11, u1, ldu1 )
543 CALL
dorglq( p, p, q, u1, ldu1, work(itaup1), work(iorglq),
546 IF( wantu2 .AND. m-p .GT. 0 )
THEN
547 CALL
dlacpy(
'U', q, m-p, x21, ldx21, u2, ldu2 )
548 CALL
dorglq( m-p, m-p, q, u2, ldu2, work(itaup2),
549 $ work(iorglq), lorglqwork, info )
551 IF( wantv1t .AND. q .GT. 0 )
THEN
552 CALL
dlacpy(
'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
559 CALL
dorgqr( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
560 $ work(iorgqr), lorgqrwork, info )
562 IF( wantv2t .AND. m-q .GT. 0 )
THEN
563 CALL
dlacpy(
'L', m-q, p, x12, ldx12, v2t, ldv2t )
564 CALL
dlacpy(
'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
565 $ v2t(p+1,p+1), ldv2t )
566 CALL
dorgqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
567 $ work(iorgqr), lorgqrwork, info )
573 CALL
dbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q, theta,
574 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
575 $ ldv2t, work(ib11d), work(ib11e), work(ib12d),
576 $ work(ib12e), work(ib21d), work(ib21e), work(ib22d),
577 $ work(ib22e), work(ibbcsd), lbbcsdwork, info )
584 IF( q .GT. 0 .AND. wantu2 )
THEN
586 iwork(i) = m - p - q + i
592 CALL
dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
594 CALL
dlapmr( .false., m-p, m-p, u2, ldu2, iwork )
597 IF( m .GT. 0 .AND. wantv2t )
THEN
599 iwork(i) = m - p - q + i
604 IF( .NOT. colmajor )
THEN
605 CALL
dlapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
607 CALL
dlapmr( .false., m-q, m-q, v2t, ldv2t, iwork )
LOGICAL function lsame(CA, CB)
LSAME
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
subroutine dorbdb(TRANS, SIGNS, M, P, Q, X11, LDX11, X12, LDX12, X21, LDX21, X22, LDX22, THETA, PHI, TAUP1, TAUP2, TAUQ1, TAUQ2, WORK, LWORK, INFO)
DORBDB
recursive subroutine dorcsd(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)
DORCSD
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dbbcsd(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)
DBBCSD
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dlapmt(FORWRD, M, N, X, LDX, K)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine dlapmr(FORWRD, M, N, X, LDX, K)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.