236 SUBROUTINE dorcsd2by1( JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11,
237 $ x21, ldx21, theta, u1, ldu1, u2, ldu2, v1t,
238 $ ldv1t, work, lwork, iwork, info )
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 )