LAPACK 3.12.0
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 230 of file dorcsd2by1.f.

233*
234* -- LAPACK computational routine (3.5.0) --
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 DOUBLE PRECISION THETA(*)
245 DOUBLE PRECISION U1(LDU1,*), U2(LDU2,*), V1T(LDV1T,*), WORK(*),
246 $ X11(LDX11,*), X21(LDX21,*)
247 INTEGER IWORK(*)
248* ..
249*
250* =====================================================================
251*
252* .. Parameters ..
253 DOUBLE PRECISION ONE, ZERO
254 parameter( one = 1.0d0, zero = 0.0d0 )
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 DOUBLE PRECISION DUM1(1), DUM2(1,1)
267* ..
268* .. External Subroutines ..
269 EXTERNAL dbbcsd, dcopy, dlacpy, dlapmr, dlapmt, 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, theta,
374 $ dum1, u1, ldu1, u2, ldu2, v1t, ldv1t,
375 $ dum2, 1, dum1, dum1, dum1,
376 $ dum1, dum1, dum1, dum1,
377 $ dum1, work(1), -1, childinfo )
378 lbbcsd = int( work(1) )
379 ELSE IF( r .EQ. p ) THEN
380 CALL dorbdb2( m, p, q, x11, ldx11, x21, ldx21, theta,
381 $ dum1, dum1, dum1, dum1,
382 $ work(1), -1, childinfo )
383 lorbdb = int( work(1) )
384 IF( wantu1 .AND. p .GT. 0 ) THEN
385 CALL dorgqr( 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 dorgqr( m-p, m-p, q, u2, ldu2, dum1, work(1),
392 $ -1, 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 dorglq( 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 dbbcsd( jobv1t, 'N', jobu1, jobu2, 'T', m, q, p, theta,
403 $ dum1, v1t, ldv1t, dum2, 1, u1, ldu1,
404 $ u2, ldu2, dum1, dum1, dum1,
405 $ dum1, dum1, dum1, dum1,
406 $ dum1, work(1), -1, childinfo )
407 lbbcsd = int( work(1) )
408 ELSE IF( r .EQ. m-p ) THEN
409 CALL dorbdb3( m, p, q, x11, ldx11, x21, ldx21, theta,
410 $ dum1, dum1, dum1, dum1,
411 $ work(1), -1, childinfo )
412 lorbdb = int( work(1) )
413 IF( wantu1 .AND. p .GT. 0 ) THEN
414 CALL dorgqr( 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 dorgqr( m-p-1, m-p-1, m-p-1, u2(2,2), ldu2,
421 $ dum1, 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 dorglq( 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 dbbcsd( 'N', jobv1t, jobu2, jobu1, 'T', m, m-q, m-p,
432 $ theta, dum1, dum2, 1, v1t, ldv1t, u2,
433 $ ldu2, u1, ldu1, dum1, dum1, dum1,
434 $ dum1, dum1, dum1, dum1,
435 $ dum1, work(1), -1, childinfo )
436 lbbcsd = int( work(1) )
437 ELSE
438 CALL dorbdb4( m, p, q, x11, ldx11, x21, ldx21, theta,
439 $ dum1, dum1, dum1, dum1,
440 $ dum1, work(1), -1, childinfo )
441 lorbdb = m + int( work(1) )
442 IF( wantu1 .AND. p .GT. 0 ) THEN
443 CALL dorgqr( 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 dorgqr( 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 dorglq( 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 dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
461 $ theta, dum1, u2, ldu2, u1, ldu1, dum2,
462 $ 1, v1t, ldv1t, dum1, dum1, dum1,
463 $ dum1, dum1, dum1, dum1,
464 $ dum1, work(1), -1, 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( 'DORCSD2BY1', -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 dorbdb1( 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 dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
506 CALL dorgqr( 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 dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
511 CALL dorgqr( 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 dlacpy( 'U', q-1, q-1, x21(1,2), ldx21, v1t(2,2),
521 $ ldv1t )
522 CALL dorglq( 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 dbbcsd( 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),
531 $ work(ib12d), work(ib12e), work(ib21d),
532 $ work(ib21e), work(ib22d), work(ib22e),
533 $ work(ibbcsd), lbbcsd, 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 dlapmt( .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 dorbdb2( 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 dlacpy( 'L', p-1, p-1, x11(2,1), ldx11, u1(2,2), ldu1 )
566 CALL dorgqr( 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 dlacpy( 'L', m-p, q, x21, ldx21, u2, ldu2 )
571 CALL dorgqr( 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 dlacpy( 'U', p, q, x11, ldx11, v1t, ldv1t )
576 CALL dorglq( 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 dbbcsd( 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 dlapmt( .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 dorbdb3( 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 dlacpy( 'L', p, q, x11, ldx11, u1, ldu1 )
615 CALL dorgqr( 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 dlacpy( 'L', m-p-1, m-p-1, x21(2,1), ldx21, u2(2,2),
625 $ ldu2 )
626 CALL dorgqr( 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 dlacpy( 'U', m-p, q, x21, ldx21, v1t, ldv1t )
631 CALL dorglq( 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 dbbcsd( '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 dlapmt( .false., p, q, u1, ldu1, iwork )
656 END IF
657 IF( wantv1t ) THEN
658 CALL dlapmr( .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 dorbdb4( 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 dcopy( m-p, work(iorbdb+p), 1, u2, 1 )
676 END IF
677 IF( wantu1 .AND. p .GT. 0 ) THEN
678 CALL dcopy( p, work(iorbdb), 1, u1, 1 )
679 DO j = 2, p
680 u1(1,j) = zero
681 END DO
682 CALL dlacpy( 'L', p-1, m-q-1, x11(2,1), ldx11, u1(2,2),
683 $ ldu1 )
684 CALL dorgqr( 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 dlacpy( 'L', m-p-1, m-q-1, x21(2,1), ldx21, u2(2,2),
692 $ ldu2 )
693 CALL dorgqr( 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 dlacpy( 'U', m-q, q, x21, ldx21, v1t, ldv1t )
698 CALL dlacpy( '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 dlacpy( 'U', -p+q, q-p, x21(m-q+1,p+1), ldx21,
701 $ v1t(p+1,p+1), ldv1t )
702 CALL dorglq( 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 dbbcsd( jobu2, jobu1, 'N', jobv1t, 'N', m, m-p, m-q,
709 $ theta, work(iphi), u2, ldu2, u1, ldu1, dum1,
710 $ 1, v1t, ldv1t, work(ib11d), work(ib11e),
711 $ work(ib12d), work(ib12e), work(ib21d),
712 $ work(ib21e), work(ib22d), work(ib22e),
713 $ work(ibbcsd), lbbcsd, 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 dlapmt( .false., p, p, u1, ldu1, iwork )
727 END IF
728 IF( wantv1t ) THEN
729 CALL dlapmr( .false., p, q, v1t, ldv1t, iwork )
730 END IF
731 END IF
732 END IF
733*
734 RETURN
735*
736* End of DORCSD2BY1
737*
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:332
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:103
dlapmr
subroutine dlapmr(forwrd, m, n, x, ldx, k)
DLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition dlapmr.f:104
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:104
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:203
dorbdb2
subroutine dorbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB2
Definition dorbdb2.f:202
dorbdb3
subroutine dorbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
DORBDB3
Definition dorbdb3.f:201
dorbdb4
subroutine dorbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
DORBDB4
Definition dorbdb4.f:213
dorglq
subroutine dorglq(m, n, k, a, lda, tau, work, lwork, info)
DORGLQ
Definition dorglq.f:127
dorgqr
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:128
Here is the call graph for this function:
Here is the caller graph for this function: