259 SUBROUTINE zuncsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
260 $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
261 $ ldv1t, work, lwork, rwork, lrwork, iwork,
270 CHARACTER jobu1, jobu2, jobv1t
271 INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
273 INTEGER lrwork, lrworkmin, lrworkopt
276 DOUBLE PRECISION rwork(*)
277 DOUBLE PRECISION theta(*)
278 COMPLEX*16 u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
279 $ x11(ldx11,*), x21(ldx21,*)
287 parameter( one = (1.0d0,0.0d0), zero = (0.0d0,0.0d0) )
290 INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
291 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
292 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
293 $
j, lbbcsd, lorbdb, lorglq, lorglqmin,
294 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
295 $ lworkmin, lworkopt, r
296 LOGICAL lquery, wantu1, wantu2, wantv1t
308 INTRINSIC int, max, min
315 wantu1 =
lsame( jobu1,
'Y' )
316 wantu2 =
lsame( jobu2,
'Y' )
317 wantv1t =
lsame( jobv1t,
'Y' )
318 lquery = lwork .EQ. -1
322 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
324 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
326 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
328 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
330 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
332 ELSE IF( wantu2 .AND. ldu2 .LT. m - p )
THEN
334 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
338 r = min( p, m-p, q, m-q )
373 IF( info .EQ. 0 )
THEN
375 ib11d = iphi + max( 1, r-1 )
376 ib11e = ib11d + max( 1, r )
377 ib12d = ib11e + max( 1, r - 1 )
378 ib12e = ib12d + max( 1, r )
379 ib21d = ib12e + max( 1, r - 1 )
380 ib21e = ib21d + max( 1, r )
381 ib22d = ib21e + max( 1, r - 1 )
382 ib22e = ib22d + max( 1, r )
383 ibbcsd = ib22e + max( 1, r - 1 )
385 itaup2 = itaup1 + max( 1, p )
386 itauq1 = itaup2 + max( 1, m-p )
387 iorbdb = itauq1 + max( 1, q )
388 iorgqr = itauq1 + max( 1, q )
389 iorglq = itauq1 + max( 1, q )
391 CALL
zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
392 $ 0, 0, work, -1, childinfo )
393 lorbdb = int( work(1) )
394 IF( p .GE. m-p )
THEN
395 CALL
zungqr( p, p, q, u1, ldu1, 0, work(1), -1,
397 lorgqrmin = max( 1, p )
398 lorgqropt = int( work(1) )
400 CALL
zungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
402 lorgqrmin = max( 1, m-p )
403 lorgqropt = int( work(1) )
405 CALL
zunglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
406 $ 0, work(1), -1, childinfo )
407 lorglqmin = max( 1, q-1 )
408 lorglqopt = int( work(1) )
409 CALL
zbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
410 $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
411 $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
412 lbbcsd = int( rwork(1) )
413 ELSE IF( r .EQ. p )
THEN
414 CALL
zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
415 $ 0, 0, work(1), -1, childinfo )
416 lorbdb = int( work(1) )
417 IF( p-1 .GE. m-p )
THEN
418 CALL
zungqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
420 lorgqrmin = max( 1, p-1 )
421 lorgqropt = int( work(1) )
423 CALL
zungqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
425 lorgqrmin = max( 1, m-p )
426 lorgqropt = int( work(1) )
428 CALL
zunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
430 lorglqmin = max( 1, q )
431 lorglqopt = int( work(1) )
432 CALL
zbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
433 $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
434 $ 0, 0, 0, 0, 0, 0, rwork(1), -1, childinfo )
435 lbbcsd = int( rwork(1) )
436 ELSE IF( r .EQ. m-p )
THEN
437 CALL
zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
438 $ 0, 0, work(1), -1, childinfo )
439 lorbdb = int( work(1) )
440 IF( p .GE. m-p-1 )
THEN
441 CALL
zungqr( p, p, q, u1, ldu1, 0, work(1), -1,
443 lorgqrmin = max( 1, p )
444 lorgqropt = int( work(1) )
446 CALL
zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
447 $ work(1), -1, childinfo )
448 lorgqrmin = max( 1, m-p-1 )
449 lorgqropt = int( work(1) )
451 CALL
zunglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
453 lorglqmin = max( 1, q )
454 lorglqopt = int( work(1) )
455 CALL
zbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
456 $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
457 $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
459 lbbcsd = int( rwork(1) )
461 CALL
zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
462 $ 0, 0, 0, work(1), -1, childinfo )
463 lorbdb = m + int( work(1) )
464 IF( p .GE. m-p )
THEN
465 CALL
zungqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
467 lorgqrmin = max( 1, p )
468 lorgqropt = int( work(1) )
470 CALL
zungqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
472 lorgqrmin = max( 1, m-p )
473 lorgqropt = int( work(1) )
475 CALL
zunglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
477 lorglqmin = max( 1, q )
478 lorglqopt = int( work(1) )
479 CALL
zbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
480 $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
481 $ 0, 0, 0, 0, 0, 0, 0, 0, rwork(1), -1,
483 lbbcsd = int( rwork(1) )
485 lrworkmin = ibbcsd+lbbcsd-1
486 lrworkopt = lrworkmin
488 lworkmin = max( iorbdb+lorbdb-1,
489 $ iorgqr+lorgqrmin-1,
490 $ iorglq+lorglqmin-1 )
491 lworkopt = max( iorbdb+lorbdb-1,
492 $ iorgqr+lorgqropt-1,
493 $ iorglq+lorglqopt-1 )
495 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
499 IF( info .NE. 0 )
THEN
500 CALL
xerbla(
'ZUNCSD2BY1', -info )
502 ELSE IF( lquery )
THEN
505 lorgqr = lwork-iorgqr+1
506 lorglq = lwork-iorglq+1
517 CALL
zunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
518 $ rwork(iphi), work(itaup1), work(itaup2),
519 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
523 IF( wantu1 .AND. p .GT. 0 )
THEN
524 CALL
zlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
525 CALL
zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
526 $ lorgqr, childinfo )
528 IF( wantu2 .AND. m-p .GT. 0 )
THEN
529 CALL
zlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
530 CALL
zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
531 $ work(iorgqr), lorgqr, childinfo )
533 IF( wantv1t .AND. q .GT. 0 )
THEN
539 CALL
zlacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
541 CALL
zunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
542 $ work(iorglq), lorglq, childinfo )
547 CALL
zbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
548 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
549 $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
550 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
551 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
557 IF( q .GT. 0 .AND. wantu2 )
THEN
559 iwork(i) = m - p - q + i
564 CALL
zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
566 ELSE IF( r .EQ. p )
THEN
572 CALL
zunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
573 $ rwork(iphi), work(itaup1), work(itaup2),
574 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
578 IF( wantu1 .AND. p .GT. 0 )
THEN
584 CALL
zlacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
585 CALL
zungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
586 $ work(iorgqr), lorgqr, childinfo )
588 IF( wantu2 .AND. m-p .GT. 0 )
THEN
589 CALL
zlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
590 CALL
zungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
591 $ work(iorgqr), lorgqr, childinfo )
593 IF( wantv1t .AND. q .GT. 0 )
THEN
594 CALL
zlacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
595 CALL
zunglq( q, q, r, v1t, ldv1t, work(itauq1),
596 $ work(iorglq), lorglq, childinfo )
601 CALL
zbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
602 $ rwork(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
603 $ rwork(ib11d), rwork(ib11e), rwork(ib12d),
604 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
605 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
611 IF( q .GT. 0 .AND. wantu2 )
THEN
613 iwork(i) = m - p - q + i
618 CALL
zlapmt( .false., m-p, m-p, u2, ldu2, iwork )
620 ELSE IF( r .EQ. m-p )
THEN
626 CALL
zunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
627 $ rwork(iphi), work(itaup1), work(itaup2),
628 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
632 IF( wantu1 .AND. p .GT. 0 )
THEN
633 CALL
zlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
634 CALL
zungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
635 $ lorgqr, childinfo )
637 IF( wantu2 .AND. m-p .GT. 0 )
THEN
643 CALL
zlacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
645 CALL
zungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
646 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
648 IF( wantv1t .AND. q .GT. 0 )
THEN
649 CALL
zlacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
650 CALL
zunglq( q, q, r, v1t, ldv1t, work(itauq1),
651 $ work(iorglq), lorglq, childinfo )
656 CALL
zbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
657 $ theta, rwork(iphi), 0, 1, v1t, ldv1t, u2, ldu2,
658 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
659 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
660 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
661 $ rwork(ibbcsd), lbbcsd, childinfo )
674 CALL
zlapmt( .false., p, q, u1, ldu1, iwork )
677 CALL
zlapmr( .false., q, q, v1t, ldv1t, iwork )
686 CALL
zunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
687 $ rwork(iphi), work(itaup1), work(itaup2),
688 $ work(itauq1), work(iorbdb), work(iorbdb+m),
689 $ lorbdb-m, childinfo )
693 IF( wantu1 .AND. p .GT. 0 )
THEN
694 CALL
zcopy( p, work(iorbdb), 1, u1, 1 )
698 CALL
zlacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
700 CALL
zungqr( p, p, m-q, u1, ldu1, work(itaup1),
701 $ work(iorgqr), lorgqr, childinfo )
703 IF( wantu2 .AND. m-p .GT. 0 )
THEN
704 CALL
zcopy( m-p, work(iorbdb+p), 1, u2, 1 )
708 CALL
zlacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
710 CALL
zungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
711 $ work(iorgqr), lorgqr, childinfo )
713 IF( wantv1t .AND. q .GT. 0 )
THEN
714 CALL
zlacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
715 CALL
zlacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
716 $ v1t(m-q+1,m-q+1), ldv1t )
717 CALL
zlacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
718 $ v1t(p+1,p+1), ldv1t )
719 CALL
zunglq( q, q, q, v1t, ldv1t, work(itauq1),
720 $ work(iorglq), lorglq, childinfo )
725 CALL
zbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
726 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
727 $ ldv1t, rwork(ib11d), rwork(ib11e), rwork(ib12d),
728 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
729 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
743 CALL
zlapmt( .false., p, p, u1, ldu1, iwork )
746 CALL
zlapmr( .false., p, q, v1t, ldv1t, iwork )