LAPACK 3.12.0
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 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,
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, 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*
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:332
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:103
slapmr
subroutine slapmr(forwrd, m, n, x, ldx, k)
SLAPMR rearranges rows of a matrix as specified by a permutation vector.
Definition slapmr.f:104
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:104
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:203
sorbdb2
subroutine sorbdb2(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
SORBDB2
Definition sorbdb2.f:201
sorbdb3
subroutine sorbdb3(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, work, lwork, info)
SORBDB3
Definition sorbdb3.f:202
sorbdb4
subroutine sorbdb4(m, p, q, x11, ldx11, x21, ldx21, theta, phi, taup1, taup2, tauq1, phantom, work, lwork, info)
SORBDB4
Definition sorbdb4.f:214
sorglq
subroutine sorglq(m, n, k, a, lda, tau, work, lwork, info)
SORGLQ
Definition sorglq.f:127
sorgqr
subroutine sorgqr(m, n, k, a, lda, tau, work, lwork, info)
SORGQR
Definition sorgqr.f:128
Here is the call graph for this function:
Here is the caller graph for this function: