293 RECURSIVE SUBROUTINE sorcsd( JOBU1, JOBU2, JOBV1T, JOBV2T,
295 $ SIGNS, M, P, Q, X11, LDX11, X12,
296 $ LDX12, X21, LDX21, X22, LDX22, THETA,
297 $ U1, LDU1, U2, LDU2, V1T, LDV1T, V2T,
298 $ LDV2T, WORK, LWORK, IWORK, INFO )
305 CHARACTER jobu1, jobu2, jobv1t, jobv2t, signs, trans
306 INTEGER info, ldu1, ldu2, ldv1t, ldv2t, ldx11, ldx12,
307 $ ldx21, ldx22, lwork, m, p, q
312 REAL u1( ldu1, * ), u2( ldu2, * ), v1t( ldv1t, * ),
313 $ v2t( ldv2t, * ), work( * ), x11( ldx11, * ),
314 $ x12( ldx12, * ), x21( ldx21, * ), x22( ldx22,
322 PARAMETER ( one = 1.0e+0,
329 CHARACTER transt, signst
330 INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
331 $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
332 $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
333 $ itauq2, j, lbbcsdwork, lbbcsdworkmin,
334 $ lbbcsdworkopt, lorbdbwork, lorbdbworkmin,
335 $ lorbdbworkopt, lorglqwork, lorglqworkmin,
336 $ lorglqworkopt, lorgqrwork, lorgqrworkmin,
337 $ lorgqrworkopt, lworkmin, lworkopt
338 LOGICAL colmajor, defaultsigns, lquery, wantu1, wantu2,
350 INTRINSIC int, max, min
357 wantu1 =
lsame( jobu1,
'Y' )
358 wantu2 =
lsame( jobu2,
'Y' )
359 wantv1t =
lsame( jobv1t,
'Y' )
360 wantv2t =
lsame( jobv2t,
'Y' )
361 colmajor = .NOT.
lsame( trans,
'T' )
362 defaultsigns = .NOT.
lsame( signs,
'O' )
363 lquery = lwork .EQ. -1
366 ELSE IF( p .LT. 0 .OR. p .GT. m )
THEN
368 ELSE IF( q .LT. 0 .OR. q .GT. m )
THEN
370 ELSE IF ( colmajor .AND. ldx11 .LT. max( 1, p ) )
THEN
372 ELSE IF (.NOT. colmajor .AND. ldx11 .LT. max( 1, q ) )
THEN
374 ELSE IF (colmajor .AND. ldx12 .LT. max( 1, p ) )
THEN
376 ELSE IF (.NOT. colmajor .AND. ldx12 .LT. max( 1, m-q ) )
THEN
378 ELSE IF (colmajor .AND. ldx21 .LT. max( 1, m-p ) )
THEN
380 ELSE IF (.NOT. colmajor .AND. ldx21 .LT. max( 1, q ) )
THEN
382 ELSE IF (colmajor .AND. ldx22 .LT. max( 1, m-p ) )
THEN
384 ELSE IF (.NOT. colmajor .AND. ldx22 .LT. max( 1, m-q ) )
THEN
386 ELSE IF( wantu1 .AND. ldu1 .LT. p )
THEN
388 ELSE IF( wantu2 .AND. ldu2 .LT. m-p )
THEN
390 ELSE IF( wantv1t .AND. ldv1t .LT. q )
THEN
392 ELSE IF( wantv2t .AND. ldv2t .LT. m-q )
THEN
398 IF( info .EQ. 0 .AND. min( p, m-p ) .LT. min( q, m-q ) )
THEN
404 IF( defaultsigns )
THEN
409 CALL sorcsd( jobv1t, jobv2t, jobu1, jobu2, transt, signst,
411 $ q, p, x11, ldx11, x21, ldx21, x12, ldx12, x22,
412 $ ldx22, theta, v1t, ldv1t, v2t, ldv2t, u1, ldu1,
413 $ u2, ldu2, work, lwork, iwork, info )
420 IF( info .EQ. 0 .AND. m-q .LT. q )
THEN
421 IF( defaultsigns )
THEN
426 CALL sorcsd( jobu2, jobu1, jobv2t, jobv1t, trans, signst, m,
427 $ m-p, m-q, x22, ldx22, x21, ldx21, x12, ldx12, x11,
428 $ ldx11, theta, u2, ldu2, u1, ldu1, v2t, ldv2t, v1t,
429 $ ldv1t, work, lwork, iwork, info )
435 IF( info .EQ. 0 )
THEN
438 itaup1 = iphi + max( 1, q - 1 )
439 itaup2 = itaup1 + max( 1, p )
440 itauq1 = itaup2 + max( 1, m - p )
441 itauq2 = itauq1 + max( 1, q )
442 iorgqr = itauq2 + max( 1, m - q )
443 CALL sorgqr( m-q, m-q, m-q, dummy, max(1,m-q), dummy, work,
446 lorgqrworkopt = int( work(1) )
447 lorgqrworkmin = max( 1, m - q )
448 iorglq = itauq2 + max( 1, m - q )
449 CALL sorglq( m-q, m-q, m-q, dummy, max(1,m-q), dummy, work,
452 lorglqworkopt = int( work(1) )
453 lorglqworkmin = max( 1, m - q )
454 iorbdb = itauq2 + max( 1, m - q )
455 CALL sorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
456 $ x21, ldx21, x22, ldx22, dummy, dummy, dummy, dummy, dummy,
457 $ dummy,work,-1,childinfo )
458 lorbdbworkopt = int( work(1) )
459 lorbdbworkmin = lorbdbworkopt
460 ib11d = itauq2 + max( 1, m - q )
461 ib11e = ib11d + max( 1, q )
462 ib12d = ib11e + max( 1, q - 1 )
463 ib12e = ib12d + max( 1, q )
464 ib21d = ib12e + max( 1, q - 1 )
465 ib21e = ib21d + max( 1, q )
466 ib22d = ib21e + max( 1, q - 1 )
467 ib22e = ib22d + max( 1, q )
468 ibbcsd = ib22e + max( 1, q - 1 )
469 CALL sbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
470 $ dummy, dummy, u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
471 $ ldv2t, dummy, dummy, dummy, dummy, dummy, dummy,
472 $ dummy, dummy, work, -1, childinfo )
473 lbbcsdworkopt = int( work(1) )
474 lbbcsdworkmin = lbbcsdworkopt
475 lworkopt = max( iorgqr + lorgqrworkopt, iorglq + lorglqworkopt,
476 $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkopt ) - 1
477 lworkmin = max( iorgqr + lorgqrworkmin, iorglq + lorglqworkmin,
478 $ iorbdb + lorbdbworkopt, ibbcsd + lbbcsdworkmin ) - 1
479 work(1) = real( max(lworkopt,lworkmin) )
481 IF( lwork .LT. lworkmin .AND. .NOT. lquery )
THEN
484 lorgqrwork = lwork - iorgqr + 1
485 lorglqwork = lwork - iorglq + 1
486 lorbdbwork = lwork - iorbdb + 1
487 lbbcsdwork = lwork - ibbcsd + 1
493 IF( info .NE. 0 )
THEN
494 CALL xerbla(
'SORCSD', -info )
496 ELSE IF( lquery )
THEN
502 CALL sorbdb( trans, signs, m, p, q, x11, ldx11, x12, ldx12,
504 $ ldx21, x22, ldx22, theta, work(iphi), work(itaup1),
505 $ work(itaup2), work(itauq1), work(itauq2),
506 $ work(iorbdb), lorbdbwork, childinfo )
511 IF( wantu1 .AND. p .GT. 0 )
THEN
512 CALL slacpy(
'L', p, q, x11, ldx11, u1, ldu1 )
513 CALL sorgqr( p, p, q, u1, ldu1, work(itaup1),
517 IF( wantu2 .AND. m-p .GT. 0 )
THEN
518 CALL slacpy(
'L', m-p, q, x21, ldx21, u2, ldu2 )
519 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
520 $ work(iorgqr), lorgqrwork, info )
522 IF( wantv1t .AND. q .GT. 0 )
THEN
523 CALL slacpy(
'U', q-1, q-1, x11(1,2), ldx11, v1t(2,2),
530 CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t,
532 $ work(iorglq), lorglqwork, info )
534 IF( wantv2t .AND. m-q .GT. 0 )
THEN
535 CALL slacpy(
'U', p, m-q, x12, ldx12, v2t, ldv2t )
536 CALL slacpy(
'U', m-p-q, m-p-q, x22(q+1,p+1), ldx22,
537 $ v2t(p+1,p+1), ldv2t )
538 CALL sorglq( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
539 $ work(iorglq), lorglqwork, info )
542 IF( wantu1 .AND. p .GT. 0 )
THEN
543 CALL slacpy(
'U', q, p, x11, ldx11, u1, ldu1 )
544 CALL sorglq( p, p, q, u1, ldu1, work(itaup1),
548 IF( wantu2 .AND. m-p .GT. 0 )
THEN
549 CALL slacpy(
'U', q, m-p, x21, ldx21, u2, ldu2 )
550 CALL sorglq( m-p, m-p, q, u2, ldu2, work(itaup2),
551 $ work(iorglq), lorglqwork, info )
553 IF( wantv1t .AND. q .GT. 0 )
THEN
554 CALL slacpy(
'L', q-1, q-1, x11(2,1), ldx11, v1t(2,2),
561 CALL sorgqr( q-1, q-1, q-1, v1t(2,2), ldv1t,
563 $ work(iorgqr), lorgqrwork, info )
565 IF( wantv2t .AND. m-q .GT. 0 )
THEN
566 CALL slacpy(
'L', m-q, p, x12, ldx12, v2t, ldv2t )
567 CALL slacpy(
'L', m-p-q, m-p-q, x22(p+1,q+1), ldx22,
568 $ v2t(p+1,p+1), ldv2t )
569 CALL sorgqr( m-q, m-q, m-q, v2t, ldv2t, work(itauq2),
570 $ work(iorgqr), lorgqrwork, info )
576 CALL sbbcsd( jobu1, jobu2, jobv1t, jobv2t, trans, m, p, q,
578 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, v2t,
579 $ ldv2t, work(ib11d), work(ib11e), work(ib12d),
580 $ work(ib12e), work(ib21d), work(ib21e), work(ib22d),
581 $ work(ib22e), work(ibbcsd), lbbcsdwork, info )
588 IF( q .GT. 0 .AND. wantu2 )
THEN
590 iwork(i) = m - p - q + i
596 CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
598 CALL slapmr( .false., m-p, m-p, u2, ldu2, iwork )
601 IF( m .GT. 0 .AND. wantv2t )
THEN
603 iwork(i) = m - p - q + i
608 IF( .NOT. colmajor )
THEN
609 CALL slapmt( .false., m-q, m-q, v2t, ldv2t, iwork )
611 CALL slapmr( .false., m-q, m-q, v2t, ldv2t, iwork )