242 CHARACTER jobu1, jobu2, jobv1t
243 INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
248 REAL u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
249 $ x11(ldx11,*), x21(ldx21,*)
257 parameter ( one = 1.0e0, zero = 0.0e0 )
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 REAL 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 sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
355 $ dum1, dum1, dum1, dum1, work, -1,
357 lorbdb = int( work(1) )
358 IF( wantu1 .AND. p .GT. 0 )
THEN
359 CALL sorgqr( 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 sorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1), -1,
367 lorgqrmin = max( lorgqrmin, m-p )
368 lorgqropt = max( lorgqropt, int( work(1) ) )
370 IF( wantv1t .AND. q .GT. 0 )
THEN
371 CALL sorglq( 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 sbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
377 $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t, dum2,
378 $ 1, dum1, dum1, dum1, dum1, dum1,
379 $ dum1, dum1, dum1, work(1), -1, childinfo
381 lbbcsd = int( work(1) )
382 ELSE IF( r .EQ. p )
THEN
383 CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
384 $ dum1, dum1, dum1, dum1, work(1), -1,
386 lorbdb = int( work(1) )
387 IF( wantu1 .AND. p .GT. 0 )
THEN
388 CALL sorgqr( 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 sorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1), -1,
396 lorgqrmin = max( lorgqrmin, m-p )
397 lorgqropt = max( lorgqropt, int( work(1) ) )
399 IF( wantv1t .AND. q .GT. 0 )
THEN
400 CALL sorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
402 lorglqmin = max( lorglqmin, q )
403 lorglqopt = max( lorglqopt, int( work(1) ) )
405 CALL sbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
406 $ dum1, v1t, ldv1t, dum2, 1, u1, ldu1, u2,
407 $ ldu2, dum1, dum1, dum1, dum1, dum1,
408 $ dum1, dum1, dum1, work(1), -1, childinfo
410 lbbcsd = int( work(1) )
411 ELSE IF( r .EQ. m-p )
THEN
412 CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
413 $ dum1, dum1, dum1, dum1, work(1), -1,
415 lorbdb = int( work(1) )
416 IF( wantu1 .AND. p .GT. 0 )
THEN
417 CALL sorgqr( 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 sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, dum1,
424 $ 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 sorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
431 lorglqmin = max( lorglqmin, q )
432 lorglqopt = max( lorglqopt, int( work(1) ) )
434 CALL sbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
435 $ theta, dum1, dum2, 1, v1t, ldv1t, u2, ldu2,
436 $ u1, ldu1, dum1, dum1, dum1, dum1,
437 $ dum1, dum1, dum1, dum1, work(1), -1,
439 lbbcsd = int( work(1) )
441 CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
442 $ dum1, dum1, dum1, dum1, dum1,
443 $ work(1), -1, childinfo )
444 lorbdb = m + int( work(1) )
445 IF( wantu1 .AND. p .GT. 0 )
THEN
446 CALL sorgqr( 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 sorgqr( 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 sorglq( q, q, q, v1t, ldv1t, dum1, work(1), -1,
460 lorglqmin = max( lorglqmin, q )
461 lorglqopt = max( lorglqopt, int( work(1) ) )
463 CALL sbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
464 $ theta, dum1, u2, ldu2, u1, ldu1, dum2, 1,
465 $ v1t, ldv1t, dum1, dum1, dum1, dum1,
466 $ dum1, dum1, dum1, dum1, work(1), -1,
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(
'SORCSD2BY1', -info )
486 ELSE IF( lquery )
THEN
489 lorgqr = lwork-iorgqr+1
490 lorglq = lwork-iorglq+1
501 CALL sorbdb1( 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 slacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
509 CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
510 $ lorgqr, childinfo )
512 IF( wantu2 .AND. m-p .GT. 0 )
THEN
513 CALL slacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
514 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
515 $ work(iorgqr), lorgqr, childinfo )
517 IF( wantv1t .AND. q .GT. 0 )
THEN
523 CALL slacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
525 CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
526 $ work(iorglq), lorglq, childinfo )
531 CALL sbbcsd( 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), work(ib12d),
534 $ work(ib12e), work(ib21d), work(ib21e),
535 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
541 IF( q .GT. 0 .AND. wantu2 )
THEN
543 iwork(i) = m - p - q + i
548 CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
550 ELSE IF( r .EQ. p )
THEN
556 CALL sorbdb2( 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 slacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
569 CALL sorgqr( 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 slacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
574 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
575 $ work(iorgqr), lorgqr, childinfo )
577 IF( wantv1t .AND. q .GT. 0 )
THEN
578 CALL slacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
579 CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
580 $ work(iorglq), lorglq, childinfo )
585 CALL sbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
586 $ work(iphi), v1t, ldv1t, dum1, 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 slapmt( .false., m-p, m-p, u2, ldu2, iwork )
604 ELSE IF( r .EQ. m-p )
THEN
610 CALL sorbdb3( 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 slacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
618 CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
619 $ lorgqr, childinfo )
621 IF( wantu2 .AND. m-p .GT. 0 )
THEN
627 CALL slacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
629 CALL sorgqr( 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 slacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
634 CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
635 $ work(iorglq), lorglq, childinfo )
640 CALL sbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
641 $ theta, work(iphi), dum1, 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 slapmt( .false., p, q, u1, ldu1, iwork )
661 CALL slapmr( .false., q, q, v1t, ldv1t, iwork )
670 CALL sorbdb4( 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 scopy( p, work(iorbdb), 1, u1, 1 )
682 CALL slacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
684 CALL sorgqr( p, p, m-q, u1, ldu1, work(itaup1),
685 $ work(iorgqr), lorgqr, childinfo )
687 IF( wantu2 .AND. m-p .GT. 0 )
THEN
688 CALL scopy( m-p, work(iorbdb+p), 1, u2, 1 )
692 CALL slacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
694 CALL sorgqr( 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 slacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
699 CALL slacpy(
'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 slacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
702 $ v1t(p+1,p+1), ldv1t )
703 CALL sorglq( q, q, q, v1t, ldv1t, work(itauq1),
704 $ work(iorglq), lorglq, childinfo )
709 CALL sbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
710 $ theta, work(iphi), u2, ldu2, u1, ldu1, dum1, 1,
711 $ v1t, ldv1t, work(ib11d), work(ib11e), work(ib12d),
712 $ work(ib12e), work(ib21d), work(ib21e),
713 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
727 CALL slapmt( .false., p, p, u1, ldu1, iwork )
730 CALL slapmr( .false., p, q, v1t, ldv1t, iwork )
subroutine sorbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
SORBDB4
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 sorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB1
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
logical function lsame(CA, CB)
LSAME
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sorbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB3
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB2
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY