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

◆ 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.
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 284 of file zunbdb.f.

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