242 CHARACTER jobu1, jobu2, jobv1t
243 INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
247 DOUBLE PRECISION theta(*)
248 DOUBLE PRECISION u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
249 $ x11(ldx11,*), x21(ldx21,*)
256 DOUBLE PRECISION one, zero
257 parameter ( one = 1.0d0, zero = 0.0d0 )
260 INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
261 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
262 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
263 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
264 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
265 $ lworkmin, lworkopt, r
266 LOGICAL lquery, wantu1, wantu2, wantv1t
269 DOUBLE PRECISION dum1(1), dum2(1,1)
281 INTRINSIC int, max, min
288 wantu1 =
lsame( jobu1,
'Y' )
289 wantu2 =
lsame( jobu2,
'Y' )
290 wantv1t =
lsame( jobv1t,
'Y' )
291 lquery = lwork .EQ. -1
295 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
297 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
299 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
301 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
303 ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) )
THEN
305 ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) )
THEN
307 ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) )
THEN
311 r = min( p, m-p, q, m-q )
332 IF( info .EQ. 0 )
THEN
334 ib11d = iphi + max( 1, r-1 )
335 ib11e = ib11d + max( 1, r )
336 ib12d = ib11e + max( 1, r - 1 )
337 ib12e = ib12d + max( 1, r )
338 ib21d = ib12e + max( 1, r - 1 )
339 ib21e = ib21d + max( 1, r )
340 ib22d = ib21e + max( 1, r - 1 )
341 ib22e = ib22d + max( 1, r )
342 ibbcsd = ib22e + max( 1, r - 1 )
343 itaup1 = iphi + max( 1, r-1 )
344 itaup2 = itaup1 + max( 1, p )
345 itauq1 = itaup2 + max( 1, m-p )
346 iorbdb = itauq1 + max( 1, q )
347 iorgqr = itauq1 + max( 1, q )
348 iorglq = itauq1 + max( 1, q )
354 CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
355 $ dum1, dum1, dum1, dum1, work,
357 lorbdb = int( work(1) )
358 IF( wantu1 .AND. p .GT. 0 )
THEN
359 CALL dorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
361 lorgqrmin = max( lorgqrmin, p )
362 lorgqropt = max( lorgqropt, int( work(1) ) )
364 IF( wantu2 .AND. m-p .GT. 0 )
THEN
365 CALL dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
367 lorgqrmin = max( lorgqrmin, m-p )
368 lorgqropt = max( lorgqropt, int( work(1) ) )
370 IF( wantv1t .AND. q .GT. 0 )
THEN
371 CALL dorglq( q-1, q-1, q-1, v1t, ldv1t,
372 $ dum1, work(1), -1, childinfo )
373 lorglqmin = max( lorglqmin, q-1 )
374 lorglqopt = max( lorglqopt, int( work(1) ) )
376 CALL dbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
377 $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t,
378 $ dum2, 1, dum1, dum1, dum1,
379 $ dum1, dum1, dum1, dum1,
380 $ dum1, work(1), -1, childinfo )
381 lbbcsd = int( work(1) )
382 ELSE IF( r .EQ. p )
THEN
383 CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
384 $ dum1, dum1, dum1, dum1,
385 $ work(1), -1, childinfo )
386 lorbdb = int( work(1) )
387 IF( wantu1 .AND. p .GT. 0 )
THEN
388 CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, dum1,
389 $ work(1), -1, childinfo )
390 lorgqrmin = max( lorgqrmin, p-1 )
391 lorgqropt = max( lorgqropt, int( work(1) ) )
393 IF( wantu2 .AND. m-p .GT. 0 )
THEN
394 CALL dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
396 lorgqrmin = max( lorgqrmin, m-p )
397 lorgqropt = max( lorgqropt, int( work(1) ) )
399 IF( wantv1t .AND. q .GT. 0 )
THEN
400 CALL dorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
402 lorglqmin = max( lorglqmin, q )
403 lorglqopt = max( lorglqopt, int( work(1) ) )
405 CALL dbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
406 $ dum1, v1t, ldv1t, dum2, 1, u1, ldu1,
407 $ u2, ldu2, dum1, dum1, dum1,
408 $ dum1, dum1, dum1, dum1,
409 $ dum1, work(1), -1, childinfo )
410 lbbcsd = int( work(1) )
411 ELSE IF( r .EQ. m-p )
THEN
412 CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
413 $ dum1, dum1, dum1, dum1,
414 $ work(1), -1, childinfo )
415 lorbdb = int( work(1) )
416 IF( wantu1 .AND. p .GT. 0 )
THEN
417 CALL dorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
419 lorgqrmin = max( lorgqrmin, p )
420 lorgqropt = max( lorgqropt, int( work(1) ) )
422 IF( wantu2 .AND. m-p .GT. 0 )
THEN
423 CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
424 $ dum1, work(1), -1, childinfo )
425 lorgqrmin = max( lorgqrmin, m-p-1 )
426 lorgqropt = max( lorgqropt, int( work(1) ) )
428 IF( wantv1t .AND. q .GT. 0 )
THEN
429 CALL dorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
431 lorglqmin = max( lorglqmin, q )
432 lorglqopt = max( lorglqopt, int( work(1) ) )
434 CALL dbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
435 $ theta, dum1, dum2, 1, v1t, ldv1t, u2,
436 $ ldu2, u1, ldu1, dum1, dum1, dum1,
437 $ dum1, dum1, dum1, dum1,
438 $ dum1, work(1), -1, childinfo )
439 lbbcsd = int( work(1) )
441 CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
442 $ dum1, dum1, dum1, dum1,
443 $ dum1, work(1), -1, childinfo )
444 lorbdb = m + int( work(1) )
445 IF( wantu1 .AND. p .GT. 0 )
THEN
446 CALL dorgqr( p, p, m-q, u1, ldu1, dum1, work(1), -1,
448 lorgqrmin = max( lorgqrmin, p )
449 lorgqropt = max( lorgqropt, int( work(1) ) )
451 IF( wantu2 .AND. m-p .GT. 0 )
THEN
452 CALL dorgqr( m-p, m-p, m-q, u2, ldu2, dum1, work(1),
454 lorgqrmin = max( lorgqrmin, m-p )
455 lorgqropt = max( lorgqropt, int( work(1) ) )
457 IF( wantv1t .AND. q .GT. 0 )
THEN
458 CALL dorglq( q, q, q, v1t, ldv1t, dum1, work(1), -1,
460 lorglqmin = max( lorglqmin, q )
461 lorglqopt = max( lorglqopt, int( work(1) ) )
463 CALL dbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
464 $ theta, dum1, u2, ldu2, u1, ldu1, dum2,
465 $ 1, v1t, ldv1t, dum1, dum1, dum1,
466 $ dum1, dum1, dum1, dum1,
467 $ dum1, work(1), -1, childinfo )
468 lbbcsd = int( work(1) )
470 lworkmin = max( iorbdb+lorbdb-1,
471 $ iorgqr+lorgqrmin-1,
472 $ iorglq+lorglqmin-1,
474 lworkopt = max( iorbdb+lorbdb-1,
475 $ iorgqr+lorgqropt-1,
476 $ iorglq+lorglqopt-1,
479 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
483 IF( info .NE. 0 )
THEN
484 CALL xerbla(
'DORCSD2BY1', -info )
486 ELSE IF( lquery )
THEN
489 lorgqr = lwork-iorgqr+1
490 lorglq = lwork-iorglq+1
501 CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
502 $ work(iphi), work(itaup1), work(itaup2),
503 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
507 IF( wantu1 .AND. p .GT. 0 )
THEN
508 CALL dlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
509 CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
510 $ lorgqr, childinfo )
512 IF( wantu2 .AND. m-p .GT. 0 )
THEN
513 CALL dlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
514 CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
515 $ work(iorgqr), lorgqr, childinfo )
517 IF( wantv1t .AND. q .GT. 0 )
THEN
523 CALL dlacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
525 CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
526 $ work(iorglq), lorglq, childinfo )
531 CALL dbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
532 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t,
533 $ dum2, 1, work(ib11d), work(ib11e),
534 $ work(ib12d), work(ib12e), work(ib21d),
535 $ work(ib21e), work(ib22d), work(ib22e),
536 $ work(ibbcsd), lbbcsd, childinfo )
541 IF( q .GT. 0 .AND. wantu2 )
THEN
543 iwork(i) = m - p - q + i
548 CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
550 ELSE IF( r .EQ. p )
THEN
556 CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
557 $ work(iphi), work(itaup1), work(itaup2),
558 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
562 IF( wantu1 .AND. p .GT. 0 )
THEN
568 CALL dlacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
569 CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
570 $ work(iorgqr), lorgqr, childinfo )
572 IF( wantu2 .AND. m-p .GT. 0 )
THEN
573 CALL dlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
574 CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
575 $ work(iorgqr), lorgqr, childinfo )
577 IF( wantv1t .AND. q .GT. 0 )
THEN
578 CALL dlacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
579 CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
580 $ work(iorglq), lorglq, childinfo )
585 CALL dbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
586 $ work(iphi), v1t, ldv1t, dum2, 1, u1, ldu1, u2,
587 $ ldu2, work(ib11d), work(ib11e), work(ib12d),
588 $ work(ib12e), work(ib21d), work(ib21e),
589 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
595 IF( q .GT. 0 .AND. wantu2 )
THEN
597 iwork(i) = m - p - q + i
602 CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
604 ELSE IF( r .EQ. m-p )
THEN
610 CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
611 $ work(iphi), work(itaup1), work(itaup2),
612 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
616 IF( wantu1 .AND. p .GT. 0 )
THEN
617 CALL dlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
618 CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
619 $ lorgqr, childinfo )
621 IF( wantu2 .AND. m-p .GT. 0 )
THEN
627 CALL dlacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
629 CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
630 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
632 IF( wantv1t .AND. q .GT. 0 )
THEN
633 CALL dlacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
634 CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
635 $ work(iorglq), lorglq, childinfo )
640 CALL dbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
641 $ theta, work(iphi), dum2, 1, v1t, ldv1t, u2,
642 $ ldu2, u1, ldu1, work(ib11d), work(ib11e),
643 $ work(ib12d), work(ib12e), work(ib21d),
644 $ work(ib21e), work(ib22d), work(ib22e),
645 $ work(ibbcsd), lbbcsd, childinfo )
658 CALL dlapmt( .false., p, q, u1, ldu1, iwork )
661 CALL dlapmr( .false., q, q, v1t, ldv1t, iwork )
670 CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
671 $ work(iphi), work(itaup1), work(itaup2),
672 $ work(itauq1), work(iorbdb), work(iorbdb+m),
673 $ lorbdb-m, childinfo )
677 IF( wantu1 .AND. p .GT. 0 )
THEN
678 CALL dcopy( p, work(iorbdb), 1, u1, 1 )
682 CALL dlacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
684 CALL dorgqr( p, p, m-q, u1, ldu1, work(itaup1),
685 $ work(iorgqr), lorgqr, childinfo )
687 IF( wantu2 .AND. m-p .GT. 0 )
THEN
688 CALL dcopy( m-p, work(iorbdb+p), 1, u2, 1 )
692 CALL dlacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
694 CALL dorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
695 $ work(iorgqr), lorgqr, childinfo )
697 IF( wantv1t .AND. q .GT. 0 )
THEN
698 CALL dlacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
699 CALL dlacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
700 $ v1t(m-q+1,m-q+1), ldv1t )
701 CALL dlacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
702 $ v1t(p+1,p+1), ldv1t )
703 CALL dorglq( q, q, q, v1t, ldv1t, work(itauq1),
704 $ work(iorglq), lorglq, childinfo )
709 CALL dbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
710 $ theta, work(iphi), u2, ldu2, u1, ldu1, dum2,
711 $ 1, v1t, ldv1t, work(ib11d), work(ib11e),
712 $ work(ib12d), work(ib12e), work(ib21d),
713 $ work(ib21e), work(ib22d), work(ib22e),
714 $ work(ibbcsd), lbbcsd, childinfo )
727 CALL dlapmt( .false., p, p, u1, ldu1, iwork )
730 CALL dlapmr( .false., p, q, v1t, ldv1t, iwork )
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 dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dlapmr(FORWRD, M, N, X, LDX, K)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.
subroutine dorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB1
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dorbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
DORBDB4
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine dorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB2
subroutine dorbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB3
subroutine dlapmt(FORWRD, M, N, X, LDX, K)
DLAPMT performs a forward or backward permutation of the columns of a matrix.