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
278 INTRINSIC int, max, min
285 wantu1 =
lsame( jobu1,
'Y' )
286 wantu2 =
lsame( jobu2,
'Y' )
287 wantv1t =
lsame( jobv1t,
'Y' )
288 lquery = lwork .EQ. -1
292 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
294 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
296 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
298 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
300 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
302 ELSE IF( wantu2 .AND. ldu2 .LT. m - p )
THEN
304 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
308 r = min( p, m-p, q, m-q )
329 IF( info .EQ. 0 )
THEN
331 ib11d = iphi + max( 1, r-1 )
332 ib11e = ib11d + max( 1, r )
333 ib12d = ib11e + max( 1, r - 1 )
334 ib12e = ib12d + max( 1, r )
335 ib21d = ib12e + max( 1, r - 1 )
336 ib21e = ib21d + max( 1, r )
337 ib22d = ib21e + max( 1, r - 1 )
338 ib22e = ib22d + max( 1, r )
339 ibbcsd = ib22e + max( 1, r - 1 )
340 itaup1 = iphi + max( 1, r-1 )
341 itaup2 = itaup1 + max( 1, p )
342 itauq1 = itaup2 + max( 1, m-p )
343 iorbdb = itauq1 + max( 1, q )
344 iorgqr = itauq1 + max( 1, q )
345 iorglq = itauq1 + max( 1, q )
347 CALL sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
348 $ 0, 0, work, -1, childinfo )
349 lorbdb = int( work(1) )
350 IF( p .GE. m-p )
THEN
351 CALL sorgqr( p, p, q, u1, ldu1, 0, work(1), -1,
353 lorgqrmin = max( 1, p )
354 lorgqropt = int( work(1) )
356 CALL sorgqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
358 lorgqrmin = max( 1, m-p )
359 lorgqropt = int( work(1) )
361 CALL sorglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
362 $ 0, work(1), -1, childinfo )
363 lorglqmin = max( 1, q-1 )
364 lorglqopt = int( work(1) )
365 CALL sbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
366 $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
367 $ 0, 0, 0, 0, 0, 0, work(1), -1, childinfo )
368 lbbcsd = int( work(1) )
369 ELSE IF( r .EQ. p )
THEN
370 CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
371 $ 0, 0, work(1), -1, childinfo )
372 lorbdb = int( work(1) )
373 IF( p-1 .GE. m-p )
THEN
374 CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
376 lorgqrmin = max( 1, p-1 )
377 lorgqropt = int( work(1) )
379 CALL sorgqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
381 lorgqrmin = max( 1, m-p )
382 lorgqropt = int( work(1) )
384 CALL sorglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
386 lorglqmin = max( 1, q )
387 lorglqopt = int( work(1) )
388 CALL sbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
389 $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
390 $ 0, 0, 0, 0, 0, 0, work(1), -1, childinfo )
391 lbbcsd = int( work(1) )
392 ELSE IF( r .EQ. m-p )
THEN
393 CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
394 $ 0, 0, work(1), -1, childinfo )
395 lorbdb = int( work(1) )
396 IF( p .GE. m-p-1 )
THEN
397 CALL sorgqr( p, p, q, u1, ldu1, 0, work(1), -1,
399 lorgqrmin = max( 1, p )
400 lorgqropt = int( work(1) )
402 CALL sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
403 $ work(1), -1, childinfo )
404 lorgqrmin = max( 1, m-p-1 )
405 lorgqropt = int( work(1) )
407 CALL sorglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
409 lorglqmin = max( 1, q )
410 lorglqopt = int( work(1) )
411 CALL sbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
412 $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
413 $ 0, 0, 0, 0, 0, 0, 0, 0, work(1), -1,
415 lbbcsd = int( work(1) )
417 CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
418 $ 0, 0, 0, work(1), -1, childinfo )
419 lorbdb = m + int( work(1) )
420 IF( p .GE. m-p )
THEN
421 CALL sorgqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
423 lorgqrmin = max( 1, p )
424 lorgqropt = int( work(1) )
426 CALL sorgqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
428 lorgqrmin = max( 1, m-p )
429 lorgqropt = int( work(1) )
431 CALL sorglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
433 lorglqmin = max( 1, q )
434 lorglqopt = int( work(1) )
435 CALL sbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
436 $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
437 $ 0, 0, 0, 0, 0, 0, 0, 0, work(1), -1,
439 lbbcsd = int( work(1) )
441 lworkmin = max( iorbdb+lorbdb-1,
442 $ iorgqr+lorgqrmin-1,
443 $ iorglq+lorglqmin-1,
445 lworkopt = max( iorbdb+lorbdb-1,
446 $ iorgqr+lorgqropt-1,
447 $ iorglq+lorglqopt-1,
450 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
454 IF( info .NE. 0 )
THEN
455 CALL xerbla(
'SORCSD2BY1', -info )
457 ELSE IF( lquery )
THEN
460 lorgqr = lwork-iorgqr+1
461 lorglq = lwork-iorglq+1
472 CALL sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
473 $ work(iphi), work(itaup1), work(itaup2),
474 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
478 IF( wantu1 .AND. p .GT. 0 )
THEN
479 CALL slacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
480 CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
481 $ lorgqr, childinfo )
483 IF( wantu2 .AND. m-p .GT. 0 )
THEN
484 CALL slacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
485 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
486 $ work(iorgqr), lorgqr, childinfo )
488 IF( wantv1t .AND. q .GT. 0 )
THEN
494 CALL slacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
496 CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
497 $ work(iorglq), lorglq, childinfo )
502 CALL sbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
503 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
504 $ work(ib11d), work(ib11e), work(ib12d),
505 $ work(ib12e), work(ib21d), work(ib21e),
506 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
512 IF( q .GT. 0 .AND. wantu2 )
THEN
514 iwork(i) = m - p - q + i
519 CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
521 ELSE IF( r .EQ. p )
THEN
527 CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
528 $ work(iphi), work(itaup1), work(itaup2),
529 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
533 IF( wantu1 .AND. p .GT. 0 )
THEN
539 CALL slacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
540 CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
541 $ work(iorgqr), lorgqr, childinfo )
543 IF( wantu2 .AND. m-p .GT. 0 )
THEN
544 CALL slacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
545 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
546 $ work(iorgqr), lorgqr, childinfo )
548 IF( wantv1t .AND. q .GT. 0 )
THEN
549 CALL slacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
550 CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
551 $ work(iorglq), lorglq, childinfo )
556 CALL sbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
557 $ work(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
558 $ work(ib11d), work(ib11e), work(ib12d),
559 $ work(ib12e), work(ib21d), work(ib21e),
560 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
566 IF( q .GT. 0 .AND. wantu2 )
THEN
568 iwork(i) = m - p - q + i
573 CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
575 ELSE IF( r .EQ. m-p )
THEN
581 CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
582 $ work(iphi), work(itaup1), work(itaup2),
583 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
587 IF( wantu1 .AND. p .GT. 0 )
THEN
588 CALL slacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
589 CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
590 $ lorgqr, childinfo )
592 IF( wantu2 .AND. m-p .GT. 0 )
THEN
598 CALL slacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
600 CALL sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
601 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
603 IF( wantv1t .AND. q .GT. 0 )
THEN
604 CALL slacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
605 CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
606 $ work(iorglq), lorglq, childinfo )
611 CALL sbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
612 $ theta, work(iphi), 0, 1, v1t, ldv1t, u2, ldu2, u1,
613 $ ldu1, work(ib11d), work(ib11e), work(ib12d),
614 $ work(ib12e), work(ib21d), work(ib21e),
615 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
629 CALL slapmt( .false., p, q, u1, ldu1, iwork )
632 CALL slapmr( .false., q, q, v1t, ldv1t, iwork )
641 CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
642 $ work(iphi), work(itaup1), work(itaup2),
643 $ work(itauq1), work(iorbdb), work(iorbdb+m),
644 $ lorbdb-m, childinfo )
648 IF( wantu1 .AND. p .GT. 0 )
THEN
649 CALL scopy( p, work(iorbdb), 1, u1, 1 )
653 CALL slacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
655 CALL sorgqr( p, p, m-q, u1, ldu1, work(itaup1),
656 $ work(iorgqr), lorgqr, childinfo )
658 IF( wantu2 .AND. m-p .GT. 0 )
THEN
659 CALL scopy( m-p, work(iorbdb+p), 1, u2, 1 )
663 CALL slacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
665 CALL sorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
666 $ work(iorgqr), lorgqr, childinfo )
668 IF( wantv1t .AND. q .GT. 0 )
THEN
669 CALL slacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
670 CALL slacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
671 $ v1t(m-q+1,m-q+1), ldv1t )
672 CALL slacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
673 $ v1t(p+1,p+1), ldv1t )
674 CALL sorglq( q, q, q, v1t, ldv1t, work(itauq1),
675 $ work(iorglq), lorglq, childinfo )
680 CALL sbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
681 $ theta, work(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
682 $ ldv1t, work(ib11d), work(ib11e), work(ib12d),
683 $ work(ib12e), work(ib21d), work(ib21e),
684 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
698 CALL slapmt( .false., p, p, u1, ldu1, iwork )
701 CALL slapmr( .false., p, q, v1t, ldv1t, iwork )
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 sorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB2
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 scopy(N, SX, INCX, SY, INCY)
SCOPY
logical function lsame(CA, CB)
LSAME
subroutine sorbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
SORBDB4
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 slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
subroutine sorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB1
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ