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