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

◆ 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 228 of file sorcsd2by1.f.

232*
233* -- LAPACK computational routine --
234* -- LAPACK is a software package provided by Univ. of Tennessee, --
235* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
236*
237* .. Scalar Arguments ..
238 CHARACTER JOBU1, JOBU2, JOBV1T
239 INTEGER INFO, LDU1, LDU2, LDV1T, LWORK, LDX11, LDX21,
240 $ M, P, Q
241* ..
242* .. Array Arguments ..
243 REAL THETA(*)
244 REAL U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
245 $ X11(LDX11,*), X21(LDX21,*)
246 INTEGER IWORK(*)
247* ..
248*
249* =====================================================================
250*
251* .. Parameters ..
252 REAL ONE, ZERO
253 parameter( one = 1.0e0, zero = 0.0e0 )
254* ..
255* .. Local Scalars ..
256 INTEGER CHILDINFO, I, IB11D, IB11E, IB12D, IB12E,
257 $ IB21D, IB21E, IB22D, IB22E, IBBCSD, IORBDB,
258 $ IORGLQ, IORGQR, IPHI, ITAUP1, ITAUP2, ITAUQ1,
259 $ J, LBBCSD, LORBDB, LORGLQ, LORGLQMIN,
260 $ LORGLQOPT, LORGQR, LORGQRMIN, LORGQROPT,
261 $ LWORKMIN, LWORKOPT, R
262 LOGICAL LQUERY, WANTU1, WANTU2, WANTV1T
263* ..
264* .. Local Arrays ..
265 REAL DUM1(1), DUM2(1,1)
266* ..
267* .. External Subroutines ..
268 EXTERNAL sbbcsd, scopy, slacpy, slapmr, slapmt,
269 $ sorbdb1,
270 $ sorbdb2, sorbdb3, sorbdb4, sorglq, sorgqr,
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,
374 $ theta,
375 $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t, dum2,
376 $ 1, dum1, dum1, dum1, dum1, dum1,
377 $ dum1, dum1, dum1, work(1), -1, childinfo
378 $ )
379 lbbcsd = int( work(1) )
380 ELSE IF( r .EQ. p ) THEN
381 CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
382 $ dum1, dum1, dum1, dum1, work(1), -1,
383 $ childinfo )
384 lorbdb = int( work(1) )
385 IF( wantu1 .AND. p .GT. 0 ) THEN
386 CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, dum1,
387 $ work(1), -1, childinfo )
388 lorgqrmin = max( lorgqrmin, p-1 )
389 lorgqropt = max( lorgqropt, int( work(1) ) )
390 END IF
391 IF( wantu2 .AND. m-p .GT. 0 ) THEN
392 CALL sorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1), -1,
393 $ childinfo )
394 lorgqrmin = max( lorgqrmin, m-p )
395 lorgqropt = max( lorgqropt, int( work(1) ) )
396 END IF
397 IF( wantv1t .AND. q .GT. 0 ) THEN
398 CALL sorglq( q, q, r, v1t, ldv1t, dum1, work(1), -1,
399 $ childinfo )
400 lorglqmin = max( lorglqmin, q )
401 lorglqopt = max( lorglqopt, int( work(1) ) )
402 END IF
403 CALL sbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p,
404 $ 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) = real( 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),
509 $ work(iorgqr),
510 $ lorgqr, childinfo )
511 END IF
512 IF( wantu2 .AND. m-p .GT. 0 ) THEN
513 CALL slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
514 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
515 $ work(iorgqr), lorgqr, childinfo )
516 END IF
517 IF( wantv1t .AND. q .GT. 0 ) THEN
518 v1t(1,1) = one
519 DO j = 2, q
520 v1t(1,j) = zero
521 v1t(j,1) = zero
522 END DO
523 CALL slacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
524 $ ldv1t )
525 CALL sorglq( q-1, q-1, q-1, v1t(2,2), ldv1t,
526 $ work(itauq1),
527 $ work(iorglq), lorglq, childinfo )
528 END IF
529*
530* Simultaneously diagonalize X11 and X21.
531*
532 CALL sbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q, theta,
533 $ work(iphi), u1, ldu1, u2, ldu2, v1t, ldv1t,
534 $ dum2, 1, work(ib11d), work(ib11e), work(ib12d),
535 $ work(ib12e), work(ib21d), work(ib21e),
536 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
537 $ childinfo )
538*
539* Permute rows and columns to place zero submatrices in
540* preferred positions
541*
542 IF( q .GT. 0 .AND. wantu2 ) THEN
543 DO i = 1, q
544 iwork(i) = m - p - q + i
545 END DO
546 DO i = q + 1, m - p
547 iwork(i) = i - q
548 END DO
549 CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
550 END IF
551 ELSE IF( r .EQ. p ) THEN
552*
553* Case 2: R = P
554*
555* Simultaneously bidiagonalize X11 and X21
556*
557 CALL sorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
558 $ work(iphi), work(itaup1), work(itaup2),
559 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
560*
561* Accumulate Householder reflectors
562*
563 IF( wantu1 .AND. p .GT. 0 ) THEN
564 u1(1,1) = one
565 DO j = 2, p
566 u1(1,j) = zero
567 u1(j,1) = zero
568 END DO
569 CALL slacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2),
570 $ ldu1 )
571 CALL sorgqr( p-1, p-1, p-1, u1(2,2), ldu1, work(itaup1),
572 $ work(iorgqr), lorgqr, childinfo )
573 END IF
574 IF( wantu2 .AND. m-p .GT. 0 ) THEN
575 CALL slacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
576 CALL sorgqr( m-p, m-p, q, u2, ldu2, work(itaup2),
577 $ work(iorgqr), lorgqr, childinfo )
578 END IF
579 IF( wantv1t .AND. q .GT. 0 ) THEN
580 CALL slacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
581 CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
582 $ work(iorglq), lorglq, childinfo )
583 END IF
584*
585* Simultaneously diagonalize X11 and X21.
586*
587 CALL sbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
588 $ work(iphi), v1t, ldv1t, dum1, 1, u1, ldu1, u2,
589 $ ldu2, work(ib11d), work(ib11e), work(ib12d),
590 $ work(ib12e), work(ib21d), work(ib21e),
591 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
592 $ childinfo )
593*
594* Permute rows and columns to place identity submatrices in
595* preferred positions
596*
597 IF( q .GT. 0 .AND. wantu2 ) THEN
598 DO i = 1, q
599 iwork(i) = m - p - q + i
600 END DO
601 DO i = q + 1, m - p
602 iwork(i) = i - q
603 END DO
604 CALL slapmt( .false., m-p, m-p, u2, ldu2, iwork )
605 END IF
606 ELSE IF( r .EQ. m-p ) THEN
607*
608* Case 3: R = M-P
609*
610* Simultaneously bidiagonalize X11 and X21
611*
612 CALL sorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
613 $ work(iphi), work(itaup1), work(itaup2),
614 $ work(itauq1), work(iorbdb), lorbdb, childinfo )
615*
616* Accumulate Householder reflectors
617*
618 IF( wantu1 .AND. p .GT. 0 ) THEN
619 CALL slacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
620 CALL sorgqr( p, p, q, u1, ldu1, work(itaup1),
621 $ work(iorgqr),
622 $ lorgqr, childinfo )
623 END IF
624 IF( wantu2 .AND. m-p .GT. 0 ) THEN
625 u2(1,1) = one
626 DO j = 2, m-p
627 u2(1,j) = zero
628 u2(j,1) = zero
629 END DO
630 CALL slacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
631 $ ldu2 )
632 CALL sorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
633 $ work(itaup2), work(iorgqr), lorgqr, childinfo )
634 END IF
635 IF( wantv1t .AND. q .GT. 0 ) THEN
636 CALL slacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
637 CALL sorglq( q, q, r, v1t, ldv1t, work(itauq1),
638 $ work(iorglq), lorglq, childinfo )
639 END IF
640*
641* Simultaneously diagonalize X11 and X21.
642*
643 CALL sbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
644 $ theta, work(iphi), dum1, 1, v1t, ldv1t, u2,
645 $ ldu2, u1, ldu1, work(ib11d), work(ib11e),
646 $ work(ib12d), work(ib12e), work(ib21d),
647 $ work(ib21e), work(ib22d), work(ib22e),
648 $ work(ibbcsd), lbbcsd, childinfo )
649*
650* Permute rows and columns to place identity submatrices in
651* preferred positions
652*
653 IF( q .GT. r ) THEN
654 DO i = 1, r
655 iwork(i) = q - r + i
656 END DO
657 DO i = r + 1, q
658 iwork(i) = i - r
659 END DO
660 IF( wantu1 ) THEN
661 CALL slapmt( .false., p, q, u1, ldu1, iwork )
662 END IF
663 IF( wantv1t ) THEN
664 CALL slapmr( .false., q, q, v1t, ldv1t, iwork )
665 END IF
666 END IF
667 ELSE
668*
669* Case 4: R = M-Q
670*
671* Simultaneously bidiagonalize X11 and X21
672*
673 CALL sorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
674 $ work(iphi), work(itaup1), work(itaup2),
675 $ work(itauq1), work(iorbdb), work(iorbdb+m),
676 $ lorbdb-m, childinfo )
677*
678* Accumulate Householder reflectors
679*
680 IF( wantu2 .AND. m-p .GT. 0 ) THEN
681 CALL scopy( m-p, work(iorbdb+p), 1, u2, 1 )
682 END IF
683 IF( wantu1 .AND. p .GT. 0 ) THEN
684 CALL scopy( p, work(iorbdb), 1, u1, 1 )
685 DO j = 2, p
686 u1(1,j) = zero
687 END DO
688 CALL slacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
689 $ ldu1 )
690 CALL sorgqr( p, p, m-q, u1, ldu1, work(itaup1),
691 $ work(iorgqr), lorgqr, childinfo )
692 END IF
693 IF( wantu2 .AND. m-p .GT. 0 ) THEN
694 DO j = 2, m-p
695 u2(1,j) = zero
696 END DO
697 CALL slacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
698 $ ldu2 )
699 CALL sorgqr( m-p, m-p, m-q, u2, ldu2, work(itaup2),
700 $ work(iorgqr), lorgqr, childinfo )
701 END IF
702 IF( wantv1t .AND. q .GT. 0 ) THEN
703 CALL slacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
704 CALL slacpy( 'U', p-(m-q), q-(m-q), x11(m-q+1,m-q+1),
705 $ ldx11,
706 $ v1t(m-q+1,m-q+1), ldv1t )
707 CALL slacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
708 $ v1t(p+1,p+1), ldv1t )
709 CALL sorglq( q, q, q, v1t, ldv1t, work(itauq1),
710 $ work(iorglq), lorglq, childinfo )
711 END IF
712*
713* Simultaneously diagonalize X11 and X21.
714*
715 CALL sbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
716 $ theta, work(iphi), u2, ldu2, u1, ldu1, dum1, 1,
717 $ v1t, ldv1t, work(ib11d), work(ib11e), work(ib12d),
718 $ work(ib12e), work(ib21d), work(ib21e),
719 $ work(ib22d), work(ib22e), work(ibbcsd), lbbcsd,
720 $ childinfo )
721*
722* Permute rows and columns to place identity submatrices in
723* preferred positions
724*
725 IF( p .GT. r ) THEN
726 DO i = 1, r
727 iwork(i) = p - r + i
728 END DO
729 DO i = r + 1, p
730 iwork(i) = i - r
731 END DO
732 IF( wantu1 ) THEN
733 CALL slapmt( .false., p, p, u1, ldu1, iwork )
734 END IF
735 IF( wantv1t ) THEN
736 CALL slapmr( .false., p, q, v1t, ldv1t, iwork )
737 END IF
738 END IF
739 END IF
740*
741 RETURN
742*
743* End of SORCSD2BY1
744*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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:331
scopy
subroutine scopy(n, sx, incx, sy, incy)
SCOPY
Definition scopy.f:82
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:101
slapmr
subroutine slapmr(forwrd, m, n, x, ldx, k)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition slapmr.f:102
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:102
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
sorbdb1
subroutine sorbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
SORBDB1
Definition sorbdb1.f:202
sorbdb2
subroutine sorbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
SORBDB2
Definition sorbdb2.f:200
sorbdb3
subroutine sorbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
SORBDB3
Definition sorbdb3.f:201
sorbdb4
subroutine sorbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
SORBDB4
Definition sorbdb4.f:213
sorglq
subroutine sorglq(m, n, k, a, lda, tau, work, lwork, info)
SORGLQ
Definition sorglq.f:125
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: