246 CHARACTER jobu1, jobu2, jobv1t
247 INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
251 DOUBLE PRECISION theta(*)
252 DOUBLE PRECISION u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
253 $ x11(ldx11,*), x21(ldx21,*)
260 DOUBLE PRECISION one, zero
261 parameter( one = 1.0d0, zero = 0.0d0 )
264 INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
265 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
266 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
267 $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
268 $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
269 $ lworkmin, lworkopt, r
270 LOGICAL lquery, wantu1, wantu2, wantv1t
282 INTRINSIC int, max, min
289 wantu1 =
lsame( jobu1,
'Y' )
290 wantu2 =
lsame( jobu2,
'Y' )
291 wantv1t =
lsame( jobv1t,
'Y' )
292 lquery = lwork .EQ. -1
296 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
298 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
300 ELSE IF( ldx11 .LT. max( 1, p ) )
THEN
302 ELSE IF( ldx21 .LT. max( 1, m-p ) )
THEN
304 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
306 ELSE IF( wantu2 .AND. ldu2 .LT. m - p )
THEN
308 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
312 r = min( p, m-p, q, m-q )
333 IF( info .EQ. 0 )
THEN
335 ib11d = iphi + max( 1, r-1 )
336 ib11e = ib11d + max( 1, r )
337 ib12d = ib11e + max( 1, r - 1 )
338 ib12e = ib12d + max( 1, r )
339 ib21d = ib12e + max( 1, r - 1 )
340 ib21e = ib21d + max( 1, r )
341 ib22d = ib21e + max( 1, r - 1 )
342 ib22e = ib22d + max( 1, r )
343 ibbcsd = ib22e + max( 1, r - 1 )
344 itaup1 = iphi + max( 1, r-1 )
345 itaup2 = itaup1 + max( 1, p )
346 itauq1 = itaup2 + max( 1, m-p )
347 iorbdb = itauq1 + max( 1, q )
348 iorgqr = itauq1 + max( 1, q )
349 iorglq = itauq1 + max( 1, q )
351 CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
352 $ 0, 0, work, -1, childinfo )
353 lorbdb = int( work(1) )
354 IF( p .GE. m-p )
THEN
355 CALL dorgqr( p, p, q, u1, ldu1, 0, work(1), -1,
357 lorgqrmin = max( 1, p )
358 lorgqropt = int( work(1) )
360 CALL dorgqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
362 lorgqrmin = max( 1, m-p )
363 lorgqropt = int( work(1) )
365 CALL dorglq( max(0,q-1), max(0,q-1), max(0,q-1), v1t, ldv1t,
366 $ 0, work(1), -1, childinfo )
367 lorglqmin = max( 1, q-1 )
368 lorglqopt = int( work(1) )
369 CALL dbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
370 $ 0, u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1, 0, 0,
371 $ 0, 0, 0, 0, 0, 0, work(1), -1, childinfo )
372 lbbcsd = int( work(1) )
373 ELSE IF( r .EQ. p )
THEN
374 CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
375 $ 0, 0, work(1), -1, childinfo )
376 lorbdb = int( work(1) )
377 IF( p-1 .GE. m-p )
THEN
378 CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, 0, work(1),
380 lorgqrmin = max( 1, p-1 )
381 lorgqropt = int( work(1) )
383 CALL dorgqr( m-p, m-p, q, u2, ldu2, 0, work(1), -1,
385 lorgqrmin = max( 1, m-p )
386 lorgqropt = int( work(1) )
388 CALL dorglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
390 lorglqmin = max( 1, q )
391 lorglqopt = int( work(1) )
392 CALL dbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
393 $ 0, v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2, 0, 0,
394 $ 0, 0, 0, 0, 0, 0, work(1), -1, childinfo )
395 lbbcsd = int( work(1) )
396 ELSE IF( r .EQ. m-p )
THEN
397 CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
398 $ 0, 0, work(1), -1, childinfo )
399 lorbdb = int( work(1) )
400 IF( p .GE. m-p-1 )
THEN
401 CALL dorgqr( p, p, q, u1, ldu1, 0, work(1), -1,
403 lorgqrmin = max( 1, p )
404 lorgqropt = int( work(1) )
406 CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, 0,
407 $ work(1), -1, childinfo )
408 lorgqrmin = max( 1, m-p-1 )
409 lorgqropt = int( work(1) )
411 CALL dorglq( q, q, r, v1t, ldv1t, 0, work(1), -1,
413 lorglqmin = max( 1, q )
414 lorglqopt = int( work(1) )
415 CALL dbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
416 $ theta, 0, 0, 1, v1t, ldv1t, u2, ldu2, u1, ldu1,
417 $ 0, 0, 0, 0, 0, 0, 0, 0, work(1), -1,
419 lbbcsd = int( work(1) )
421 CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, 0, 0,
422 $ 0, 0, 0, work(1), -1, childinfo )
423 lorbdb = m + int( work(1) )
424 IF( p .GE. m-p )
THEN
425 CALL dorgqr( p, p, m-q, u1, ldu1, 0, work(1), -1,
427 lorgqrmin = max( 1, p )
428 lorgqropt = int( work(1) )
430 CALL dorgqr( m-p, m-p, m-q, u2, ldu2, 0, work(1), -1,
432 lorgqrmin = max( 1, m-p )
433 lorgqropt = int( work(1) )
435 CALL dorglq( q, q, q, v1t, ldv1t, 0, work(1), -1,
437 lorglqmin = max( 1, q )
438 lorglqopt = int( work(1) )
439 CALL dbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
440 $ theta, 0, u2, ldu2, u1, ldu1, 0, 1, v1t, ldv1t,
441 $ 0, 0, 0, 0, 0, 0, 0, 0, work(1), -1,
443 lbbcsd = int( work(1) )
445 lworkmin = max( iorbdb+lorbdb-1,
446 $ iorgqr+lorgqrmin-1,
447 $ iorglq+lorglqmin-1,
449 lworkopt = max( iorbdb+lorbdb-1,
450 $ iorgqr+lorgqropt-1,
451 $ iorglq+lorglqopt-1,
454 IF( lwork .LT. lworkmin .AND. .NOT.lquery )
THEN
458 IF( info .NE. 0 )
THEN
459 CALL xerbla(
'DORCSD2BY1', -info )
461 ELSE IF( lquery )
THEN
464 lorgqr = lwork-iorgqr+1
465 lorglq = lwork-iorglq+1
476 CALL dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
477 $ work(iphi), work(itaup1), work(itaup2),
478 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
482 IF( wantu1 .AND. p .GT. 0 )
THEN
483 CALL dlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
484 CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
485 $ lorgqr, childinfo )
487 IF( wantu2 .AND. m-p .GT. 0 )
THEN
488 CALL dlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
489 CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
490 $ work(iorgqr), lorgqr, childinfo )
492 IF( wantv1t .AND. q .GT. 0 )
THEN
498 CALL dlacpy(
'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
500 CALL dorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
501 $ work(iorglq), lorglq, childinfo )
506 CALL dbbcsd( jobu1, jobu2, jobv1t,
'N',
'N', m, p, q, theta,
507 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, 0, 1,
508 $ work(ib11d), work(ib11e), work(ib12d),
509 $ work(ib12e), work(ib21d), work(ib21e),
510 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
516 IF( q .GT. 0 .AND. wantu2 )
THEN
518 iwork(i) = m - p - q + i
523 CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
525 ELSE IF( r .EQ. p )
THEN
531 CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
532 $ work(iphi), work(itaup1), work(itaup2),
533 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
537 IF( wantu1 .AND. p .GT. 0 )
THEN
543 CALL dlacpy(
'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
544 CALL dorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
545 $ work(iorgqr), lorgqr, childinfo )
547 IF( wantu2 .AND. m-p .GT. 0 )
THEN
548 CALL dlacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
549 CALL dorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
550 $ work(iorgqr), lorgqr, childinfo )
552 IF( wantv1t .AND. q .GT. 0 )
THEN
553 CALL dlacpy(
'U', p, q, x11, ldx11, v1t, ldv1t )
554 CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
555 $ work(iorglq), lorglq, childinfo )
560 CALL dbbcsd( jobv1t,
'N', jobu1, jobu2,
'T', m, q, p, theta,
561 $ work(iphi), v1t, ldv1t, 0, 1, u1, ldu1, u2, ldu2,
562 $ work(ib11d), work(ib11e), work(ib12d),
563 $ work(ib12e), work(ib21d), work(ib21e),
564 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
570 IF( q .GT. 0 .AND. wantu2 )
THEN
572 iwork(i) = m - p - q + i
577 CALL dlapmt( .false., m-p, m-p, u2, ldu2, iwork )
579 ELSE IF( r .EQ. m-p )
THEN
585 CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
586 $ work(iphi), work(itaup1), work(itaup2),
587 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
591 IF( wantu1 .AND. p .GT. 0 )
THEN
592 CALL dlacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
593 CALL dorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
594 $ lorgqr, childinfo )
596 IF( wantu2 .AND. m-p .GT. 0 )
THEN
602 CALL dlacpy(
'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
604 CALL dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
605 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
607 IF( wantv1t .AND. q .GT. 0 )
THEN
608 CALL dlacpy(
'U', m-p, q, x21, ldx21, v1t, ldv1t )
609 CALL dorglq( q, q, r, v1t, ldv1t, work(itauq1),
610 $ work(iorglq), lorglq, childinfo )
615 CALL dbbcsd(
'N', jobv1t, jobu2, jobu1,
'T', m, m-q, m-p,
616 $ theta, work(iphi), 0, 1, v1t, ldv1t, u2, ldu2, u1,
617 $ ldu1, work(ib11d), work(ib11e), work(ib12d),
618 $ work(ib12e), work(ib21d), work(ib21e),
619 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
633 CALL dlapmt( .false., p, q, u1, ldu1, iwork )
636 CALL dlapmr( .false., q, q, v1t, ldv1t, iwork )
645 CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
646 $ work(iphi), work(itaup1), work(itaup2),
647 $ work(itauq1), work(iorbdb), work(iorbdb+m),
648 $ lorbdb-m, childinfo )
652 IF( wantu1 .AND. p .GT. 0 )
THEN
653 CALL dcopy( p, work(iorbdb), 1, u1, 1 )
657 CALL dlacpy(
'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
659 CALL dorgqr( p, p, m-q, u1, ldu1, work(itaup1),
660 $ work(iorgqr), lorgqr, childinfo )
662 IF( wantu2 .AND. m-p .GT. 0 )
THEN
663 CALL dcopy( m-p, work(iorbdb+p), 1, u2, 1 )
667 CALL dlacpy(
'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
669 CALL dorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
670 $ work(iorgqr), lorgqr, childinfo )
672 IF( wantv1t .AND. q .GT. 0 )
THEN
673 CALL dlacpy(
'U', m-q, q, x21, ldx21, v1t, ldv1t )
674 CALL dlacpy(
'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
675 $ v1t(m-q+1,m-q+1), ldv1t )
676 CALL dlacpy(
'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
677 $ v1t(p+1,p+1), ldv1t )
678 CALL dorglq( q, q, q, v1t, ldv1t, work(itauq1),
679 $ work(iorglq), lorglq, childinfo )
684 CALL dbbcsd( jobu2, jobu1,
'N', jobv1t,
'N', m, m-p, m-q,
685 $ theta, work(iphi), u2, ldu2, u1, ldu1, 0, 1, v1t,
686 $ ldv1t, work(ib11d), work(ib11e), work(ib12d),
687 $ work(ib12e), work(ib21d), work(ib21e),
688 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
702 CALL dlapmt( .false., p, p, u1, ldu1, iwork )
705 CALL dlapmr( .false., p, q, v1t, ldv1t, iwork )
subroutine dorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGLQ
subroutine dorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB1
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dorbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
DORBDB4
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB2
logical function lsame(CA, CB)
LSAME
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dorbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB3
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 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.