LAPACK 3.11.0
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 253 of file cuncsd2by1.f.

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