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 )
LOGICAL function lsame(CA, CB)
LSAME
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 dorcsd2by1(JOBU1, JOBU2, JOBV1T, M, P, Q, X11, LDX11, X21, LDX21, THETA, U1, LDU1, U2, LDU2, V1T, LDV1T, WORK, LWORK, IWORK, INFO)
DORCSD2BY1
subroutine dorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
DORBDB2
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
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
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.