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

◆ dorcsd2by1()

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

DORCSD2BY1

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

Purpose:
!>
!> DORCSD2BY1 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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:  DBBCSD 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 dorcsd2by1.f.

232*
233* -- LAPACK computational routine (3.5.0) --
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 DOUBLE PRECISION THETA(*)
244 DOUBLE PRECISION U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
245 $ X11(LDX11,*), X21(LDX21,*)
246 INTEGER IWORK(*)
247* ..
248*
249* =====================================================================
250*
251* .. Parameters ..
252 DOUBLE PRECISION ONE, ZERO
253 parameter( one = 1.0d0, zero = 0.0d0 )
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 DOUBLE PRECISION DUM1(1), DUM2(1,1)
266* ..
267* .. External Subroutines ..
268 EXTERNAL dbbcsd, dcopy, dlacpy, dlapmr, dlapmt,
269 $ dorbdb1,
270 $ dorbdb2, dorbdb3, dorbdb4, dorglq, dorgqr,
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* | DORBDB WORK | DORGQR WORK | DORGLQ WORK | B21D (R) |
323* | | | | B21E (R-1) |
324* | | | | B22D (R) |
325* | | | | B22E (R-1) |
326* | | | | DBBCSD 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 dorbdb1( m, p, q, x11, ldx11, x21, ldx21, theta,
352 $ dum1, dum1, dum1, dum1, work,
353 $ -1, childinfo )
354 lorbdb = int( work(1) )
355 IF( wantu1 .AND. p .GT. 0 ) THEN
356 CALL dorgqr( 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 dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
363 $ -1, 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 dorglq( 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 dbbcsd( jobu1, jobu2, jobv1t, 'N', 'N', m, p, q,
374 $ theta,
375 $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t,
376 $ dum2, 1, dum1, dum1, dum1,
377 $ dum1, dum1, dum1, dum1,
378 $ dum1, work(1), -1, childinfo )
379 lbbcsd = int( work(1) )
380 ELSE IF( r .EQ. p ) THEN
381 CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
382 $ dum1, dum1, dum1, dum1,
383 $ work(1), -1, childinfo )
384 lorbdb = int( work(1) )
385 IF( wantu1 .AND. p .GT. 0 ) THEN
386 CALL dorgqr( 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 dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
393 $ -1, 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 dorglq( 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 dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p,
404 $ theta,
405 $ dum1, v1t, ldv1t, dum2, 1, u1, ldu1,
406 $ u2, ldu2, dum1, dum1, dum1,
407 $ dum1, dum1, dum1, dum1,
408 $ dum1, work(1), -1, childinfo )
409 lbbcsd = int( work(1) )
410 ELSE IF( r .EQ. m-p ) THEN
411 CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
412 $ dum1, dum1, dum1, dum1,
413 $ work(1), -1, childinfo )
414 lorbdb = int( work(1) )
415 IF( wantu1 .AND. p .GT. 0 ) THEN
416 CALL dorgqr( 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 dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
423 $ dum1, 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 dorglq( 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 dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
434 $ theta, dum1, dum2, 1, v1t, ldv1t, u2,
435 $ ldu2, u1, ldu1, dum1, dum1, dum1,
436 $ dum1, dum1, dum1, dum1,
437 $ dum1, work(1), -1, childinfo )
438 lbbcsd = int( work(1) )
439 ELSE
440 CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
441 $ dum1, dum1, dum1, dum1,
442 $ dum1, work(1), -1, childinfo )
443 lorbdb = m + int( work(1) )
444 IF( wantu1 .AND. p .GT. 0 ) THEN
445 CALL dorgqr( 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 dorgqr( 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 dorglq( 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 dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
463 $ theta, dum1, u2, ldu2, u1, ldu1, dum2,
464 $ 1, v1t, ldv1t, dum1, dum1, dum1,
465 $ dum1, dum1, dum1, dum1,
466 $ dum1, work(1), -1, 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( 'DORCSD2BY1', -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 dorbdb1( 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 dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
508 CALL dorgqr( 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 dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
514 CALL dorgqr( 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 dlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
524 $ ldv1t )
525 CALL dorglq( 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 dbbcsd( 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),
535 $ work(ib12d), work(ib12e), work(ib21d),
536 $ work(ib21e), work(ib22d), work(ib22e),
537 $ work(ibbcsd), lbbcsd, 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 dlapmt( .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 dorbdb2( 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 dlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2),
570 $ ldu1 )
571 CALL dorgqr( 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 dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
576 CALL dorgqr( 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 dlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
581 CALL dorglq( 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 dbbcsd( 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 dlapmt( .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 dorbdb3( 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 dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
620 CALL dorgqr( 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 dlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
631 $ ldu2 )
632 CALL dorgqr( 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 dlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
637 CALL dorglq( 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 dbbcsd( '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 dlapmt( .false., p, q, u1, ldu1, iwork )
662 END IF
663 IF( wantv1t ) THEN
664 CALL dlapmr( .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 dorbdb4( 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 dcopy( m-p, work(iorbdb+p), 1, u2, 1 )
682 END IF
683 IF( wantu1 .AND. p .GT. 0 ) THEN
684 CALL dcopy( p, work(iorbdb), 1, u1, 1 )
685 DO j = 2, p
686 u1(1,j) = zero
687 END DO
688 CALL dlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
689 $ ldu1 )
690 CALL dorgqr( 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 dlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
698 $ ldu2 )
699 CALL dorgqr( 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 dlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
704 CALL dlacpy( '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 dlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
708 $ v1t(p+1,p+1), ldv1t )
709 CALL dorglq( 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 dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
716 $ theta, work(iphi), u2, ldu2, u1, ldu1, dum1,
717 $ 1, v1t, ldv1t, work(ib11d), work(ib11e),
718 $ work(ib12d), work(ib12e), work(ib21d),
719 $ work(ib21e), work(ib22d), work(ib22e),
720 $ work(ibbcsd), lbbcsd, 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 dlapmt( .false., p, p, u1, ldu1, iwork )
734 END IF
735 IF( wantv1t ) THEN
736 CALL dlapmr( .false., p, q, v1t, ldv1t, iwork )
737 END IF
738 END IF
739 END IF
740*
741 RETURN
742*
743* End of DORCSD2BY1
744*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dbbcsd
subroutine dbbcsd(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)
DBBCSD
Definition dbbcsd.f:331
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlapmr
subroutine dlapmr(forwrd, m, n, x, ldx, k)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition dlapmr.f:102
dlapmt
subroutine dlapmt(forwrd, m, n, x, ldx, k)
DLAPMT performs a forward or backward permutation of the columns of a matrix.
Definition dlapmt.f:102
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dorbdb1
subroutine dorbdb1(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB1
Definition dorbdb1.f:202
dorbdb2
subroutine dorbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB2
Definition dorbdb2.f:201
dorbdb3
subroutine dorbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB3
Definition dorbdb3.f:200
dorbdb4
subroutine dorbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
DORBDB4
Definition dorbdb4.f:212
dorglq
subroutine dorglq(m, n, k, a, lda, tau, work, lwork, info)
DORGLQ
Definition dorglq.f:125
dorgqr
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: