LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine sorcsd2by1 ( character  JOBU1,
character  JOBU2,
character  JOBV1T,
integer  M,
integer  P,
integer  Q,
real, dimension(ldx11,*)  X11,
integer  LDX11,
real, dimension(ldx21,*)  X21,
integer  LDX21,
real, dimension(*)  THETA,
real, dimension(ldu1,*)  U1,
integer  LDU1,
real, dimension(ldu2,*)  U2,
integer  LDU2,
real, dimension(ldv1t,*)  V1T,
integer  LDV1T,
real, dimension(*)  WORK,
integer  LWORK,
integer, dimension(*)  IWORK,
integer  INFO 
)

SORCSD2BY1

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

Purpose: 

 SORCSD2BY1 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:

                                [  I  0  0 ]
                                [  0  C  0 ]
          [ X11 ]   [ U1 |    ] [  0  0  0 ]
      X = [-----] = [---------] [----------] V1**T .
          [ X21 ]   [    | U2 ] [  0  0  0 ]
                                [  0  S  0 ]
                                [  0  0  I ]
 
 X11 is P-by-Q. The orthogonal 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).
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 REAL array, dimension (LDX11,Q)
          On entry, part of the orthogonal matrix whose CSD is desired.
[in]LDX11
          LDX11 is INTEGER
          The leading dimension of X11. LDX11 >= MAX(1,P).
[in,out]X21
          X21 is REAL array, dimension (LDX21,Q)
          On entry, part of the orthogonal 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 REAL array, dimension (P)
          If JOBU1 = 'Y', U1 contains the P-by-P orthogonal matrix U1.
[in]LDU1
          LDU1 is INTEGER
          The leading dimension of U1. If JOBU1 = 'Y', LDU1 >=
          MAX(1,P).
[out]U2
          U2 is REAL array, dimension (M-P)
          If JOBU2 = 'Y', U2 contains the (M-P)-by-(M-P) orthogonal
          matrix U2.
[in]LDU2
          LDU2 is INTEGER
          The leading dimension of U2. If JOBU2 = 'Y', LDU2 >=
          MAX(1,M-P).
[out]V1T
          V1T is REAL array, dimension (Q)
          If JOBV1T = 'Y', V1T contains the Q-by-Q matrix orthogonal
          matrix V1**T.
[in]LDV1T
          LDV1T is INTEGER
          The leading dimension of V1T. If JOBV1T = 'Y', LDV1T >=
          MAX(1,Q).
[out]WORK
          WORK is REAL array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
          If INFO > 0 on exit, WORK(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]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 array, returns
          this value as the first entry of the work array, and no error
          message related to LWORK 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:  SBBCSD 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.
Date
July 2012

Definition at line 234 of file sorcsd2by1.f.

234 *
235 * -- LAPACK computational routine (version 3.6.1) --
236 * -- LAPACK is a software package provided by Univ. of Tennessee, --
237 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
238 * July 2012
239 *
240 * .. Scalar Arguments ..
241  CHARACTER jobu1, jobu2, jobv1t
242  INTEGER info, ldu1, ldu2, ldv1t, lwork, ldx11, ldx21,
243  $ m, p, q
244 * ..
245 * .. Array Arguments ..
246  REAL theta(*)
247  REAL u1(ldu1,*), u2(ldu2,*), v1t(ldv1t,*), work(*),
248  $ x11(ldx11,*), x21(ldx21,*)
249  INTEGER iwork(*)
250 * ..
251 *
252 * =====================================================================
253 *
254 * .. Parameters ..
255  REAL one, zero
256  parameter ( one = 1.0e0, zero = 0.0e0 )
257 * ..
258 * .. Local Scalars ..
259  INTEGER childinfo, i, ib11d, ib11e, ib12d, ib12e,
260  $ ib21d, ib21e, ib22d, ib22e, ibbcsd, iorbdb,
261  $ iorglq, iorgqr, iphi, itaup1, itaup2, itauq1,
262  $ j, lbbcsd, lorbdb, lorglq, lorglqmin,
263  $ lorglqopt, lorgqr, lorgqrmin, lorgqropt,
264  $ lworkmin, lworkopt, r
265  LOGICAL lquery, wantu1, wantu2, wantv1t
266 * ..
267 * .. Local Arrays ..
268  REAL dum1(1), dum2(1,1)
269 * ..
270 * .. External Subroutines ..
271  EXTERNAL sbbcsd, scopy, slacpy, slapmr, slapmt, sorbdb1,
272  $ sorbdb2, sorbdb3, sorbdb4, sorglq, sorgqr,
273  $ xerbla
274 * ..
275 * .. External Functions ..
276  LOGICAL lsame
277  EXTERNAL lsame
278 * ..
279 * .. Intrinsic Function ..
280  INTRINSIC int, max, min
281 * ..
282 * .. Executable Statements ..
283 *
284 * Test input arguments
285 *
286  info = 0
287  wantu1 = lsame( jobu1, 'Y' )
288  wantu2 = lsame( jobu2, 'Y' )
289  wantv1t = lsame( jobv1t, 'Y' )
290  lquery = lwork .EQ. -1
291 *
292  IF( m .LT. 0 ) THEN
293  info = -4
294  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
295  info = -5
296  ELSE IF( q .LT. 0 .OR. q .GT. m ) THEN
297  info = -6
298  ELSE IF( ldx11 .LT. max( 1, p ) ) THEN
299  info = -8
300  ELSE IF( ldx21 .LT. max( 1, m-p ) ) THEN
301  info = -10
302  ELSE IF( wantu1 .AND. ldu1 .LT. max( 1, p ) ) THEN
303  info = -13
304  ELSE IF( wantu2 .AND. ldu2 .LT. max( 1, m - p ) ) THEN
305  info = -15
306  ELSE IF( wantv1t .AND. ldv1t .LT. max( 1, q ) ) THEN
307  info = -17
308  END IF
309 *
310  r = min( p, m-p, q, m-q )
311 *
312 * Compute workspace
313 *
314 * WORK layout:
315 * |-------------------------------------------------------|
316 * | LWORKOPT (1) |
317 * |-------------------------------------------------------|
318 * | PHI (MAX(1,R-1)) |
319 * |-------------------------------------------------------|
320 * | TAUP1 (MAX(1,P)) | B11D (R) |
321 * | TAUP2 (MAX(1,M-P)) | B11E (R-1) |
322 * | TAUQ1 (MAX(1,Q)) | B12D (R) |
323 * |-----------------------------------------| B12E (R-1) |
324 * | SORBDB WORK | SORGQR WORK | SORGLQ WORK | B21D (R) |
325 * | | | | B21E (R-1) |
326 * | | | | B22D (R) |
327 * | | | | B22E (R-1) |
328 * | | | | SBBCSD WORK |
329 * |-------------------------------------------------------|
330 *
331  IF( info .EQ. 0 ) THEN
332  iphi = 2
333  ib11d = iphi + max( 1, r-1 )
334  ib11e = ib11d + max( 1, r )
335  ib12d = ib11e + max( 1, r - 1 )
336  ib12e = ib12d + max( 1, r )
337  ib21d = ib12e + max( 1, r - 1 )
338  ib21e = ib21d + max( 1, r )
339  ib22d = ib21e + max( 1, r - 1 )
340  ib22e = ib22d + max( 1, r )
341  ibbcsd = ib22e + max( 1, r - 1 )
342  itaup1 = iphi + max( 1, r-1 )
343  itaup2 = itaup1 + max( 1, p )
344  itauq1 = itaup2 + max( 1, m-p )
345  iorbdb = itauq1 + max( 1, q )
346  iorgqr = itauq1 + max( 1, q )
347  iorglq = itauq1 + max( 1, q )
348  lorgqrmin = 1
349  lorgqropt = 1
350  lorglqmin = 1
351  lorglqopt = 1
352  IF( r .EQ. q ) THEN
353  CALL sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
354  $ dum1, dum1, dum1, dum1, work, -1,
355  $ childinfo )
356  lorbdb = int( work(1) )
357  IF( wantu1 .AND. p .GT. 0 ) THEN
358  CALL sorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
359  $ childinfo )
360  lorgqrmin = max( lorgqrmin, p )
361  lorgqropt = max( lorgqropt, int( work(1) ) )
362  ENDIF
363  IF( wantu2 .AND. m-p .GT. 0 ) THEN
364  CALL sorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1), -1,
365  $ childinfo )
366  lorgqrmin = max( lorgqrmin, m-p )
367  lorgqropt = max( lorgqropt, int( work(1) ) )
368  END IF
369  IF( wantv1t .AND. q .GT. 0 ) THEN
370  CALL sorglq( q-1, q-1, q-1, v1t, ldv1t,
371  $ dum1, work(1), -1, childinfo )
372  lorglqmin = max( lorglqmin, q-1 )
373  lorglqopt = max( lorglqopt, int( work(1) ) )
374  END IF
375  CALL sbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
376  $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t, dum2,
377  $ 1, dum1, dum1, dum1, dum1, dum1,
378  $ dum1, dum1, dum1, work(1), -1, childinfo
379  $ )
380  lbbcsd = int( work(1) )
381  ELSE IF( r .EQ. p ) THEN
382  CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
383  $ dum1, dum1, dum1, dum1, work(1), -1,
384  $ childinfo )
385  lorbdb = int( work(1) )
386  IF( wantu1 .AND. p .GT. 0 ) THEN
387  CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, dum1,
388  $ work(1), -1, childinfo )
389  lorgqrmin = max( lorgqrmin, p-1 )
390  lorgqropt = max( lorgqropt, int( work(1) ) )
391  END IF
392  IF( wantu2 .AND. m-p .GT. 0 ) THEN
393  CALL sorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1), -1,
394  $ childinfo )
395  lorgqrmin = max( lorgqrmin, m-p )
396  lorgqropt = max( lorgqropt, int( work(1) ) )
397  END IF
398  IF( wantv1t .AND. q .GT. 0 ) THEN
399  CALL sorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
400  $ childinfo )
401  lorglqmin = max( lorglqmin, q )
402  lorglqopt = max( lorglqopt, int( work(1) ) )
403  END IF
404  CALL sbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
405  $ dum1, v1t, ldv1t, dum2, 1, u1, ldu1, u2,
406  $ ldu2, dum1, dum1, dum1, dum1, dum1,
407  $ dum1, dum1, dum1, work(1), -1, childinfo
408  $ )
409  lbbcsd = int( work(1) )
410  ELSE IF( r .EQ. m-p ) THEN
411  CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
412  $ dum1, dum1, dum1, dum1, work(1), -1,
413  $ childinfo )
414  lorbdb = int( work(1) )
415  IF( wantu1 .AND. p .GT. 0 ) THEN
416  CALL sorgqr( p, p, q, u1, ldu1, dum1, work(1), -1,
417  $ childinfo )
418  lorgqrmin = max( lorgqrmin, p )
419  lorgqropt = max( lorgqropt, int( work(1) ) )
420  END IF
421  IF( wantu2 .AND. m-p .GT. 0 ) THEN
422  CALL sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2, dum1,
423  $ work(1), -1, childinfo )
424  lorgqrmin = max( lorgqrmin, m-p-1 )
425  lorgqropt = max( lorgqropt, int( work(1) ) )
426  END IF
427  IF( wantv1t .AND. q .GT. 0 ) THEN
428  CALL sorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
429  $ childinfo )
430  lorglqmin = max( lorglqmin, q )
431  lorglqopt = max( lorglqopt, int( work(1) ) )
432  END IF
433  CALL sbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
434  $ theta, dum1, dum2, 1, v1t, ldv1t, u2, ldu2,
435  $ u1, ldu1, dum1, dum1, dum1, dum1,
436  $ dum1, dum1, dum1, dum1, work(1), -1,
437  $ childinfo )
438  lbbcsd = int( work(1) )
439  ELSE
440  CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
441  $ dum1, dum1, dum1, dum1, dum1,
442  $ work(1), -1, childinfo )
443  lorbdb = m + int( work(1) )
444  IF( wantu1 .AND. p .GT. 0 ) THEN
445  CALL sorgqr( p, p, m-q, u1, ldu1, dum1, work(1), -1,
446  $ childinfo )
447  lorgqrmin = max( lorgqrmin, p )
448  lorgqropt = max( lorgqropt, int( work(1) ) )
449  END IF
450  IF( wantu2 .AND. m-p .GT. 0 ) THEN
451  CALL sorgqr( m-p, m-p, m-q, u2, ldu2, dum1, work(1),
452  $ -1, childinfo )
453  lorgqrmin = max( lorgqrmin, m-p )
454  lorgqropt = max( lorgqropt, int( work(1) ) )
455  END IF
456  IF( wantv1t .AND. q .GT. 0 ) THEN
457  CALL sorglq( q, q, q, v1t, ldv1t, dum1, work(1), -1,
458  $ childinfo )
459  lorglqmin = max( lorglqmin, q )
460  lorglqopt = max( lorglqopt, int( work(1) ) )
461  END IF
462  CALL sbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
463  $ theta, dum1, u2, ldu2, u1, ldu1, dum2, 1,
464  $ v1t, ldv1t, dum1, dum1, dum1, dum1,
465  $ dum1, dum1, dum1, dum1, work(1), -1,
466  $ childinfo )
467  lbbcsd = int( work(1) )
468  END IF
469  lworkmin = max( iorbdb+lorbdb-1,
470  $ iorgqr+lorgqrmin-1,
471  $ iorglq+lorglqmin-1,
472  $ ibbcsd+lbbcsd-1 )
473  lworkopt = max( iorbdb+lorbdb-1,
474  $ iorgqr+lorgqropt-1,
475  $ iorglq+lorglqopt-1,
476  $ ibbcsd+lbbcsd-1 )
477  work(1) = lworkopt
478  IF( lwork .LT. lworkmin .AND. .NOT.lquery ) THEN
479  info = -19
480  END IF
481  END IF
482  IF( info .NE. 0 ) THEN
483  CALL xerbla( 'SORCSD2BY1', -info )
484  RETURN
485  ELSE IF( lquery ) THEN
486  RETURN
487  END IF
488  lorgqr = lwork-iorgqr+1
489  lorglq = lwork-iorglq+1
490 *
491 * Handle four cases separately: R = Q, R = P, R = M-P, and R = M-Q,
492 * in which R = MIN(P,M-P,Q,M-Q)
493 *
494  IF( r .EQ. q ) THEN
495 *
496 * Case 1: R = Q
497 *
498 * Simultaneously bidiagonalize X11 and X21
499 *
500  CALL sorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
501  $ work(iphi), work(itaup1), work(itaup2),
502  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
503 *
504 * Accumulate Householder reflectors
505 *
506  IF( wantu1 .AND. p .GT. 0 ) THEN
507  CALL slacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
508  CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
509  $ lorgqr, childinfo )
510  END IF
511  IF( wantu2 .AND. m-p .GT. 0 ) THEN
512  CALL slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
513  CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
514  $ work(iorgqr), lorgqr, childinfo )
515  END IF
516  IF( wantv1t .AND. q .GT. 0 ) THEN
517  v1t(1,1) = one
518  DO j = 2, q
519  v1t(1,j) = zero
520  v1t(j,1) = zero
521  END DO
522  CALL slacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
523  $ ldv1t )
524  CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t, work(itauq1),
525  $ work(iorglq), lorglq, childinfo )
526  END IF
527 *
528 * Simultaneously diagonalize X11 and X21.
529 *
530  CALL sbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
531  $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t,
532  $ dum2, 1, work(ib11d), work(ib11e), work(ib12d),
533  $ work(ib12e), work(ib21d), work(ib21e),
534  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
535  $ childinfo )
536 *
537 * Permute rows and columns to place zero submatrices in
538 * preferred positions
539 *
540  IF( q .GT. 0 .AND. wantu2 ) THEN
541  DO i = 1, q
542  iwork(i) = m - p - q + i
543  END DO
544  DO i = q + 1, m - p
545  iwork(i) = i - q
546  END DO
547  CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
548  END IF
549  ELSE IF( r .EQ. p ) THEN
550 *
551 * Case 2: R = P
552 *
553 * Simultaneously bidiagonalize X11 and X21
554 *
555  CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
556  $ work(iphi), work(itaup1), work(itaup2),
557  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
558 *
559 * Accumulate Householder reflectors
560 *
561  IF( wantu1 .AND. p .GT. 0 ) THEN
562  u1(1,1) = one
563  DO j = 2, p
564  u1(1,j) = zero
565  u1(j,1) = zero
566  END DO
567  CALL slacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
568  CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
569  $ work(iorgqr), lorgqr, childinfo )
570  END IF
571  IF( wantu2 .AND. m-p .GT. 0 ) THEN
572  CALL slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
573  CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
574  $ work(iorgqr), lorgqr, childinfo )
575  END IF
576  IF( wantv1t .AND. q .GT. 0 ) THEN
577  CALL slacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
578  CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
579  $ work(iorglq), lorglq, childinfo )
580  END IF
581 *
582 * Simultaneously diagonalize X11 and X21.
583 *
584  CALL sbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
585  $ work(iphi), v1t, ldv1t, dum1, 1, u1, ldu1, u2,
586  $ ldu2, work(ib11d), work(ib11e), work(ib12d),
587  $ work(ib12e), work(ib21d), work(ib21e),
588  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
589  $ childinfo )
590 *
591 * Permute rows and columns to place identity submatrices in
592 * preferred positions
593 *
594  IF( q .GT. 0 .AND. wantu2 ) THEN
595  DO i = 1, q
596  iwork(i) = m - p - q + i
597  END DO
598  DO i = q + 1, m - p
599  iwork(i) = i - q
600  END DO
601  CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
602  END IF
603  ELSE IF( r .EQ. m-p ) THEN
604 *
605 * Case 3: R = M-P
606 *
607 * Simultaneously bidiagonalize X11 and X21
608 *
609  CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
610  $ work(iphi), work(itaup1), work(itaup2),
611  $ work(itauq1), work(iorbdb), lorbdb, childinfo )
612 *
613 * Accumulate Householder reflectors
614 *
615  IF( wantu1 .AND. p .GT. 0 ) THEN
616  CALL slacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
617  CALL sorgqr( p, p, q, u1, ldu1, work(itaup1), work(iorgqr),
618  $ lorgqr, childinfo )
619  END IF
620  IF( wantu2 .AND. m-p .GT. 0 ) THEN
621  u2(1,1) = one
622  DO j = 2, m-p
623  u2(1,j) = zero
624  u2(j,1) = zero
625  END DO
626  CALL slacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
627  $ ldu2 )
628  CALL sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
629  $ work(itaup2), work(iorgqr), lorgqr, childinfo )
630  END IF
631  IF( wantv1t .AND. q .GT. 0 ) THEN
632  CALL slacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
633  CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
634  $ work(iorglq), lorglq, childinfo )
635  END IF
636 *
637 * Simultaneously diagonalize X11 and X21.
638 *
639  CALL sbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
640  $ theta, work(iphi), dum1, 1, v1t, ldv1t, u2,
641  $ ldu2, u1, ldu1, work(ib11d), work(ib11e),
642  $ work(ib12d), work(ib12e), work(ib21d),
643  $ work(ib21e), work(ib22d), work(ib22e),
644  $ work(ibbcsd), lbbcsd, childinfo )
645 *
646 * Permute rows and columns to place identity submatrices in
647 * preferred positions
648 *
649  IF( q .GT. r ) THEN
650  DO i = 1, r
651  iwork(i) = q - r + i
652  END DO
653  DO i = r + 1, q
654  iwork(i) = i - r
655  END DO
656  IF( wantu1 ) THEN
657  CALL slapmt( .false., p, q, u1, ldu1, iwork )
658  END IF
659  IF( wantv1t ) THEN
660  CALL slapmr( .false., q, q, v1t, ldv1t, iwork )
661  END IF
662  END IF
663  ELSE
664 *
665 * Case 4: R = M-Q
666 *
667 * Simultaneously bidiagonalize X11 and X21
668 *
669  CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
670  $ work(iphi), work(itaup1), work(itaup2),
671  $ work(itauq1), work(iorbdb), work(iorbdb+m),
672  $ lorbdb-m, childinfo )
673 *
674 * Accumulate Householder reflectors
675 *
676  IF( wantu1 .AND. p .GT. 0 ) THEN
677  CALL scopy( p, work(iorbdb), 1, u1, 1 )
678  DO j = 2, p
679  u1(1,j) = zero
680  END DO
681  CALL slacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
682  $ ldu1 )
683  CALL sorgqr( p, p, m-q, u1, ldu1, work(itaup1),
684  $ work(iorgqr), lorgqr, childinfo )
685  END IF
686  IF( wantu2 .AND. m-p .GT. 0 ) THEN
687  CALL scopy( m-p, work(iorbdb+p), 1, u2, 1 )
688  DO j = 2, m-p
689  u2(1,j) = zero
690  END DO
691  CALL slacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
692  $ ldu2 )
693  CALL sorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
694  $ work(iorgqr), lorgqr, childinfo )
695  END IF
696  IF( wantv1t .AND. q .GT. 0 ) THEN
697  CALL slacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
698  CALL slacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1), ldx11,
699  $ v1t(m-q+1,m-q+1), ldv1t )
700  CALL slacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
701  $ v1t(p+1,p+1), ldv1t )
702  CALL sorglq( q, q, q, v1t, ldv1t, work(itauq1),
703  $ work(iorglq), lorglq, childinfo )
704  END IF
705 *
706 * Simultaneously diagonalize X11 and X21.
707 *
708  CALL sbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
709  $ theta, work(iphi), u2, ldu2, u1, ldu1, dum1, 1,
710  $ v1t, ldv1t, work(ib11d), work(ib11e), work(ib12d),
711  $ work(ib12e), work(ib21d), work(ib21e),
712  $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
713  $ childinfo )
714 *
715 * Permute rows and columns to place identity submatrices in
716 * preferred positions
717 *
718  IF( p .GT. r ) THEN
719  DO i = 1, r
720  iwork(i) = p - r + i
721  END DO
722  DO i = r + 1, p
723  iwork(i) = i - r
724  END DO
725  IF( wantu1 ) THEN
726  CALL slapmt( .false., p, p, u1, ldu1, iwork )
727  END IF
728  IF( wantv1t ) THEN
729  CALL slapmr( .false., p, q, v1t, ldv1t, iwork )
730  END IF
731  END IF
732  END IF
733 *
734  RETURN
735 *
736 * End of SORCSD2BY1
737 *
sorglq
subroutine sorglq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGLQ
Definition: sorglq.f:129
slapmt
subroutine slapmt(FORWRD, M, N, X, LDX, K)
SLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition: slapmt.f:106
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
slacpy
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
sorbdb4
subroutine sorbdb4(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, PHANTOM, WORK, LWORK, INFO)
SORBDB4
Definition: sorbdb4.f:216
slapmr
subroutine slapmr(FORWRD, M, N, X, LDX, K)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition: slapmr.f:106
sorbdb2
subroutine sorbdb2(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB2
Definition: sorbdb2.f:203
sorbdb1
subroutine sorbdb1(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB1
Definition: sorbdb1.f:205
sorbdb3
subroutine sorbdb3(M, P, Q, X11, LDX11, X21, LDX21, THETA, PHI, TAUP1, TAUP2, TAUQ1, WORK, LWORK, INFO)
SORBDB3
Definition: sorbdb3.f:204
sbbcsd
subroutine sbbcsd(JOBU1, JOBU2, JOBV1T, JOBV2T, TRANS, M, P, Q, THETA, PHI, U1, LDU1, U2, LDU2, V1T, LDV1T, V2T, LDV2T, B11D, B11E, B12D, B12E, B21D, B21E, B22D, B22E, WORK, LWORK, INFO)
SBBCSD
Definition: sbbcsd.f:334
sorgqr
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:130
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: