LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine cunbdb ( character  TRANS,
character  SIGNS,
integer  M,
integer  P,
integer  Q,
complex, dimension( ldx11, * )  X11,
integer  LDX11,
complex, dimension( ldx12, * )  X12,
integer  LDX12,
complex, dimension( ldx21, * )  X21,
integer  LDX21,
complex, dimension( ldx22, * )  X22,
integer  LDX22,
real, dimension( * )  THETA,
real, dimension( * )  PHI,
complex, dimension( * )  TAUP1,
complex, dimension( * )  TAUP2,
complex, dimension( * )  TAUQ1,
complex, dimension( * )  TAUQ2,
complex, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

CUNBDB

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

Purpose:
 CUNBDB simultaneously bidiagonalizes the blocks of an M-by-M
 partitioned unitary matrix X:

                                 [ B11 | B12 0  0 ]
     [ X11 | X12 ]   [ P1 |    ] [  0  |  0 -I  0 ] [ Q1 |    ]**H
 X = [-----------] = [---------] [----------------] [---------]   .
     [ X21 | X22 ]   [    | P2 ] [ B21 | B22 0  0 ] [    | Q2 ]
                                 [  0  |  0  0  I ]

 X11 is P-by-Q. Q must be no larger than P, M-P, or M-Q. (If this is
 not the case, then X must be transposed and/or permuted. This can be
 done in constant time using the TRANS and SIGNS options. See CUNCSD
 for details.)

 The unitary matrices P1, P2, Q1, and Q2 are P-by-P, (M-P)-by-
 (M-P), Q-by-Q, and (M-Q)-by-(M-Q), respectively. They are
 represented implicitly by Householder vectors.

 B11, B12, B21, and B22 are Q-by-Q bidiagonal matrices represented
 implicitly by angles THETA, PHI.
Parameters
[in]TRANS
          TRANS is CHARACTER
          = 'T':      X, U1, U2, V1T, and V2T are stored in row-major
                      order;
          otherwise:  X, U1, U2, V1T, and V2T are stored in column-
                      major order.
[in]SIGNS
          SIGNS is CHARACTER
          = 'O':      The lower-left block is made nonpositive (the
                      "other" convention);
          otherwise:  The upper-right block is made nonpositive (the
                      "default" convention).
[in]M
          M is INTEGER
          The number of rows and columns in X.
[in]P
          P is INTEGER
          The number of rows in X11 and X12. 0 <= P <= M.
[in]Q
          Q is INTEGER
          The number of columns in X11 and X21. 0 <= Q <=
          MIN(P,M-P,M-Q).
[in,out]X11
          X11 is COMPLEX array, dimension (LDX11,Q)
          On entry, the top-left block of the unitary matrix to be
          reduced. On exit, the form depends on TRANS:
          If TRANS = 'N', then
             the columns of tril(X11) specify reflectors for P1,
             the rows of triu(X11,1) specify reflectors for Q1;
          else TRANS = 'T', and
             the rows of triu(X11) specify reflectors for P1,
             the columns of tril(X11,-1) specify reflectors for Q1.
[in]LDX11
          LDX11 is INTEGER
          The leading dimension of X11. If TRANS = 'N', then LDX11 >=
          P; else LDX11 >= Q.
[in,out]X12
          X12 is COMPLEX array, dimension (LDX12,M-Q)
          On entry, the top-right block of the unitary matrix to
          be reduced. On exit, the form depends on TRANS:
          If TRANS = 'N', then
             the rows of triu(X12) specify the first P reflectors for
             Q2;
          else TRANS = 'T', and
             the columns of tril(X12) specify the first P reflectors
             for Q2.
[in]LDX12
          LDX12 is INTEGER
          The leading dimension of X12. If TRANS = 'N', then LDX12 >=
          P; else LDX11 >= M-Q.
[in,out]X21
          X21 is COMPLEX array, dimension (LDX21,Q)
          On entry, the bottom-left block of the unitary matrix to
          be reduced. On exit, the form depends on TRANS:
          If TRANS = 'N', then
             the columns of tril(X21) specify reflectors for P2;
          else TRANS = 'T', and
             the rows of triu(X21) specify reflectors for P2.
[in]LDX21
          LDX21 is INTEGER
          The leading dimension of X21. If TRANS = 'N', then LDX21 >=
          M-P; else LDX21 >= Q.
[in,out]X22
          X22 is COMPLEX array, dimension (LDX22,M-Q)
          On entry, the bottom-right block of the unitary matrix to
          be reduced. On exit, the form depends on TRANS:
          If TRANS = 'N', then
             the rows of triu(X22(Q+1:M-P,P+1:M-Q)) specify the last
             M-P-Q reflectors for Q2,
          else TRANS = 'T', and
             the columns of tril(X22(P+1:M-Q,Q+1:M-P)) specify the last
             M-P-Q reflectors for P2.
[in]LDX22
          LDX22 is INTEGER
          The leading dimension of X22. If TRANS = 'N', then LDX22 >=
          M-P; else LDX22 >= M-Q.
[out]THETA
          THETA is REAL array, dimension (Q)
          The entries of the bidiagonal blocks B11, B12, B21, B22 can
          be computed from the angles THETA and PHI. See Further
          Details.
[out]PHI
          PHI is REAL array, dimension (Q-1)
          The entries of the bidiagonal blocks B11, B12, B21, B22 can
          be computed from the angles THETA and PHI. See Further
          Details.
[out]TAUP1
          TAUP1 is COMPLEX array, dimension (P)
          The scalar factors of the elementary reflectors that define
          P1.
[out]TAUP2
          TAUP2 is COMPLEX array, dimension (M-P)
          The scalar factors of the elementary reflectors that define
          P2.
[out]TAUQ1
          TAUQ1 is COMPLEX array, dimension (Q)
          The scalar factors of the elementary reflectors that define
          Q1.
[out]TAUQ2
          TAUQ2 is COMPLEX array, dimension (M-Q)
          The scalar factors of the elementary reflectors that define
          Q2.
[out]WORK
          WORK is COMPLEX array, dimension (LWORK)
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK. LWORK >= M-Q.

          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]INFO
          INFO is INTEGER
          = 0:  successful exit.
          < 0:  if INFO = -i, the i-th argument had an illegal value.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2013
Further Details:
  The bidiagonal blocks B11, B12, B21, and B22 are represented
  implicitly by angles THETA(1), ..., THETA(Q) and PHI(1), ...,
  PHI(Q-1). B11 and B21 are upper bidiagonal, while B21 and B22 are
  lower bidiagonal. Every entry in each bidiagonal band is a product
  of a sine or cosine of a THETA with a sine or cosine of a PHI. See
  [1] or CUNCSD for details.

  P1, P2, Q1, and Q2 are represented as products of elementary
  reflectors. See CUNCSD for details on generating P1, P2, Q1, and Q2
  using CUNGQR and CUNGLQ.
References:
[1] Brian D. Sutton. Computing the complete CS decomposition. Numer. Algorithms, 50(1):33-65, 2009.

Definition at line 289 of file cunbdb.f.

289 *
290 * -- LAPACK computational routine (version 3.5.0) --
291 * -- LAPACK is a software package provided by Univ. of Tennessee, --
292 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
293 * November 2013
294 *
295 * .. Scalar Arguments ..
296  CHARACTER signs, trans
297  INTEGER info, ldx11, ldx12, ldx21, ldx22, lwork, m, p,
298  $ q
299 * ..
300 * .. Array Arguments ..
301  REAL phi( * ), theta( * )
302  COMPLEX taup1( * ), taup2( * ), tauq1( * ), tauq2( * ),
303  $ work( * ), x11( ldx11, * ), x12( ldx12, * ),
304  $ x21( ldx21, * ), x22( ldx22, * )
305 * ..
306 *
307 * ====================================================================
308 *
309 * .. Parameters ..
310  REAL realone
311  parameter ( realone = 1.0e0 )
312  COMPLEX one
313  parameter ( one = (1.0e0,0.0e0) )
314 * ..
315 * .. Local Scalars ..
316  LOGICAL colmajor, lquery
317  INTEGER i, lworkmin, lworkopt
318  REAL z1, z2, z3, z4
319 * ..
320 * .. External Subroutines ..
321  EXTERNAL caxpy, clarf, clarfgp, cscal, xerbla
322  EXTERNAL clacgv
323 *
324 * ..
325 * .. External Functions ..
326  REAL scnrm2
327  LOGICAL lsame
328  EXTERNAL scnrm2, lsame
329 * ..
330 * .. Intrinsic Functions
331  INTRINSIC atan2, cos, max, min, sin
332  INTRINSIC cmplx, conjg
333 * ..
334 * .. Executable Statements ..
335 *
336 * Test input arguments
337 *
338  info = 0
339  colmajor = .NOT. lsame( trans, 'T' )
340  IF( .NOT. lsame( signs, 'O' ) ) THEN
341  z1 = realone
342  z2 = realone
343  z3 = realone
344  z4 = realone
345  ELSE
346  z1 = realone
347  z2 = -realone
348  z3 = realone
349  z4 = -realone
350  END IF
351  lquery = lwork .EQ. -1
352 *
353  IF( m .LT. 0 ) THEN
354  info = -3
355  ELSE IF( p .LT. 0 .OR. p .GT. m ) THEN
356  info = -4
357  ELSE IF( q .LT. 0 .OR. q .GT. p .OR. q .GT. m-p .OR.
358  $ q .GT. m-q ) THEN
359  info = -5
360  ELSE IF( colmajor .AND. ldx11 .LT. max( 1, p ) ) THEN
361  info = -7
362  ELSE IF( .NOT.colmajor .AND. ldx11 .LT. max( 1, q ) ) THEN
363  info = -7
364  ELSE IF( colmajor .AND. ldx12 .LT. max( 1, p ) ) THEN
365  info = -9
366  ELSE IF( .NOT.colmajor .AND. ldx12 .LT. max( 1, m-q ) ) THEN
367  info = -9
368  ELSE IF( colmajor .AND. ldx21 .LT. max( 1, m-p ) ) THEN
369  info = -11
370  ELSE IF( .NOT.colmajor .AND. ldx21 .LT. max( 1, q ) ) THEN
371  info = -11
372  ELSE IF( colmajor .AND. ldx22 .LT. max( 1, m-p ) ) THEN
373  info = -13
374  ELSE IF( .NOT.colmajor .AND. ldx22 .LT. max( 1, m-q ) ) THEN
375  info = -13
376  END IF
377 *
378 * Compute workspace
379 *
380  IF( info .EQ. 0 ) THEN
381  lworkopt = m - q
382  lworkmin = m - q
383  work(1) = lworkopt
384  IF( lwork .LT. lworkmin .AND. .NOT. lquery ) THEN
385  info = -21
386  END IF
387  END IF
388  IF( info .NE. 0 ) THEN
389  CALL xerbla( 'xORBDB', -info )
390  RETURN
391  ELSE IF( lquery ) THEN
392  RETURN
393  END IF
394 *
395 * Handle column-major and row-major separately
396 *
397  IF( colmajor ) THEN
398 *
399 * Reduce columns 1, ..., Q of X11, X12, X21, and X22
400 *
401  DO i = 1, q
402 *
403  IF( i .EQ. 1 ) THEN
404  CALL cscal( p-i+1, cmplx( z1, 0.0e0 ), x11(i,i), 1 )
405  ELSE
406  CALL cscal( p-i+1, cmplx( z1*cos(phi(i-1)), 0.0e0 ),
407  $ x11(i,i), 1 )
408  CALL caxpy( p-i+1, cmplx( -z1*z3*z4*sin(phi(i-1)),
409  $ 0.0e0 ), x12(i,i-1), 1, x11(i,i), 1 )
410  END IF
411  IF( i .EQ. 1 ) THEN
412  CALL cscal( m-p-i+1, cmplx( z2, 0.0e0 ), x21(i,i), 1 )
413  ELSE
414  CALL cscal( m-p-i+1, cmplx( z2*cos(phi(i-1)), 0.0e0 ),
415  $ x21(i,i), 1 )
416  CALL caxpy( m-p-i+1, cmplx( -z2*z3*z4*sin(phi(i-1)),
417  $ 0.0e0 ), x22(i,i-1), 1, x21(i,i), 1 )
418  END IF
419 *
420  theta(i) = atan2( scnrm2( m-p-i+1, x21(i,i), 1 ),
421  $ scnrm2( p-i+1, x11(i,i), 1 ) )
422 *
423  IF( p .GT. i ) THEN
424  CALL clarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
425  ELSE IF ( p .EQ. i ) THEN
426  CALL clarfgp( p-i+1, x11(i,i), x11(i,i), 1, taup1(i) )
427  END IF
428  x11(i,i) = one
429  IF ( m-p .GT. i ) THEN
430  CALL clarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
431  $ taup2(i) )
432  ELSE IF ( m-p .EQ. i ) THEN
433  CALL clarfgp( m-p-i+1, x21(i,i), x21(i,i), 1,
434  $ taup2(i) )
435  END IF
436  x21(i,i) = one
437 *
438  IF ( q .GT. i ) THEN
439  CALL clarf( 'L', p-i+1, q-i, x11(i,i), 1,
440  $ conjg(taup1(i)), x11(i,i+1), ldx11, work )
441  CALL clarf( 'L', m-p-i+1, q-i, x21(i,i), 1,
442  $ conjg(taup2(i)), x21(i,i+1), ldx21, work )
443  END IF
444  IF ( m-q+1 .GT. i ) THEN
445  CALL clarf( 'L', p-i+1, m-q-i+1, x11(i,i), 1,
446  $ conjg(taup1(i)), x12(i,i), ldx12, work )
447  CALL clarf( 'L', m-p-i+1, m-q-i+1, x21(i,i), 1,
448  $ conjg(taup2(i)), x22(i,i), ldx22, work )
449  END IF
450 *
451  IF( i .LT. q ) THEN
452  CALL cscal( q-i, cmplx( -z1*z3*sin(theta(i)), 0.0e0 ),
453  $ x11(i,i+1), ldx11 )
454  CALL caxpy( q-i, cmplx( z2*z3*cos(theta(i)), 0.0e0 ),
455  $ x21(i,i+1), ldx21, x11(i,i+1), ldx11 )
456  END IF
457  CALL cscal( m-q-i+1, cmplx( -z1*z4*sin(theta(i)), 0.0e0 ),
458  $ x12(i,i), ldx12 )
459  CALL caxpy( m-q-i+1, cmplx( z2*z4*cos(theta(i)), 0.0e0 ),
460  $ x22(i,i), ldx22, x12(i,i), ldx12 )
461 *
462  IF( i .LT. q )
463  $ phi(i) = atan2( scnrm2( q-i, x11(i,i+1), ldx11 ),
464  $ scnrm2( m-q-i+1, x12(i,i), ldx12 ) )
465 *
466  IF( i .LT. q ) THEN
467  CALL clacgv( q-i, x11(i,i+1), ldx11 )
468  IF ( i .EQ. q-1 ) THEN
469  CALL clarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
470  $ tauq1(i) )
471  ELSE
472  CALL clarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
473  $ tauq1(i) )
474  END IF
475  x11(i,i+1) = one
476  END IF
477  IF ( m-q+1 .GT. i ) THEN
478  CALL clacgv( m-q-i+1, x12(i,i), ldx12 )
479  IF ( m-q .EQ. i ) THEN
480  CALL clarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
481  $ tauq2(i) )
482  ELSE
483  CALL clarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
484  $ tauq2(i) )
485  END IF
486  END IF
487  x12(i,i) = one
488 *
489  IF( i .LT. q ) THEN
490  CALL clarf( 'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
491  $ x11(i+1,i+1), ldx11, work )
492  CALL clarf( 'R', m-p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
493  $ x21(i+1,i+1), ldx21, work )
494  END IF
495  IF ( p .GT. i ) THEN
496  CALL clarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
497  $ x12(i+1,i), ldx12, work )
498  END IF
499  IF ( m-p .GT. i ) THEN
500  CALL clarf( 'R', m-p-i, m-q-i+1, x12(i,i), ldx12,
501  $ tauq2(i), x22(i+1,i), ldx22, work )
502  END IF
503 *
504  IF( i .LT. q )
505  $ CALL clacgv( q-i, x11(i,i+1), ldx11 )
506  CALL clacgv( m-q-i+1, x12(i,i), ldx12 )
507 *
508  END DO
509 *
510 * Reduce columns Q + 1, ..., P of X12, X22
511 *
512  DO i = q + 1, p
513 *
514  CALL cscal( m-q-i+1, cmplx( -z1*z4, 0.0e0 ), x12(i,i),
515  $ ldx12 )
516  CALL clacgv( m-q-i+1, x12(i,i), ldx12 )
517  IF ( i .GE. m-q ) THEN
518  CALL clarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
519  $ tauq2(i) )
520  ELSE
521  CALL clarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
522  $ tauq2(i) )
523  END IF
524  x12(i,i) = one
525 *
526  IF ( p .GT. i ) THEN
527  CALL clarf( 'R', p-i, m-q-i+1, x12(i,i), ldx12, tauq2(i),
528  $ x12(i+1,i), ldx12, work )
529  END IF
530  IF( m-p-q .GE. 1 )
531  $ CALL clarf( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
532  $ tauq2(i), x22(q+1,i), ldx22, work )
533 *
534  CALL clacgv( m-q-i+1, x12(i,i), ldx12 )
535 *
536  END DO
537 *
538 * Reduce columns P + 1, ..., M - Q of X12, X22
539 *
540  DO i = 1, m - p - q
541 *
542  CALL cscal( m-p-q-i+1, cmplx( z2*z4, 0.0e0 ),
543  $ x22(q+i,p+i), ldx22 )
544  CALL clacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
545  CALL clarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
546  $ ldx22, tauq2(p+i) )
547  x22(q+i,p+i) = one
548  CALL clarf( 'R', m-p-q-i, m-p-q-i+1, x22(q+i,p+i), ldx22,
549  $ tauq2(p+i), x22(q+i+1,p+i), ldx22, work )
550 *
551  CALL clacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
552 *
553  END DO
554 *
555  ELSE
556 *
557 * Reduce columns 1, ..., Q of X11, X12, X21, X22
558 *
559  DO i = 1, q
560 *
561  IF( i .EQ. 1 ) THEN
562  CALL cscal( p-i+1, cmplx( z1, 0.0e0 ), x11(i,i),
563  $ ldx11 )
564  ELSE
565  CALL cscal( p-i+1, cmplx( z1*cos(phi(i-1)), 0.0e0 ),
566  $ x11(i,i), ldx11 )
567  CALL caxpy( p-i+1, cmplx( -z1*z3*z4*sin(phi(i-1)),
568  $ 0.0e0 ), x12(i-1,i), ldx12, x11(i,i), ldx11 )
569  END IF
570  IF( i .EQ. 1 ) THEN
571  CALL cscal( m-p-i+1, cmplx( z2, 0.0e0 ), x21(i,i),
572  $ ldx21 )
573  ELSE
574  CALL cscal( m-p-i+1, cmplx( z2*cos(phi(i-1)), 0.0e0 ),
575  $ x21(i,i), ldx21 )
576  CALL caxpy( m-p-i+1, cmplx( -z2*z3*z4*sin(phi(i-1)),
577  $ 0.0e0 ), x22(i-1,i), ldx22, x21(i,i), ldx21 )
578  END IF
579 *
580  theta(i) = atan2( scnrm2( m-p-i+1, x21(i,i), ldx21 ),
581  $ scnrm2( p-i+1, x11(i,i), ldx11 ) )
582 *
583  CALL clacgv( p-i+1, x11(i,i), ldx11 )
584  CALL clacgv( m-p-i+1, x21(i,i), ldx21 )
585 *
586  CALL clarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11, taup1(i) )
587  x11(i,i) = one
588  IF ( i .EQ. m-p ) THEN
589  CALL clarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
590  $ taup2(i) )
591  ELSE
592  CALL clarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
593  $ taup2(i) )
594  END IF
595  x21(i,i) = one
596 *
597  CALL clarf( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
598  $ x11(i+1,i), ldx11, work )
599  CALL clarf( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11, taup1(i),
600  $ x12(i,i), ldx12, work )
601  CALL clarf( 'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
602  $ x21(i+1,i), ldx21, work )
603  CALL clarf( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
604  $ taup2(i), x22(i,i), ldx22, work )
605 *
606  CALL clacgv( p-i+1, x11(i,i), ldx11 )
607  CALL clacgv( m-p-i+1, x21(i,i), ldx21 )
608 *
609  IF( i .LT. q ) THEN
610  CALL cscal( q-i, cmplx( -z1*z3*sin(theta(i)), 0.0e0 ),
611  $ x11(i+1,i), 1 )
612  CALL caxpy( q-i, cmplx( z2*z3*cos(theta(i)), 0.0e0 ),
613  $ x21(i+1,i), 1, x11(i+1,i), 1 )
614  END IF
615  CALL cscal( m-q-i+1, cmplx( -z1*z4*sin(theta(i)), 0.0e0 ),
616  $ x12(i,i), 1 )
617  CALL caxpy( m-q-i+1, cmplx( z2*z4*cos(theta(i)), 0.0e0 ),
618  $ x22(i,i), 1, x12(i,i), 1 )
619 *
620  IF( i .LT. q )
621  $ phi(i) = atan2( scnrm2( q-i, x11(i+1,i), 1 ),
622  $ scnrm2( m-q-i+1, x12(i,i), 1 ) )
623 *
624  IF( i .LT. q ) THEN
625  CALL clarfgp( q-i, x11(i+1,i), x11(i+2,i), 1, tauq1(i) )
626  x11(i+1,i) = one
627  END IF
628  CALL clarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
629  x12(i,i) = one
630 *
631  IF( i .LT. q ) THEN
632  CALL clarf( 'L', q-i, p-i, x11(i+1,i), 1,
633  $ conjg(tauq1(i)), x11(i+1,i+1), ldx11, work )
634  CALL clarf( 'L', q-i, m-p-i, x11(i+1,i), 1,
635  $ conjg(tauq1(i)), x21(i+1,i+1), ldx21, work )
636  END IF
637  CALL clarf( 'L', m-q-i+1, p-i, x12(i,i), 1, conjg(tauq2(i)),
638  $ x12(i,i+1), ldx12, work )
639 
640  IF ( m-p .GT. i ) THEN
641  CALL clarf( 'L', m-q-i+1, m-p-i, x12(i,i), 1,
642  $ conjg(tauq2(i)), x22(i,i+1), ldx22, work )
643  END IF
644  END DO
645 *
646 * Reduce columns Q + 1, ..., P of X12, X22
647 *
648  DO i = q + 1, p
649 *
650  CALL cscal( m-q-i+1, cmplx( -z1*z4, 0.0e0 ), x12(i,i), 1 )
651  CALL clarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1, tauq2(i) )
652  x12(i,i) = one
653 *
654  IF ( p .GT. i ) THEN
655  CALL clarf( 'L', m-q-i+1, p-i, x12(i,i), 1,
656  $ conjg(tauq2(i)), x12(i,i+1), ldx12, work )
657  END IF
658  IF( m-p-q .GE. 1 )
659  $ CALL clarf( 'L', m-q-i+1, m-p-q, x12(i,i), 1,
660  $ conjg(tauq2(i)), x22(i,q+1), ldx22, work )
661 *
662  END DO
663 *
664 * Reduce columns P + 1, ..., M - Q of X12, X22
665 *
666  DO i = 1, m - p - q
667 *
668  CALL cscal( m-p-q-i+1, cmplx( z2*z4, 0.0e0 ),
669  $ x22(p+i,q+i), 1 )
670  CALL clarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i), 1,
671  $ tauq2(p+i) )
672  x22(p+i,q+i) = one
673  IF ( m-p-q .NE. i ) THEN
674  CALL clarf( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
675  $ conjg(tauq2(p+i)), x22(p+i,q+i+1), ldx22,
676  $ work )
677  END IF
678  END DO
679 *
680  END IF
681 *
682  RETURN
683 *
684 * End of CUNBDB
685 *
real function scnrm2(N, X, INCX)
SCNRM2
Definition: scnrm2.f:56
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:54
subroutine clarfgp(N, ALPHA, X, INCX, TAU)
CLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: clarfgp.f:106
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:130
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:76
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:53
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: