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

◆ zuncsd2by1()

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

ZUNCSD2BY1

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

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