LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ zunbdb()

subroutine zunbdb ( character  TRANS,
character  SIGNS,
integer  M,
integer  P,
integer  Q,
complex*16, dimension( ldx11, * )  X11,
integer  LDX11,
complex*16, dimension( ldx12, * )  X12,
integer  LDX12,
complex*16, dimension( ldx21, * )  X21,
integer  LDX21,
complex*16, dimension( ldx22, * )  X22,
integer  LDX22,
double precision, dimension( * )  THETA,
double precision, dimension( * )  PHI,
complex*16, dimension( * )  TAUP1,
complex*16, dimension( * )  TAUP2,
complex*16, dimension( * )  TAUQ1,
complex*16, dimension( * )  TAUQ2,
complex*16, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

ZUNBDB

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

Purpose:
 ZUNBDB 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 ZUNCSD
 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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*16 array, dimension (P)
          The scalar factors of the elementary reflectors that define
          P1.
[out]TAUP2
          TAUP2 is COMPLEX*16 array, dimension (M-P)
          The scalar factors of the elementary reflectors that define
          P2.
[out]TAUQ1
          TAUQ1 is COMPLEX*16 array, dimension (Q)
          The scalar factors of the elementary reflectors that define
          Q1.
[out]TAUQ2
          TAUQ2 is COMPLEX*16 array, dimension (M-Q)
          The scalar factors of the elementary reflectors that define
          Q2.
[out]WORK
          WORK is COMPLEX*16 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
December 2016
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 ZUNCSD for details.

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

Definition at line 289 of file zunbdb.f.

289 *
290 * -- LAPACK computational routine (version 3.7.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 * December 2016
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  DOUBLE PRECISION phi( * ), theta( * )
302  COMPLEX*16 taup1( * ), taup2( * ), tauq1( * ), tauq2( * ),
303  $ work( * ), x11( ldx11, * ), x12( ldx12, * ),
304  $ x21( ldx21, * ), x22( ldx22, * )
305 * ..
306 *
307 * ====================================================================
308 *
309 * .. Parameters ..
310  DOUBLE PRECISION realone
311  parameter( realone = 1.0d0 )
312  COMPLEX*16 one
313  parameter( one = (1.0d0,0.0d0) )
314 * ..
315 * .. Local Scalars ..
316  LOGICAL colmajor, lquery
317  INTEGER i, lworkmin, lworkopt
318  DOUBLE PRECISION z1, z2, z3, z4
319 * ..
320 * .. External Subroutines ..
321  EXTERNAL zaxpy, zlarf, zlarfgp, zscal, xerbla
322  EXTERNAL zlacgv
323 *
324 * ..
325 * .. External Functions ..
326  DOUBLE PRECISION dznrm2
327  LOGICAL lsame
328  EXTERNAL dznrm2, lsame
329 * ..
330 * .. Intrinsic Functions
331  INTRINSIC atan2, cos, max, min, sin
332  INTRINSIC dcmplx, dconjg
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 zscal( p-i+1, dcmplx( z1, 0.0d0 ), x11(i,i), 1 )
405  ELSE
406  CALL zscal( p-i+1, dcmplx( z1*cos(phi(i-1)), 0.0d0 ),
407  $ x11(i,i), 1 )
408  CALL zaxpy( p-i+1, dcmplx( -z1*z3*z4*sin(phi(i-1)),
409  $ 0.0d0 ), x12(i,i-1), 1, x11(i,i), 1 )
410  END IF
411  IF( i .EQ. 1 ) THEN
412  CALL zscal( m-p-i+1, dcmplx( z2, 0.0d0 ), x21(i,i), 1 )
413  ELSE
414  CALL zscal( m-p-i+1, dcmplx( z2*cos(phi(i-1)), 0.0d0 ),
415  $ x21(i,i), 1 )
416  CALL zaxpy( m-p-i+1, dcmplx( -z2*z3*z4*sin(phi(i-1)),
417  $ 0.0d0 ), x22(i,i-1), 1, x21(i,i), 1 )
418  END IF
419 *
420  theta(i) = atan2( dznrm2( m-p-i+1, x21(i,i), 1 ),
421  $ dznrm2( p-i+1, x11(i,i), 1 ) )
422 *
423  IF( p .GT. i ) THEN
424  CALL zlarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
425  ELSE IF ( p .EQ. i ) THEN
426  CALL zlarfgp( 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 zlarfgp( 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 zlarfgp( 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 zlarf( 'L', p-i+1, q-i, x11(i,i), 1,
440  $ dconjg(taup1(i)), x11(i,i+1), ldx11, work )
441  CALL zlarf( 'L', m-p-i+1, q-i, x21(i,i), 1,
442  $ dconjg(taup2(i)), x21(i,i+1), ldx21, work )
443  END IF
444  IF ( m-q+1 .GT. i ) THEN
445  CALL zlarf( 'L', p-i+1, m-q-i+1, x11(i,i), 1,
446  $ dconjg(taup1(i)), x12(i,i), ldx12, work )
447  CALL zlarf( 'L', m-p-i+1, m-q-i+1, x21(i,i), 1,
448  $ dconjg(taup2(i)), x22(i,i), ldx22, work )
449  END IF
450 *
451  IF( i .LT. q ) THEN
452  CALL zscal( q-i, dcmplx( -z1*z3*sin(theta(i)), 0.0d0 ),
453  $ x11(i,i+1), ldx11 )
454  CALL zaxpy( q-i, dcmplx( z2*z3*cos(theta(i)), 0.0d0 ),
455  $ x21(i,i+1), ldx21, x11(i,i+1), ldx11 )
456  END IF
457  CALL zscal( m-q-i+1, dcmplx( -z1*z4*sin(theta(i)), 0.0d0 ),
458  $ x12(i,i), ldx12 )
459  CALL zaxpy( m-q-i+1, dcmplx( z2*z4*cos(theta(i)), 0.0d0 ),
460  $ x22(i,i), ldx22, x12(i,i), ldx12 )
461 *
462  IF( i .LT. q )
463  $ phi(i) = atan2( dznrm2( q-i, x11(i,i+1), ldx11 ),
464  $ dznrm2( m-q-i+1, x12(i,i), ldx12 ) )
465 *
466  IF( i .LT. q ) THEN
467  CALL zlacgv( q-i, x11(i,i+1), ldx11 )
468  IF ( i .EQ. q-1 ) THEN
469  CALL zlarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
470  $ tauq1(i) )
471  ELSE
472  CALL zlarfgp( 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 zlacgv( m-q-i+1, x12(i,i), ldx12 )
479  IF ( m-q .EQ. i ) THEN
480  CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
481  $ tauq2(i) )
482  ELSE
483  CALL zlarfgp( 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 zlarf( 'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
491  $ x11(i+1,i+1), ldx11, work )
492  CALL zlarf( '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 zlarf( '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 zlarf( '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 zlacgv( q-i, x11(i,i+1), ldx11 )
506  CALL zlacgv( 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 zscal( m-q-i+1, dcmplx( -z1*z4, 0.0d0 ), x12(i,i),
515  $ ldx12 )
516  CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
517  IF ( i .GE. m-q ) THEN
518  CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
519  $ tauq2(i) )
520  ELSE
521  CALL zlarfgp( 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 zlarf( '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 zlarf( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
532  $ tauq2(i), x22(q+1,i), ldx22, work )
533 *
534  CALL zlacgv( 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 zscal( m-p-q-i+1, dcmplx( z2*z4, 0.0d0 ),
543  $ x22(q+i,p+i), ldx22 )
544  CALL zlacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
545  CALL zlarfgp( 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 zlarf( '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 zlacgv( 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 zscal( p-i+1, dcmplx( z1, 0.0d0 ), x11(i,i),
563  $ ldx11 )
564  ELSE
565  CALL zscal( p-i+1, dcmplx( z1*cos(phi(i-1)), 0.0d0 ),
566  $ x11(i,i), ldx11 )
567  CALL zaxpy( p-i+1, dcmplx( -z1*z3*z4*sin(phi(i-1)),
568  $ 0.0d0 ), x12(i-1,i), ldx12, x11(i,i), ldx11 )
569  END IF
570  IF( i .EQ. 1 ) THEN
571  CALL zscal( m-p-i+1, dcmplx( z2, 0.0d0 ), x21(i,i),
572  $ ldx21 )
573  ELSE
574  CALL zscal( m-p-i+1, dcmplx( z2*cos(phi(i-1)), 0.0d0 ),
575  $ x21(i,i), ldx21 )
576  CALL zaxpy( m-p-i+1, dcmplx( -z2*z3*z4*sin(phi(i-1)),
577  $ 0.0d0 ), x22(i-1,i), ldx22, x21(i,i), ldx21 )
578  END IF
579 *
580  theta(i) = atan2( dznrm2( m-p-i+1, x21(i,i), ldx21 ),
581  $ dznrm2( p-i+1, x11(i,i), ldx11 ) )
582 *
583  CALL zlacgv( p-i+1, x11(i,i), ldx11 )
584  CALL zlacgv( m-p-i+1, x21(i,i), ldx21 )
585 *
586  CALL zlarfgp( 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 zlarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
590  $ taup2(i) )
591  ELSE
592  CALL zlarfgp( 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 zlarf( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
598  $ x11(i+1,i), ldx11, work )
599  CALL zlarf( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11, taup1(i),
600  $ x12(i,i), ldx12, work )
601  CALL zlarf( 'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
602  $ x21(i+1,i), ldx21, work )
603  CALL zlarf( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
604  $ taup2(i), x22(i,i), ldx22, work )
605 *
606  CALL zlacgv( p-i+1, x11(i,i), ldx11 )
607  CALL zlacgv( m-p-i+1, x21(i,i), ldx21 )
608 *
609  IF( i .LT. q ) THEN
610  CALL zscal( q-i, dcmplx( -z1*z3*sin(theta(i)), 0.0d0 ),
611  $ x11(i+1,i), 1 )
612  CALL zaxpy( q-i, dcmplx( z2*z3*cos(theta(i)), 0.0d0 ),
613  $ x21(i+1,i), 1, x11(i+1,i), 1 )
614  END IF
615  CALL zscal( m-q-i+1, dcmplx( -z1*z4*sin(theta(i)), 0.0d0 ),
616  $ x12(i,i), 1 )
617  CALL zaxpy( m-q-i+1, dcmplx( z2*z4*cos(theta(i)), 0.0d0 ),
618  $ x22(i,i), 1, x12(i,i), 1 )
619 *
620  IF( i .LT. q )
621  $ phi(i) = atan2( dznrm2( q-i, x11(i+1,i), 1 ),
622  $ dznrm2( m-q-i+1, x12(i,i), 1 ) )
623 *
624  IF( i .LT. q ) THEN
625  CALL zlarfgp( q-i, x11(i+1,i), x11(i+2,i), 1, tauq1(i) )
626  x11(i+1,i) = one
627  END IF
628  CALL zlarfgp( 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 zlarf( 'L', q-i, p-i, x11(i+1,i), 1,
633  $ dconjg(tauq1(i)), x11(i+1,i+1), ldx11, work )
634  CALL zlarf( 'L', q-i, m-p-i, x11(i+1,i), 1,
635  $ dconjg(tauq1(i)), x21(i+1,i+1), ldx21, work )
636  END IF
637  CALL zlarf( 'L', m-q-i+1, p-i, x12(i,i), 1,
638  $ dconjg(tauq2(i)), x12(i,i+1), ldx12, work )
639  IF ( m-p .GT. i ) THEN
640  CALL zlarf( 'L', m-q-i+1, m-p-i, x12(i,i), 1,
641  $ dconjg(tauq2(i)), x22(i,i+1), ldx22, work )
642  END IF
643 *
644  END DO
645 *
646 * Reduce columns Q + 1, ..., P of X12, X22
647 *
648  DO i = q + 1, p
649 *
650  CALL zscal( m-q-i+1, dcmplx( -z1*z4, 0.0d0 ), x12(i,i), 1 )
651  CALL zlarfgp( 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 zlarf( 'L', m-q-i+1, p-i, x12(i,i), 1,
656  $ dconjg(tauq2(i)), x12(i,i+1), ldx12, work )
657  END IF
658  IF( m-p-q .GE. 1 )
659  $ CALL zlarf( 'L', m-q-i+1, m-p-q, x12(i,i), 1,
660  $ dconjg(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 zscal( m-p-q-i+1, dcmplx( z2*z4, 0.0d0 ),
669  $ x22(p+i,q+i), 1 )
670  CALL zlarfgp( 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 *
674  IF ( m-p-q .NE. i ) THEN
675  CALL zlarf( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
676  $ dconjg(tauq2(p+i)), x22(p+i,q+i+1), ldx22,
677  $ work )
678  END IF
679 *
680  END DO
681 *
682  END IF
683 *
684  RETURN
685 *
686 * End of ZUNBDB
687 *
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
Definition: zlarf.f:130
double precision function dznrm2(N, X, INCX)
DZNRM2
Definition: dznrm2.f:77
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55
subroutine zlacgv(N, X, INCX)
ZLACGV conjugates a complex vector.
Definition: zlacgv.f:76
subroutine zlarfgp(N, ALPHA, X, INCX, TAU)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: zlarfgp.f:106
subroutine zaxpy(N, ZA, ZX, INCX, ZY, INCY)
ZAXPY
Definition: zaxpy.f:90
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
Definition: zscal.f:80
Here is the call graph for this function:
Here is the caller graph for this function: