LAPACK 3.12.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 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, theta,
416 $ dum(1), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
417 $ 1, dum, dum, dum, dum, dum, dum, dum, dum,
418 $ rwork(1), -1, childinfo )
419 lbbcsd = int( rwork(1) )
420 ELSE IF( r .EQ. p ) THEN
421 CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
422 $ cdum, cdum, cdum, work(1), -1, childinfo )
423 lorbdb = int( work(1) )
424 IF( wantu1 .AND. p .GT. 0 ) THEN
425 CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, cdum, work(1),
426 $ -1, childinfo )
427 lorgqrmin = max( lorgqrmin, p-1 )
428 lorgqropt = max( lorgqropt, int( work(1) ) )
429 END IF
430 IF( wantu2 .AND. m-p .GT. 0 ) THEN
431 CALL cungqr( m-p, m-p, q, u2, ldu2, cdum, work(1), -1,
432 $ childinfo )
433 lorgqrmin = max( lorgqrmin, m-p )
434 lorgqropt = max( lorgqropt, int( work(1) ) )
435 END IF
436 IF( wantv1t .AND. q .GT. 0 ) THEN
437 CALL cunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
438 $ childinfo )
439 lorglqmin = max( lorglqmin, q )
440 lorglqopt = max( lorglqopt, int( work(1) ) )
441 END IF
442 CALL cbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
443 $ dum, v1t, ldv1t, cdum, 1, u1, ldu1, u2, ldu2,
444 $ dum, dum, dum, dum, dum, dum, dum, dum,
445 $ rwork(1), -1, childinfo )
446 lbbcsd = int( rwork(1) )
447 ELSE IF( r .EQ. m-p ) THEN
448 CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
449 $ cdum, cdum, cdum, work(1), -1, childinfo )
450 lorbdb = int( work(1) )
451 IF( wantu1 .AND. p .GT. 0 ) THEN
452 CALL cungqr( p, p, q, u1, ldu1, cdum, work(1), -1,
453 $ childinfo )
454 lorgqrmin = max( lorgqrmin, p )
455 lorgqropt = max( lorgqropt, int( work(1) ) )
456 END IF
457 IF( wantu2 .AND. m-p .GT. 0 ) THEN
458 CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, cdum,
459 $ work(1), -1, childinfo )
460 lorgqrmin = max( lorgqrmin, m-p-1 )
461 lorgqropt = max( lorgqropt, int( work(1) ) )
462 END IF
463 IF( wantv1t .AND. q .GT. 0 ) THEN
464 CALL cunglq( q, q, r, v1t, ldv1t, cdum, work(1), -1,
465 $ childinfo )
466 lorglqmin = max( lorglqmin, q )
467 lorglqopt = max( lorglqopt, int( work(1) ) )
468 END IF
469 CALL cbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
470 $ theta, dum, cdum, 1, v1t, ldv1t, u2, ldu2, u1,
471 $ ldu1, dum, dum, dum, dum, dum, dum, dum, dum,
472 $ rwork(1), -1, childinfo )
473 lbbcsd = int( rwork(1) )
474 ELSE
475 CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta, dum,
476 $ cdum, cdum, cdum, cdum, work(1), -1, childinfo
477 $ )
478 lorbdb = m + int( work(1) )
479 IF( wantu1 .AND. p .GT. 0 ) THEN
480 CALL cungqr( p, p, m-q, u1, ldu1, cdum, work(1), -1,
481 $ childinfo )
482 lorgqrmin = max( lorgqrmin, p )
483 lorgqropt = max( lorgqropt, int( work(1) ) )
484 END IF
485 IF( wantu2 .AND. m-p .GT. 0 ) THEN
486 CALL cungqr( m-p, m-p, m-q, u2, ldu2, cdum, work(1), -1,
487 $ childinfo )
488 lorgqrmin = max( lorgqrmin, m-p )
489 lorgqropt = max( lorgqropt, int( work(1) ) )
490 END IF
491 IF( wantv1t .AND. q .GT. 0 ) THEN
492 CALL cunglq( q, q, q, v1t, ldv1t, cdum, work(1), -1,
493 $ childinfo )
494 lorglqmin = max( lorglqmin, q )
495 lorglqopt = max( lorglqopt, int( work(1) ) )
496 END IF
497 CALL cbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
498 $ theta, dum, u2, ldu2, u1, ldu1, cdum, 1, v1t,
499 $ ldv1t, dum, dum, dum, dum, dum, dum, dum, dum,
500 $ rwork(1), -1, childinfo )
501 lbbcsd = int( rwork(1) )
502 END IF
503 lrworkmin = ibbcsd+lbbcsd-1
504 lrworkopt = lrworkmin
505 rwork(1) = lrworkopt
506 lworkmin = max( iorbdb+lorbdb-1,
507 $ iorgqr+lorgqrmin-1,
508 $ iorglq+lorglqmin-1 )
509 lworkopt = max( iorbdb+lorbdb-1,
510 $ iorgqr+lorgqropt-1,
511 $ iorglq+lorglqopt-1 )
512 work(1) = sroundup_lwork(lworkopt)
513 IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
514 info = -19
515 END IF
516 IF( lrwork .LT. lrworkmin .AND. .NOT.lquery ) THEN
517 info = -21
518 END IF
519 END IF
520 IF( info .NE. 0 ) THEN
521 CALL xerbla( 'CUNCSD2BY1', -info )
522 RETURN
523 ELSE IF( lquery ) THEN
524 RETURN
525 END IF
526 lorgqr = lwork-iorgqr+1
527 lorglq = lwork-iorglq+1
528*
529* Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
530* in which R = MIN(P,M-P,Q,M-Q)
531*
532 IF( r .EQ. q ) THEN
533*
534* Case 1: R = Q
535*
536* Simultaneously bidiagonalize X11 and X21
537*
538 CALL cunbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
539 $ rwork(iphi), work(itaup1), work(itaup2),
540 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
541*
542* Accumulate Householder reflectors
543*
544 IF( wantu1 .AND. p .GT. 0 ) THEN
545 CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
546 CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
547 $ lorgqr, childinfo )
548 END IF
549 IF( wantu2 .AND. m-p .GT. 0 ) THEN
550 CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
551 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
552 $ work(iorgqr), lorgqr, childinfo )
553 END IF
554 IF( wantv1t .AND. q .GT. 0 ) THEN
555 v1t(1,1) = one
556 DO j = 2, q
557 v1t(1,j) = zero
558 v1t(j,1) = zero
559 END DO
560 CALL clacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
561 $ ldv1t )
562 CALL cunglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
563 $ work(iorglq), lorglq, childinfo )
564 END IF
565*
566* Simultaneously diagonalize X11 and X21.
567*
568 CALL cbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
569 $ rwork(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t, cdum,
570 $ 1, rwork(ib11d), rwork(ib11e), rwork(ib12d),
571 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
572 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd),
573 $ lrwork-ibbcsd+1, childinfo )
574*
575* Permute rows and columns to place zero submatrices in
576* preferred positions
577*
578 IF( q .GT. 0 .AND. wantu2 ) THEN
579 DO i = 1, q
580 iwork(i) = m - p - q + i
581 END DO
582 DO i = q + 1, m - p
583 iwork(i) = i - q
584 END DO
585 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
586 END IF
587 ELSE IF( r .EQ. p ) THEN
588*
589* Case 2: R = P
590*
591* Simultaneously bidiagonalize X11 and X21
592*
593 CALL cunbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
594 $ rwork(iphi), work(itaup1), work(itaup2),
595 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
596*
597* Accumulate Householder reflectors
598*
599 IF( wantu1 .AND. p .GT. 0 ) THEN
600 u1(1,1) = one
601 DO j = 2, p
602 u1(1,j) = zero
603 u1(j,1) = zero
604 END DO
605 CALL clacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
606 CALL cungqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
607 $ work(iorgqr), lorgqr, childinfo )
608 END IF
609 IF( wantu2 .AND. m-p .GT. 0 ) THEN
610 CALL clacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
611 CALL cungqr( m-p, m-p, q, u2, ldu2, work(itaup2),
612 $ work(iorgqr), lorgqr, childinfo )
613 END IF
614 IF( wantv1t .AND. q .GT. 0 ) THEN
615 CALL clacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
616 CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
617 $ work(iorglq), lorglq, childinfo )
618 END IF
619*
620* Simultaneously diagonalize X11 and X21.
621*
622 CALL cbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
623 $ rwork(iphi), v1t, ldv1t, cdum, 1, u1, ldu1, u2,
624 $ ldu2, rwork(ib11d), rwork(ib11e), rwork(ib12d),
625 $ rwork(ib12e), rwork(ib21d), rwork(ib21e),
626 $ rwork(ib22d), rwork(ib22e), rwork(ibbcsd), lbbcsd,
627 $ childinfo )
628*
629* Permute rows and columns to place identity submatrices in
630* preferred positions
631*
632 IF( q .GT. 0 .AND. wantu2 ) THEN
633 DO i = 1, q
634 iwork(i) = m - p - q + i
635 END DO
636 DO i = q + 1, m - p
637 iwork(i) = i - q
638 END DO
639 CALL clapmt( .false., m-p, m-p, u2, ldu2, iwork )
640 END IF
641 ELSE IF( r .EQ. m-p ) THEN
642*
643* Case 3: R = M-P
644*
645* Simultaneously bidiagonalize X11 and X21
646*
647 CALL cunbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
648 $ rwork(iphi), work(itaup1), work(itaup2),
649 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
650*
651* Accumulate Householder reflectors
652*
653 IF( wantu1 .AND. p .GT. 0 ) THEN
654 CALL clacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
655 CALL cungqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
656 $ lorgqr, childinfo )
657 END IF
658 IF( wantu2 .AND. m-p .GT. 0 ) THEN
659 u2(1,1) = one
660 DO j = 2, m-p
661 u2(1,j) = zero
662 u2(j,1) = zero
663 END DO
664 CALL clacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
665 $ ldu2 )
666 CALL cungqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
667 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
668 END IF
669 IF( wantv1t .AND. q .GT. 0 ) THEN
670 CALL clacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
671 CALL cunglq( q, q, r, v1t, ldv1t, work(itauq1),
672 $ work(iorglq), lorglq, childinfo )
673 END IF
674*
675* Simultaneously diagonalize X11 and X21.
676*
677 CALL cbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
678 $ theta, rwork(iphi), cdum, 1, v1t, ldv1t, u2, ldu2,
679 $ u1, ldu1, rwork(ib11d), rwork(ib11e),
680 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
681 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
682 $ rwork(ibbcsd), lbbcsd, childinfo )
683*
684* Permute rows and columns to place identity submatrices in
685* preferred positions
686*
687 IF( q .GT. r ) THEN
688 DO i = 1, r
689 iwork(i) = q - r + i
690 END DO
691 DO i = r + 1, q
692 iwork(i) = i - r
693 END DO
694 IF( wantu1 ) THEN
695 CALL clapmt( .false., p, q, u1, ldu1, iwork )
696 END IF
697 IF( wantv1t ) THEN
698 CALL clapmr( .false., q, q, v1t, ldv1t, iwork )
699 END IF
700 END IF
701 ELSE
702*
703* Case 4: R = M-Q
704*
705* Simultaneously bidiagonalize X11 and X21
706*
707 CALL cunbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
708 $ rwork(iphi), work(itaup1), work(itaup2),
709 $ work(itauq1), work(iorbdb), work(iorbdb+m),
710 $ lorbdb-m, childinfo )
711*
712* Accumulate Householder reflectors
713*
714
715 IF( wantu2 .AND. m-p .GT. 0 ) THEN
716 CALL ccopy( m-p, work(iorbdb+p), 1, u2, 1 )
717 END IF
718 IF( wantu1 .AND. p .GT. 0 ) THEN
719 CALL ccopy( p, work(iorbdb), 1, u1, 1 )
720 DO j = 2, p
721 u1(1,j) = zero
722 END DO
723 CALL clacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
724 $ ldu1 )
725 CALL cungqr( p, p, m-q, u1, ldu1, work(itaup1),
726 $ work(iorgqr), lorgqr, childinfo )
727 END IF
728 IF( wantu2 .AND. m-p .GT. 0 ) THEN
729 DO j = 2, m-p
730 u2(1,j) = zero
731 END DO
732 CALL clacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
733 $ ldu2 )
734 CALL cungqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
735 $ work(iorgqr), lorgqr, childinfo )
736 END IF
737 IF( wantv1t .AND. q .GT. 0 ) THEN
738 CALL clacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
739 CALL clacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
740 $ v1t(m-q+1,m-q+1), ldv1t )
741 CALL clacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
742 $ v1t(p+1,p+1), ldv1t )
743 CALL cunglq( q, q, q, v1t, ldv1t, work(itauq1),
744 $ work(iorglq), lorglq, childinfo )
745 END IF
746*
747* Simultaneously diagonalize X11 and X21.
748*
749 CALL cbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
750 $ theta, rwork(iphi), u2, ldu2, u1, ldu1, cdum, 1,
751 $ v1t, ldv1t, rwork(ib11d), rwork(ib11e),
752 $ rwork(ib12d), rwork(ib12e), rwork(ib21d),
753 $ rwork(ib21e), rwork(ib22d), rwork(ib22e),
754 $ rwork(ibbcsd), lbbcsd, childinfo )
755*
756* Permute rows and columns to place identity submatrices in
757* preferred positions
758*
759 IF( p .GT. r ) THEN
760 DO i = 1, r
761 iwork(i) = p - r + i
762 END DO
763 DO i = r + 1, p
764 iwork(i) = i - r
765 END DO
766 IF( wantu1 ) THEN
767 CALL clapmt( .false., p, p, u1, ldu1, iwork )
768 END IF
769 IF( wantv1t ) THEN
770 CALL clapmr( .false., p, q, v1t, ldv1t, iwork )
771 END IF
772 END IF
773 END IF
774*
775 RETURN
776*
777* End of CUNCSD2BY1
778*
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:332
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:103
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
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
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:202
cunbdb2
subroutine cunbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
CUNBDB2
Definition cunbdb2.f:202
cunbdb3
subroutine cunbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
CUNBDB3
Definition cunbdb3.f:202
cunbdb4
subroutine cunbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
CUNBDB4
Definition cunbdb4.f:213
cunglq
subroutine cunglq(m, n, k, a, lda, tau, work, lwork, info)
CUNGLQ
Definition cunglq.f:127
cungqr
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:128
Here is the call graph for this function:
Here is the caller graph for this function: