LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sorcsd2by1()

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:

                                [  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 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). 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 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.

Definition at line 230 of file sorcsd2by1.f.

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