LAPACK 3.12.1
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
!>                       convention);
!>          otherwise:  The upper-right block is made nonpositive (the
!>                       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 282 of file zunbdb.f.

286*
287* -- LAPACK computational routine --
288* -- LAPACK is a software package provided by Univ. of Tennessee, --
289* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
290*
291* .. Scalar Arguments ..
292 CHARACTER SIGNS, TRANS
293 INTEGER INFO, LDX11, LDX12, LDX21, LDX22, LWORK, M, P,
294 $ Q
295* ..
296* .. Array Arguments ..
297 DOUBLE PRECISION PHI( * ), THETA( * )
298 COMPLEX*16 TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
299 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
300 $ X21( LDX21, * ), X22( LDX22, * )
301* ..
302*
303* ====================================================================
304*
305* .. Parameters ..
306 DOUBLE PRECISION REALONE
307 parameter( realone = 1.0d0 )
308 COMPLEX*16 ONE
309 parameter( one = (1.0d0,0.0d0) )
310* ..
311* .. Local Scalars ..
312 LOGICAL COLMAJOR, LQUERY
313 INTEGER I, LWORKMIN, LWORKOPT
314 DOUBLE PRECISION Z1, Z2, Z3, Z4
315* ..
316* .. External Subroutines ..
317 EXTERNAL zaxpy, zlarf1f, zlarfgp, zscal,
318 $ 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),
410 $ 1 )
411 ELSE
412 CALL zscal( m-p-i+1, dcmplx( z2*cos(phi(i-1)),
413 $ 0.0d0 ),
414 $ x21(i,i), 1 )
415 CALL zaxpy( m-p-i+1, dcmplx( -z2*z3*z4*sin(phi(i-1)),
416 $ 0.0d0 ), x22(i,i-1), 1, x21(i,i), 1 )
417 END IF
418*
419 theta(i) = atan2( dznrm2( m-p-i+1, x21(i,i), 1 ),
420 $ dznrm2( p-i+1, x11(i,i), 1 ) )
421*
422 IF( p .GT. i ) THEN
423 CALL zlarfgp( p-i+1, x11(i,i), x11(i+1,i), 1,
424 $ 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 IF ( m-p .GT. i ) THEN
429 CALL zlarfgp( m-p-i+1, x21(i,i), x21(i+1,i), 1,
430 $ taup2(i) )
431 ELSE IF ( m-p .EQ. i ) THEN
432 CALL zlarfgp( m-p-i+1, x21(i,i), x21(i,i), 1,
433 $ taup2(i) )
434 END IF
435*
436 IF ( q .GT. i ) THEN
437 CALL zlarf1f( 'L', p-i+1, q-i, x11(i,i), 1,
438 $ conjg(taup1(i)), x11(i,i+1), ldx11,
439 $ work )
440 CALL zlarf1f( 'L', m-p-i+1, q-i, x21(i,i), 1,
441 $ conjg(taup2(i)), x21(i,i+1), ldx21,
442 $ work )
443 END IF
444 IF ( m-q+1 .GT. i ) THEN
445 CALL zlarf1f( 'L', p-i+1, m-q-i+1, x11(i,i), 1,
446 $ conjg(taup1(i)), x12(i,i), ldx12, work )
447 CALL zlarf1f( '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 zscal( q-i, dcmplx( -z1*z3*sin(theta(i)),
453 $ 0.0d0 ),
454 $ x11(i,i+1), ldx11 )
455 CALL zaxpy( q-i, dcmplx( z2*z3*cos(theta(i)), 0.0d0 ),
456 $ x21(i,i+1), ldx21, x11(i,i+1), ldx11 )
457 END IF
458 CALL zscal( m-q-i+1, dcmplx( -z1*z4*sin(theta(i)),
459 $ 0.0d0 ),
460 $ x12(i,i), ldx12 )
461 CALL zaxpy( m-q-i+1, dcmplx( z2*z4*cos(theta(i)),
462 $ 0.0d0 ),
463 $ x22(i,i), ldx22, x12(i,i), ldx12 )
464*
465 IF( i .LT. q )
466 $ phi(i) = atan2( dznrm2( q-i, x11(i,i+1), ldx11 ),
467 $ dznrm2( m-q-i+1, x12(i,i), ldx12 ) )
468*
469 IF( i .LT. q ) THEN
470 CALL zlacgv( q-i, x11(i,i+1), ldx11 )
471 IF ( i .EQ. q-1 ) THEN
472 CALL zlarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
473 $ tauq1(i) )
474 ELSE
475 CALL zlarfgp( q-i, x11(i,i+1), x11(i,i+2), ldx11,
476 $ tauq1(i) )
477 END IF
478 END IF
479 IF ( m-q+1 .GT. i ) THEN
480 CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
481 IF ( m-q .EQ. i ) THEN
482 CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
483 $ tauq2(i) )
484 ELSE
485 CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
486 $ tauq2(i) )
487 END IF
488 END IF
489*
490 IF( i .LT. q ) THEN
491 CALL zlarf1f( 'R', p-i, q-i, x11(i,i+1), ldx11,
492 $ tauq1(i),
493 $ x11(i+1,i+1), ldx11, work )
494 CALL zlarf1f( 'R', m-p-i, q-i, x11(i,i+1), ldx11,
495 $ tauq1(i),
496 $ x21(i+1,i+1), ldx21, work )
497 END IF
498 IF ( p .GT. i ) THEN
499 CALL zlarf1f( 'R', p-i, m-q-i+1, x12(i,i), ldx12,
500 $ tauq2(i),
501 $ x12(i+1,i), ldx12, work )
502 END IF
503 IF ( m-p .GT. i ) THEN
504 CALL zlarf1f( 'R', m-p-i, m-q-i+1, x12(i,i), ldx12,
505 $ tauq2(i), x22(i+1,i), ldx22, work )
506 END IF
507*
508 IF( i .LT. q )
509 $ CALL zlacgv( q-i, x11(i,i+1), ldx11 )
510 CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
511*
512 END DO
513*
514* Reduce columns Q + 1, ..., P of X12, X22
515*
516 DO i = q + 1, p
517*
518 CALL zscal( m-q-i+1, dcmplx( -z1*z4, 0.0d0 ), x12(i,i),
519 $ ldx12 )
520 CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
521 IF ( i .GE. m-q ) THEN
522 CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
523 $ tauq2(i) )
524 ELSE
525 CALL zlarfgp( m-q-i+1, x12(i,i), x12(i,i+1), ldx12,
526 $ tauq2(i) )
527 END IF
528*
529 IF ( p .GT. i ) THEN
530 CALL zlarf1f( 'R', p-i, m-q-i+1, x12(i,i), ldx12,
531 $ tauq2(i),
532 $ x12(i+1,i), ldx12, work )
533 END IF
534 IF( m-p-q .GE. 1 )
535 $ CALL zlarf1f( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
536 $ tauq2(i), x22(q+1,i), ldx22, work )
537*
538 CALL zlacgv( m-q-i+1, x12(i,i), ldx12 )
539*
540 END DO
541*
542* Reduce columns P + 1, ..., M - Q of X12, X22
543*
544 DO i = 1, m - p - q
545*
546 CALL zscal( m-p-q-i+1, dcmplx( z2*z4, 0.0d0 ),
547 $ x22(q+i,p+i), ldx22 )
548 CALL zlacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
549 CALL zlarfgp( m-p-q-i+1, x22(q+i,p+i), x22(q+i,p+i+1),
550 $ ldx22, tauq2(p+i) )
551 CALL zlarf1f( 'R', m-p-q-i, m-p-q-i+1, x22(q+i,p+i),
552 $ ldx22, tauq2(p+i), x22(q+i+1,p+i), ldx22,
553 $ work )
554*
555 CALL zlacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
556*
557 END DO
558*
559 ELSE
560*
561* Reduce columns 1, ..., Q of X11, X12, X21, X22
562*
563 DO i = 1, q
564*
565 IF( i .EQ. 1 ) THEN
566 CALL zscal( p-i+1, dcmplx( z1, 0.0d0 ), x11(i,i),
567 $ ldx11 )
568 ELSE
569 CALL zscal( p-i+1, dcmplx( z1*cos(phi(i-1)), 0.0d0 ),
570 $ x11(i,i), ldx11 )
571 CALL zaxpy( p-i+1, dcmplx( -z1*z3*z4*sin(phi(i-1)),
572 $ 0.0d0 ), x12(i-1,i), ldx12, x11(i,i), ldx11 )
573 END IF
574 IF( i .EQ. 1 ) THEN
575 CALL zscal( m-p-i+1, dcmplx( z2, 0.0d0 ), x21(i,i),
576 $ ldx21 )
577 ELSE
578 CALL zscal( m-p-i+1, dcmplx( z2*cos(phi(i-1)),
579 $ 0.0d0 ),
580 $ x21(i,i), ldx21 )
581 CALL zaxpy( m-p-i+1, dcmplx( -z2*z3*z4*sin(phi(i-1)),
582 $ 0.0d0 ), x22(i-1,i), ldx22, x21(i,i), ldx21 )
583 END IF
584*
585 theta(i) = atan2( dznrm2( m-p-i+1, x21(i,i), ldx21 ),
586 $ dznrm2( p-i+1, x11(i,i), ldx11 ) )
587*
588 CALL zlacgv( p-i+1, x11(i,i), ldx11 )
589 CALL zlacgv( m-p-i+1, x21(i,i), ldx21 )
590*
591 CALL zlarfgp( p-i+1, x11(i,i), x11(i,i+1), ldx11,
592 $ taup1(i) )
593 IF ( i .EQ. m-p ) THEN
594 CALL zlarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
595 $ taup2(i) )
596 ELSE
597 CALL zlarfgp( m-p-i+1, x21(i,i), x21(i,i+1), ldx21,
598 $ taup2(i) )
599 END IF
600*
601 CALL zlarf1f( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
602 $ x11(i+1,i), ldx11, work )
603 CALL zlarf1f( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11,
604 $ taup1(i),
605 $ x12(i,i), ldx12, work )
606 CALL zlarf1f( 'R', q-i, m-p-i+1, x21(i,i), ldx21,
607 $ taup2(i), x21(i+1,i), ldx21, work )
608 CALL zlarf1f( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
609 $ taup2(i), x22(i,i), ldx22, work )
610*
611 CALL zlacgv( p-i+1, x11(i,i), ldx11 )
612 CALL zlacgv( m-p-i+1, x21(i,i), ldx21 )
613*
614 IF( i .LT. q ) THEN
615 CALL zscal( q-i, dcmplx( -z1*z3*sin(theta(i)),
616 $ 0.0d0 ),
617 $ x11(i+1,i), 1 )
618 CALL zaxpy( q-i, dcmplx( z2*z3*cos(theta(i)), 0.0d0 ),
619 $ x21(i+1,i), 1, x11(i+1,i), 1 )
620 END IF
621 CALL zscal( m-q-i+1, dcmplx( -z1*z4*sin(theta(i)),
622 $ 0.0d0 ),
623 $ x12(i,i), 1 )
624 CALL zaxpy( m-q-i+1, dcmplx( z2*z4*cos(theta(i)),
625 $ 0.0d0 ),
626 $ x22(i,i), 1, x12(i,i), 1 )
627*
628 IF( i .LT. q )
629 $ phi(i) = atan2( dznrm2( q-i, x11(i+1,i), 1 ),
630 $ dznrm2( m-q-i+1, x12(i,i), 1 ) )
631*
632 IF( i .LT. q ) THEN
633 CALL zlarfgp( q-i, x11(i+1,i), x11(i+2,i), 1,
634 $ tauq1(i) )
635 END IF
636 CALL zlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1,
637 $ tauq2(i) )
638*
639 IF( i .LT. q ) THEN
640 CALL zlarf1f( 'L', q-i, p-i, x11(i+1,i), 1,
641 $ conjg(tauq1(i)), x11(i+1,i+1), ldx11,
642 $ work )
643 CALL zlarf1f( 'L', q-i, m-p-i, x11(i+1,i), 1,
644 $ conjg(tauq1(i)), x21(i+1,i+1), ldx21,
645 $ work )
646 END IF
647 CALL zlarf1f( 'L', m-q-i+1, p-i, x12(i,i), 1,
648 $ conjg(tauq2(i)),
649 $ x12(i,i+1), ldx12, work )
650
651 IF ( m-p .GT. i ) THEN
652 CALL zlarf1f( 'L', m-q-i+1, m-p-i, x12(i,i), 1,
653 $ conjg(tauq2(i)), x22(i,i+1), ldx22,
654 $ work )
655 END IF
656*
657 END DO
658*
659* Reduce columns Q + 1, ..., P of X12, X22
660*
661 DO i = q + 1, p
662*
663 CALL zscal( m-q-i+1, dcmplx( -z1*z4, 0.0d0 ), x12(i,i),
664 $ 1 )
665 CALL zlarfgp( m-q-i+1, x12(i,i), x12(i+1,i), 1,
666 $ tauq2(i) )
667*
668 IF ( p .GT. i ) THEN
669 CALL zlarf1f( 'L', m-q-i+1, p-i, x12(i,i), 1,
670 $ conjg(tauq2(i)), x12(i,i+1), ldx12,
671 $ work )
672 END IF
673 IF( m-p-q .GE. 1 )
674 $ CALL zlarf1f( 'L', m-q-i+1, m-p-q, x12(i,i), 1,
675 $ conjg(tauq2(i)), x22(i,q+1), ldx22,
676 $ work )
677*
678 END DO
679*
680* Reduce columns P + 1, ..., M - Q of X12, X22
681*
682 DO i = 1, m - p - q
683*
684 CALL zscal( m-p-q-i+1, dcmplx( z2*z4, 0.0d0 ),
685 $ x22(p+i,q+i), 1 )
686 CALL zlarfgp( m-p-q-i+1, x22(p+i,q+i), x22(p+i+1,q+i), 1,
687 $ tauq2(p+i) )
688 IF ( m-p-q .NE. i ) THEN
689 CALL zlarf1f( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i),
690 $ 1, conjg(tauq2(p+i)), x22(p+i,q+i+1),
691 $ ldx22, work )
692 END IF
693*
694 END DO
695*
696 END IF
697*
698 RETURN
699*
700* End of ZUNBDB
701*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
zaxpy
subroutine zaxpy(n, za, zx, incx, zy, incy)
ZAXPY
Definition zaxpy.f:88
zlacgv
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
Definition zlacgv.f:72
zlarf1f
subroutine zlarf1f(side, m, n, v, incv, tau, c, ldc, work)
ZLARF1F applies an elementary reflector to a general rectangular
Definition zlarf1f.f:157
zlarfgp
subroutine zlarfgp(n, alpha, x, incx, tau)
ZLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition zlarfgp.f:102
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dznrm2
real(wp) function dznrm2(n, x, incx)
DZNRM2
Definition dznrm2.f90:90
zscal
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: