LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ cuncsd2by1()

subroutine cuncsd2by1 ( character jobu1,
character jobu2,
character jobv1t,
integer m,
integer p,
integer q,
complex, dimension(ldx11,*) x11,
integer ldx11,
complex, dimension(ldx21,*) x21,
integer ldx21,
real, dimension(*) theta,
complex, dimension(ldu1,*) u1,
integer ldu1,
complex, dimension(ldu2,*) u2,
integer ldu2,
complex, dimension(ldv1t,*) v1t,
integer ldv1t,
complex, dimension(*) work,
integer lwork,
real, dimension(*) rwork,
integer lrwork,
integer, dimension(*) iwork,
integer info )

CUNCSD2BY1

Download CUNCSD2BY1 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> CUNCSD2BY1 computes the CS decomposition of an M-by-Q matrix X with
!> orthonormal columns that has been partitioned into a 2-by-1 block
!> structure:
!>
!>                                [  I1 0  0 ]
!>                                [  0  C  0 ]
!>          [ X11 ]   [ U1 |    ] [  0  0  0 ]
!>      X = [-----] = [---------] [----------] V1**T .
!>          [ X21 ]   [    | U2 ] [  0  0  0 ]
!>                                [  0  S  0 ]
!>                                [  0  0  I2]
!>
!> X11 is P-by-Q. The unitary matrices U1, U2, and V1 are P-by-P,
!> (M-P)-by-(M-P), and Q-by-Q, respectively. C and S are R-by-R
!> nonnegative diagonal matrices satisfying C^2 + S^2 = I, in which
!> R = MIN(P,M-P,Q,M-Q). I1 is a K1-by-K1 identity matrix and I2 is a
!> K2-by-K2 identity matrix, where K1 = MAX(Q+P-M,0), K2 = MAX(Q-P,0).
!>
!>
Parameters
[in]JOBU1
!>          JOBU1 is CHARACTER
!>          = 'Y':      U1 is computed;
!>          otherwise:  U1 is not computed.
!>
[in]JOBU2
!>          JOBU2 is CHARACTER
!>          = 'Y':      U2 is computed;
!>          otherwise:  U2 is not computed.
!>
[in]JOBV1T
!>          JOBV1T is CHARACTER
!>          = 'Y':      V1T is computed;
!>          otherwise:  V1T is not computed.
!>
[in]M
!>          M is INTEGER
!>          The number of rows in X.
!>
[in]P
!>          P is INTEGER
!>          The number of rows in X11. 0 <= P <= M.
!>
[in]Q
!>          Q is INTEGER
!>          The number of columns in X11 and X21. 0 <= Q <= M.
!>
[in,out]X11
!>          X11 is COMPLEX array, dimension (LDX11,Q)
!>          On entry, part of the unitary matrix whose CSD is desired.
!>
[in]LDX11
!>          LDX11 is INTEGER
!>          The leading dimension of X11. LDX11 >= MAX(1,P).
!>
[in,out]X21
!>          X21 is COMPLEX array, dimension (LDX21,Q)
!>          On entry, part of the unitary matrix whose CSD is desired.
!>
[in]LDX21
!>          LDX21 is INTEGER
!>          The leading dimension of X21. LDX21 >= MAX(1,M-P).
!>
[out]THETA
!>          THETA is REAL array, dimension (R), in which R =
!>          MIN(P,M-P,Q,M-Q).
!>          C = DIAG( COS(THETA(1)), ... , COS(THETA(R)) ) and
!>          S = DIAG( SIN(THETA(1)), ... , SIN(THETA(R)) ).
!>
[out]U1
!>          U1 is COMPLEX array, dimension (P)
!>          If JOBU1 = 'Y', U1 contains the P-by-P unitary matrix U1.
!>
[in]LDU1
!>          LDU1 is INTEGER
!>          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
!>          MAX(1,P).
!>
[out]U2
!>          U2 is COMPLEX array, dimension (M-P)
!>          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) unitary
!>          matrix U2.
!>
[in]LDU2
!>          LDU2 is INTEGER
!>          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
!>          MAX(1,M-P).
!>
[out]V1T
!>          V1T is COMPLEX array, dimension (Q)
!>          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix unitary
!>          matrix V1**T.
!>
[in]LDV1T
!>          LDV1T is INTEGER
!>          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
!>          MAX(1,Q).
!>
[out]WORK
!>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!>
[in]LWORK
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK and RWORK
!>          arrays, returns this value as the first entry of the WORK
!>          and RWORK array, respectively, and no error message related
!>          to LWORK or LRWORK is issued by XERBLA.
!>
[out]RWORK
!>          RWORK is REAL array, dimension (MAX(1,LRWORK))
!>          On exit, if INFO = 0, RWORK(1) returns the optimal LRWORK.
!>          If INFO > 0 on exit, RWORK(2:R) contains the values PHI(1),
!>          ..., PHI(R-1) that, together with THETA(1), ..., THETA(R),
!>          define the matrix in intermediate bidiagonal-block form
!>          remaining after nonconvergence. INFO specifies the number
!>          of nonzero PHI's.
!>
[in]LRWORK
!>          LRWORK is INTEGER
!>          The dimension of the array RWORK.
!>
!>          If LRWORK=-1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK and RWORK
!>          arrays, returns this value as the first entry of the WORK
!>          and RWORK array, respectively, and no error message related
!>          to LWORK or LRWORK is issued by XERBLA.
!>
[out]IWORK
!>          IWORK is INTEGER array, dimension (M-MIN(P,M-P,Q,M-Q))
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  CBBCSD did not converge. See the description of WORK
!>                above for details.
!>
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 251 of file cuncsd2by1.f.

256*
257* -- LAPACK computational routine --
258* -- LAPACK is a software package provided by Univ. of Tennessee, --
259* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
260*
261* .. Scalar Arguments ..
262 CHARACTER JOBU1, JOBU2, JOBV1T
263 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
264 $ M, P, Q
265 INTEGER LRWORK, LRWORKMIN, LRWORKOPT
266* ..
267* .. Array Arguments ..
268 REAL RWORK(*)
269 REAL THETA(*)
270 COMPLEX U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
271 $ X11(LDX11,*), X21(LDX21,*)
272 INTEGER IWORK(*)
273* ..
274*
275* =====================================================================
276*
277* .. Parameters ..
278 COMPLEX ONE, ZERO
279 parameter( one = (1.0e0,0.0e0), zero = (0.0e0,0.0e0) )
280* ..
281* .. Local Scalars ..
282 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
283 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
284 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
285 $ J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN,
286 $ LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT,
287 $ LWORKMIN, LWORKOPT, R
288 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
289* ..
290* .. Local Arrays ..
291 REAL DUM( 1 )
292 COMPLEX CDUM( 1, 1 )
293* ..
294* .. External Subroutines ..
295 EXTERNAL cbbcsd, ccopy, clacpy, clapmr, clapmt,
296 $ cunbdb1,
297 $ cunbdb2, cunbdb3, cunbdb4, cunglq, cungqr,
298 $ xerbla
299* ..
300* .. External Functions ..
301 LOGICAL LSAME
302 REAL SROUNDUP_LWORK
303 EXTERNAL lsame, sroundup_lwork
304* ..
305* .. Intrinsic Function ..
306 INTRINSIC int, max, min
307* ..
308* .. Executable Statements ..
309*
310* Test input arguments
311*
312 info = 0
313 wantu1 = lsame( jobu1, 'Y' )
314 wantu2 = lsame( jobu2, 'Y' )
315 wantv1t = lsame( jobv1t, 'Y' )
316 lquery = ( lwork.EQ.-1 ) .OR. ( lrwork.EQ.-1 )
317*
318 IF( m .LT. 0 ) THEN
319 info = -4
320 ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
321 info = -5
322 ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
323 info = -6
324 ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
325 info = -8
326 ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
327 info = -10
328 ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) ) THEN
329 info = -13
330 ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) ) THEN
331 info = -15
332 ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) ) THEN
333 info = -17
334 END IF
335*
336 r = min( p, m-p, q, m-q )
337*
338* Compute workspace
339*
340* WORK layout:
341* |-----------------------------------------|
342* | LWORKOPT (1) |
343* |-----------------------------------------|
344* | TAUP1 (MAX(1,P)) |
345* | TAUP2 (MAX(1,M-P)) |
346* | TAUQ1 (MAX(1,Q)) |
347* |-----------------------------------------|
348* | CUNBDB WORK | CUNGQR WORK | CUNGLQ WORK |
349* | | | |
350* | | | |
351* | | | |
352* | | | |
353* |-----------------------------------------|
354* RWORK layout:
355* |------------------|
356* | LRWORKOPT (1) |
357* |------------------|
358* | PHI (MAX(1,R-1)) |
359* |------------------|
360* | B11D (R) |
361* | B11E (R-1) |
362* | B12D (R) |
363* | B12E (R-1) |
364* | B21D (R) |
365* | B21E (R-1) |
366* | B22D (R) |
367* | B22E (R-1) |
368* | CBBCSD RWORK |
369* |------------------|
370*
371 IF( info .EQ. 0 ) THEN
372 iphi = 2
373 ib11d = iphi + max( 1, r-1 )
374 ib11e = ib11d + max( 1, r )
375 ib12d = ib11e + max( 1, r - 1 )
376 ib12e = ib12d + max( 1, r )
377 ib21d = ib12e + max( 1, r - 1 )
378 ib21e = ib21d + max( 1, r )
379 ib22d = ib21e + max( 1, r - 1 )
380 ib22e = ib22d + max( 1, r )
381 ibbcsd = ib22e + max( 1, r - 1 )
382 itaup1 = 2
383 itaup2 = itaup1 + max( 1, p )
384 itauq1 = itaup2 + max( 1, m-p )
385 iorbdb = itauq1 + max( 1, q )
386 iorgqr = itauq1 + max( 1, q )
387 iorglq = itauq1 + max( 1, q )
388 lorgqrmin = 1
389 lorgqropt = 1
390 lorglqmin = 1
391 lorglqopt = 1
392 IF( r .EQ. q ) THEN
393 CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
394 $ dum, cdum, cdum, cdum, work, -1,
395 $ childinfo )
396 lorbdb = int( work(1) )
397 IF( wantu1 .AND. p .GT. 0 ) THEN
398 CALL cungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
399 $ childinfo )
400 lorgqrmin = max( lorgqrmin, p )
401 lorgqropt = max( lorgqropt, int( work(1) ) )
402 ENDIF
403 IF( wantu2 .AND. m-p .GT. 0 ) THEN
404 CALL cungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
405 $ childinfo )
406 lorgqrmin = max( lorgqrmin, m-p )
407 lorgqropt = max( lorgqropt, int( work(1) ) )
408 END IF
409 IF( wantv1t .AND. q .GT. 0 ) THEN
410 CALL cunglq( q-1, q-1, q-1, v1t, ldv1t,
411 $ cdum, work(1), -1, childinfo )
412 lorglqmin = max( lorglqmin, q-1 )
413 lorglqopt = max( lorglqopt, int( work(1) ) )
414 END IF
415 CALL cbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q,
416 $ theta,
417 $ dum(1), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
418 $ 1, dum, dum, dum, dum, dum, dum, dum, dum,
419 $ rwork(1), -1, childinfo )
420 lbbcsd = int( rwork(1) )
421 ELSE IF( r .EQ. p ) THEN
422 CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
423 $ dum,
424 $ cdum, cdum, cdum, work(1), -1, childinfo )
425 lorbdb = int( work(1) )
426 IF( wantu1 .AND. p .GT. 0 ) THEN
427 CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, cdum,
428 $ work(1),
429 $ -1, childinfo )
430 lorgqrmin = max( lorgqrmin, p-1 )
431 lorgqropt = max( lorgqropt, int( work(1) ) )
432 END IF
433 IF( wantu2 .AND. m-p .GT. 0 ) THEN
434 CALL cungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
435 $ childinfo )
436 lorgqrmin = max( lorgqrmin, m-p )
437 lorgqropt = max( lorgqropt, int( work(1) ) )
438 END IF
439 IF( wantv1t .AND. q .GT. 0 ) THEN
440 CALL cunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
441 $ childinfo )
442 lorglqmin = max( lorglqmin, q )
443 lorglqopt = max( lorglqopt, int( work(1) ) )
444 END IF
445 CALL cbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p,
446 $ theta,
447 $ dum, v1t, ldv1t, cdum, 1, u1, ldu1, u2, ldu2,
448 $ dum, dum, dum, dum, dum, dum, dum, dum,
449 $ rwork(1), -1, childinfo )
450 lbbcsd = int( rwork(1) )
451 ELSE IF( r .EQ. m-p ) THEN
452 CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
453 $ dum,
454 $ cdum, cdum, cdum, work(1), -1, childinfo )
455 lorbdb = int( work(1) )
456 IF( wantu1 .AND. p .GT. 0 ) THEN
457 CALL cungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
458 $ childinfo )
459 lorgqrmin = max( lorgqrmin, p )
460 lorgqropt = max( lorgqropt, int( work(1) ) )
461 END IF
462 IF( wantu2 .AND. m-p .GT. 0 ) THEN
463 CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, cdum,
464 $ work(1), -1, childinfo )
465 lorgqrmin = max( lorgqrmin, m-p-1 )
466 lorgqropt = max( lorgqropt, int( work(1) ) )
467 END IF
468 IF( wantv1t .AND. q .GT. 0 ) THEN
469 CALL cunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
470 $ childinfo )
471 lorglqmin = max( lorglqmin, q )
472 lorglqopt = max( lorglqopt, int( work(1) ) )
473 END IF
474 CALL cbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
475 $ theta, dum, cdum, 1, v1t, ldv1t, u2, ldu2, u1,
476 $ ldu1, dum, dum, dum, dum, dum, dum, dum, dum,
477 $ rwork(1), -1, childinfo )
478 lbbcsd = int( rwork(1) )
479 ELSE
480 CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
481 $ dum,
482 $ cdum, cdum, cdum, cdum, work(1), -1, childinfo
483 $ )
484 lorbdb = m + int( work(1) )
485 IF( wantu1 .AND. p .GT. 0 ) THEN
486 CALL cungqr( p, p, m-q, u1, ldu1, cdum, work(1), -1,
487 $ childinfo )
488 lorgqrmin = max( lorgqrmin, p )
489 lorgqropt = max( lorgqropt, int( work(1) ) )
490 END IF
491 IF( wantu2 .AND. m-p .GT. 0 ) THEN
492 CALL cungqr( m-p, m-p, m-q, u2, ldu2, cdum, work(1),
493 $ -1,
494 $ childinfo )
495 lorgqrmin = max( lorgqrmin, m-p )
496 lorgqropt = max( lorgqropt, int( work(1) ) )
497 END IF
498 IF( wantv1t .AND. q .GT. 0 ) THEN
499 CALL cunglq( q, q, q, v1t, ldv1t, cdum, work(1), -1,
500 $ childinfo )
501 lorglqmin = max( lorglqmin, q )
502 lorglqopt = max( lorglqopt, int( work(1) ) )
503 END IF
504 CALL cbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
505 $ theta, dum, u2, ldu2, u1, ldu1, cdum, 1, v1t,
506 $ ldv1t, dum, dum, dum, dum, dum, dum, dum, dum,
507 $ rwork(1), -1, childinfo )
508 lbbcsd = int( rwork(1) )
509 END IF
510 lrworkmin = ibbcsd+lbbcsd-1
511 lrworkopt = lrworkmin
512 rwork(1) = real( lrworkopt )
513 lworkmin = max( iorbdb+lorbdb-1,
514 $ iorgqr+lorgqrmin-1,
515 $ iorglq+lorglqmin-1 )
516 lworkopt = max( iorbdb+lorbdb-1,
517 $ iorgqr+lorgqropt-1,
518 $ iorglq+lorglqopt-1 )
519 work(1) = sroundup_lwork(lworkopt)
520 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
521 info = -19
522 END IF
523 IF( lrwork .LT. lrworkmin .AND. .NOT.lquery ) THEN
524 info = -21
525 END IF
526 END IF
527 IF( info .NE. 0 ) THEN
528 CALL xerbla( 'CUNCSD2BY1', -info )
529 RETURN
530 ELSE IF( lquery ) THEN
531 RETURN
532 END IF
533 lorgqr = lwork-iorgqr+1
534 lorglq = lwork-iorglq+1
535*
536* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
537* in which R = MIN(P,M-P,Q,M-Q)
538*
539 IF( r .EQ. q ) THEN
540*
541* Case 1: R = Q
542*
543* Simultaneously bidiagonalize X11 and X21
544*
545 CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
546 $ rwork(iphi), work(itaup1), work(itaup2),
547 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
548*
549* Accumulate Householder reflectors
550*
551 IF( wantu1 .AND. p .GT. 0 ) THEN
552 CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
553 CALL cungqr( p, p, q, u1, ldu1, work(itaup1),
554 $ work(iorgqr),
555 $ lorgqr, childinfo )
556 END IF
557 IF( wantu2 .AND. m-p .GT. 0 ) THEN
558 CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
559 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
560 $ work(iorgqr), lorgqr, childinfo )
561 END IF
562 IF( wantv1t .AND. q .GT. 0 ) THEN
563 v1t(1,1) = one
564 DO j = 2, q
565 v1t(1,j) = zero
566 v1t(j,1) = zero
567 END DO
568 CALL clacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
569 $ ldv1t )
570 CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t,
571 $ work(itauq1),
572 $ work(iorglq), lorglq, childinfo )
573 END IF
574*
575* Simultaneously diagonalize X11 and X21.
576*
577 CALL cbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
578 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
579 $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
580 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
581 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
582 $ lrwork-ibbcsd+1, childinfo )
583*
584* Permute rows and columns to place zero submatrices in
585* preferred positions
586*
587 IF( q .GT. 0 .AND. wantu2 ) THEN
588 DO i = 1, q
589 iwork(i) = m - p - q + i
590 END DO
591 DO i = q + 1, m - p
592 iwork(i) = i - q
593 END DO
594 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
595 END IF
596 ELSE IF( r .EQ. p ) THEN
597*
598* Case 2: R = P
599*
600* Simultaneously bidiagonalize X11 and X21
601*
602 CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
603 $ rwork(iphi), work(itaup1), work(itaup2),
604 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
605*
606* Accumulate Householder reflectors
607*
608 IF( wantu1 .AND. p .GT. 0 ) THEN
609 u1(1,1) = one
610 DO j = 2, p
611 u1(1,j) = zero
612 u1(j,1) = zero
613 END DO
614 CALL clacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2),
615 $ ldu1 )
616 CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
617 $ work(iorgqr), lorgqr, childinfo )
618 END IF
619 IF( wantu2 .AND. m-p .GT. 0 ) THEN
620 CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
621 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
622 $ work(iorgqr), lorgqr, childinfo )
623 END IF
624 IF( wantv1t .AND. q .GT. 0 ) THEN
625 CALL clacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
626 CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
627 $ work(iorglq), lorglq, childinfo )
628 END IF
629*
630* Simultaneously diagonalize X11 and X21.
631*
632 CALL cbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
633 $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
634 $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
635 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
636 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
637 $ childinfo )
638*
639* Permute rows and columns to place identity submatrices in
640* preferred positions
641*
642 IF( q .GT. 0 .AND. wantu2 ) THEN
643 DO i = 1, q
644 iwork(i) = m - p - q + i
645 END DO
646 DO i = q + 1, m - p
647 iwork(i) = i - q
648 END DO
649 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
650 END IF
651 ELSE IF( r .EQ. m-p ) THEN
652*
653* Case 3: R = M-P
654*
655* Simultaneously bidiagonalize X11 and X21
656*
657 CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
658 $ rwork(iphi), work(itaup1), work(itaup2),
659 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
660*
661* Accumulate Householder reflectors
662*
663 IF( wantu1 .AND. p .GT. 0 ) THEN
664 CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
665 CALL cungqr( p, p, q, u1, ldu1, work(itaup1),
666 $ work(iorgqr),
667 $ lorgqr, childinfo )
668 END IF
669 IF( wantu2 .AND. m-p .GT. 0 ) THEN
670 u2(1,1) = one
671 DO j = 2, m-p
672 u2(1,j) = zero
673 u2(j,1) = zero
674 END DO
675 CALL clacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
676 $ ldu2 )
677 CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
678 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
679 END IF
680 IF( wantv1t .AND. q .GT. 0 ) THEN
681 CALL clacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
682 CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
683 $ work(iorglq), lorglq, childinfo )
684 END IF
685*
686* Simultaneously diagonalize X11 and X21.
687*
688 CALL cbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
689 $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
690 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
691 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
692 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
693 $ rwork(ibbcsd), lbbcsd, childinfo )
694*
695* Permute rows and columns to place identity submatrices in
696* preferred positions
697*
698 IF( q .GT. r ) THEN
699 DO i = 1, r
700 iwork(i) = q - r + i
701 END DO
702 DO i = r + 1, q
703 iwork(i) = i - r
704 END DO
705 IF( wantu1 ) THEN
706 CALL clapmt( .false., p, q, u1, ldu1, iwork )
707 END IF
708 IF( wantv1t ) THEN
709 CALL clapmr( .false., q, q, v1t, ldv1t, iwork )
710 END IF
711 END IF
712 ELSE
713*
714* Case 4: R = M-Q
715*
716* Simultaneously bidiagonalize X11 and X21
717*
718 CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
719 $ rwork(iphi), work(itaup1), work(itaup2),
720 $ work(itauq1), work(iorbdb), work(iorbdb+m),
721 $ lorbdb-m, childinfo )
722*
723* Accumulate Householder reflectors
724*
725
726 IF( wantu2 .AND. m-p .GT. 0 ) THEN
727 CALL ccopy( m-p, work(iorbdb+p), 1, u2, 1 )
728 END IF
729 IF( wantu1 .AND. p .GT. 0 ) THEN
730 CALL ccopy( p, work(iorbdb), 1, u1, 1 )
731 DO j = 2, p
732 u1(1,j) = zero
733 END DO
734 CALL clacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
735 $ ldu1 )
736 CALL cungqr( p, p, m-q, u1, ldu1, work(itaup1),
737 $ work(iorgqr), lorgqr, childinfo )
738 END IF
739 IF( wantu2 .AND. m-p .GT. 0 ) THEN
740 DO j = 2, m-p
741 u2(1,j) = zero
742 END DO
743 CALL clacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
744 $ ldu2 )
745 CALL cungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
746 $ work(iorgqr), lorgqr, childinfo )
747 END IF
748 IF( wantv1t .AND. q .GT. 0 ) THEN
749 CALL clacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
750 CALL clacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1),
751 $ ldx11,
752 $ v1t(m-q+1,m-q+1), ldv1t )
753 CALL clacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
754 $ v1t(p+1,p+1), ldv1t )
755 CALL cunglq( q, q, q, v1t, ldv1t, work(itauq1),
756 $ work(iorglq), lorglq, childinfo )
757 END IF
758*
759* Simultaneously diagonalize X11 and X21.
760*
761 CALL cbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
762 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
763 $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
764 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
765 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
766 $ rwork(ibbcsd), lbbcsd, childinfo )
767*
768* Permute rows and columns to place identity submatrices in
769* preferred positions
770*
771 IF( p .GT. r ) THEN
772 DO i = 1, r
773 iwork(i) = p - r + i
774 END DO
775 DO i = r + 1, p
776 iwork(i) = i - r
777 END DO
778 IF( wantu1 ) THEN
779 CALL clapmt( .false., p, p, u1, ldu1, iwork )
780 END IF
781 IF( wantv1t ) THEN
782 CALL clapmr( .false., p, q, v1t, ldv1t, iwork )
783 END IF
784 END IF
785 END IF
786*
787 RETURN
788*
789* End of CUNCSD2BY1
790*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cbbcsd
subroutine cbbcsd(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, rwork, lrwork, info)
CBBCSD
Definition cbbcsd.f:331
ccopy
subroutine ccopy(n, cx, incx, cy, incy)
CCOPY
Definition ccopy.f:81
clacpy
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
Definition clacpy.f:101
clapmr
subroutine clapmr(forwrd, m, n, x, ldx, k)
CLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition clapmr.f:102
clapmt
subroutine clapmt(forwrd, m, n, x, ldx, k)
CLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition clapmt.f:102
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sroundup_lwork
real function sroundup_lwork(lwork)
SROUNDUP_LWORK
Definition sroundup_lwork.f:59
cunbdb1
subroutine cunbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
CUNBDB1
Definition cunbdb1.f:201
cunbdb2
subroutine cunbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
CUNBDB2
Definition cunbdb2.f:201
cunbdb3
subroutine cunbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
CUNBDB3
Definition cunbdb3.f:201
cunbdb4
subroutine cunbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
CUNBDB4
Definition cunbdb4.f:212
cunglq
subroutine cunglq(m, n, k, a, lda, tau, work, lwork, info)
CUNGLQ
Definition cunglq.f:125
cungqr
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: