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 )
subroutine slapmr(FORWRD, M, N, X, LDX, K)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
LOGICAL function lsame(CA, CB)
LSAME
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
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.
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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
subroutine sorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
SORCSD2BY1