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

◆ cunbdb()

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.
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 284 of file cunbdb.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 REAL PHI( * ), THETA( * )
299 COMPLEX TAUP1( * ), TAUP2( * ), TAUQ1( * ), TAUQ2( * ),
300 $ WORK( * ), X11( LDX11, * ), X12( LDX12, * ),
301 $ X21( LDX21, * ), X22( LDX22, * )
302* ..
303*
304* ====================================================================
305*
306* .. Parameters ..
307 REAL REALONE
308 parameter( realone = 1.0e0 )
309 COMPLEX ONE
310 parameter( one = (1.0e0,0.0e0) )
311* ..
312* .. Local Scalars ..
313 LOGICAL COLMAJOR, LQUERY
314 INTEGER I, LWORKMIN, LWORKOPT
315 REAL Z1, Z2, Z3, Z4
316* ..
317* .. External Subroutines ..
318 EXTERNAL caxpy, clarf, clarfgp, cscal, xerbla
319 EXTERNAL clacgv
320*
321* ..
322* .. External Functions ..
323 REAL SCNRM2
324 LOGICAL LSAME
325 EXTERNAL scnrm2, lsame
326* ..
327* .. Intrinsic Functions
328 INTRINSIC atan2, cos, max, min, sin
329 INTRINSIC cmplx, conjg
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 cscal( p-i+1, cmplx( z1, 0.0e0 ), x11(i,i), 1 )
402 ELSE
403 CALL cscal( p-i+1, cmplx( z1*cos(phi(i-1)), 0.0e0 ),
404 $ x11(i,i), 1 )
405 CALL caxpy( p-i+1, cmplx( -z1*z3*z4*sin(phi(i-1)),
406 $ 0.0e0 ), x12(i,i-1), 1, x11(i,i), 1 )
407 END IF
408 IF( i .EQ. 1 ) THEN
409 CALL cscal( m-p-i+1, cmplx( z2, 0.0e0 ), x21(i,i), 1 )
410 ELSE
411 CALL cscal( m-p-i+1, cmplx( z2*cos(phi(i-1)), 0.0e0 ),
412 $ x21(i,i), 1 )
413 CALL caxpy( m-p-i+1, cmplx( -z2*z3*z4*sin(phi(i-1)),
414 $ 0.0e0 ), x22(i,i-1), 1, x21(i,i), 1 )
415 END IF
416*
417 theta(i) = atan2( scnrm2( m-p-i+1, x21(i,i), 1 ),
418 $ scnrm2( p-i+1, x11(i,i), 1 ) )
419*
420 IF( p .GT. i ) THEN
421 CALL clarfgp( p-i+1, x11(i,i), x11(i+1,i), 1, taup1(i) )
422 ELSE IF ( p .EQ. i ) THEN
423 CALL clarfgp( 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 clarfgp( 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 clarfgp( 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 clarf( 'L', p-i+1, q-i, x11(i,i), 1,
437 $ conjg(taup1(i)), x11(i,i+1), ldx11, work )
438 CALL clarf( 'L', m-p-i+1, q-i, x21(i,i), 1,
439 $ conjg(taup2(i)), x21(i,i+1), ldx21, work )
440 END IF
441 IF ( m-q+1 .GT. i ) THEN
442 CALL clarf( 'L', p-i+1, m-q-i+1, x11(i,i), 1,
443 $ conjg(taup1(i)), x12(i,i), ldx12, work )
444 CALL clarf( 'L', m-p-i+1, m-q-i+1, x21(i,i), 1,
445 $ conjg(taup2(i)), x22(i,i), ldx22, work )
446 END IF
447*
448 IF( i .LT. q ) THEN
449 CALL cscal( q-i, cmplx( -z1*z3*sin(theta(i)), 0.0e0 ),
450 $ x11(i,i+1), ldx11 )
451 CALL caxpy( q-i, cmplx( z2*z3*cos(theta(i)), 0.0e0 ),
452 $ x21(i,i+1), ldx21, x11(i,i+1), ldx11 )
453 END IF
454 CALL cscal( m-q-i+1, cmplx( -z1*z4*sin(theta(i)), 0.0e0 ),
455 $ x12(i,i), ldx12 )
456 CALL caxpy( m-q-i+1, cmplx( z2*z4*cos(theta(i)), 0.0e0 ),
457 $ x22(i,i), ldx22, x12(i,i), ldx12 )
458*
459 IF( i .LT. q )
460 $ phi(i) = atan2( scnrm2( q-i, x11(i,i+1), ldx11 ),
461 $ scnrm2( m-q-i+1, x12(i,i), ldx12 ) )
462*
463 IF( i .LT. q ) THEN
464 CALL clacgv( q-i, x11(i,i+1), ldx11 )
465 IF ( i .EQ. q-1 ) THEN
466 CALL clarfgp( q-i, x11(i,i+1), x11(i,i+1), ldx11,
467 $ tauq1(i) )
468 ELSE
469 CALL clarfgp( 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 clacgv( m-q-i+1, x12(i,i), ldx12 )
476 IF ( m-q .EQ. i ) THEN
477 CALL clarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
478 $ tauq2(i) )
479 ELSE
480 CALL clarfgp( 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 clarf( 'R', p-i, q-i, x11(i,i+1), ldx11, tauq1(i),
488 $ x11(i+1,i+1), ldx11, work )
489 CALL clarf( '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 clarf( '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 clarf( '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 clacgv( q-i, x11(i,i+1), ldx11 )
503 CALL clacgv( 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 cscal( m-q-i+1, cmplx( -z1*z4, 0.0e0 ), x12(i,i),
512 $ ldx12 )
513 CALL clacgv( m-q-i+1, x12(i,i), ldx12 )
514 IF ( i .GE. m-q ) THEN
515 CALL clarfgp( m-q-i+1, x12(i,i), x12(i,i), ldx12,
516 $ tauq2(i) )
517 ELSE
518 CALL clarfgp( 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 clarf( '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 clarf( 'R', m-p-q, m-q-i+1, x12(i,i), ldx12,
529 $ tauq2(i), x22(q+1,i), ldx22, work )
530*
531 CALL clacgv( 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 cscal( m-p-q-i+1, cmplx( z2*z4, 0.0e0 ),
540 $ x22(q+i,p+i), ldx22 )
541 CALL clacgv( m-p-q-i+1, x22(q+i,p+i), ldx22 )
542 CALL clarfgp( 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 clarf( '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 clacgv( 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 cscal( p-i+1, cmplx( z1, 0.0e0 ), x11(i,i),
560 $ ldx11 )
561 ELSE
562 CALL cscal( p-i+1, cmplx( z1*cos(phi(i-1)), 0.0e0 ),
563 $ x11(i,i), ldx11 )
564 CALL caxpy( p-i+1, cmplx( -z1*z3*z4*sin(phi(i-1)),
565 $ 0.0e0 ), x12(i-1,i), ldx12, x11(i,i), ldx11 )
566 END IF
567 IF( i .EQ. 1 ) THEN
568 CALL cscal( m-p-i+1, cmplx( z2, 0.0e0 ), x21(i,i),
569 $ ldx21 )
570 ELSE
571 CALL cscal( m-p-i+1, cmplx( z2*cos(phi(i-1)), 0.0e0 ),
572 $ x21(i,i), ldx21 )
573 CALL caxpy( m-p-i+1, cmplx( -z2*z3*z4*sin(phi(i-1)),
574 $ 0.0e0 ), x22(i-1,i), ldx22, x21(i,i), ldx21 )
575 END IF
576*
577 theta(i) = atan2( scnrm2( m-p-i+1, x21(i,i), ldx21 ),
578 $ scnrm2( p-i+1, x11(i,i), ldx11 ) )
579*
580 CALL clacgv( p-i+1, x11(i,i), ldx11 )
581 CALL clacgv( m-p-i+1, x21(i,i), ldx21 )
582*
583 CALL clarfgp( 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 clarfgp( m-p-i+1, x21(i,i), x21(i,i), ldx21,
587 $ taup2(i) )
588 ELSE
589 CALL clarfgp( 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 clarf( 'R', q-i, p-i+1, x11(i,i), ldx11, taup1(i),
595 $ x11(i+1,i), ldx11, work )
596 CALL clarf( 'R', m-q-i+1, p-i+1, x11(i,i), ldx11, taup1(i),
597 $ x12(i,i), ldx12, work )
598 CALL clarf( 'R', q-i, m-p-i+1, x21(i,i), ldx21, taup2(i),
599 $ x21(i+1,i), ldx21, work )
600 CALL clarf( 'R', m-q-i+1, m-p-i+1, x21(i,i), ldx21,
601 $ taup2(i), x22(i,i), ldx22, work )
602*
603 CALL clacgv( p-i+1, x11(i,i), ldx11 )
604 CALL clacgv( m-p-i+1, x21(i,i), ldx21 )
605*
606 IF( i .LT. q ) THEN
607 CALL cscal( q-i, cmplx( -z1*z3*sin(theta(i)), 0.0e0 ),
608 $ x11(i+1,i), 1 )
609 CALL caxpy( q-i, cmplx( z2*z3*cos(theta(i)), 0.0e0 ),
610 $ x21(i+1,i), 1, x11(i+1,i), 1 )
611 END IF
612 CALL cscal( m-q-i+1, cmplx( -z1*z4*sin(theta(i)), 0.0e0 ),
613 $ x12(i,i), 1 )
614 CALL caxpy( m-q-i+1, cmplx( z2*z4*cos(theta(i)), 0.0e0 ),
615 $ x22(i,i), 1, x12(i,i), 1 )
616*
617 IF( i .LT. q )
618 $ phi(i) = atan2( scnrm2( q-i, x11(i+1,i), 1 ),
619 $ scnrm2( m-q-i+1, x12(i,i), 1 ) )
620*
621 IF( i .LT. q ) THEN
622 CALL clarfgp( q-i, x11(i+1,i), x11(i+2,i), 1, tauq1(i) )
623 x11(i+1,i) = one
624 END IF
625 CALL clarfgp( 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 clarf( 'L', q-i, p-i, x11(i+1,i), 1,
630 $ conjg(tauq1(i)), x11(i+1,i+1), ldx11, work )
631 CALL clarf( 'L', q-i, m-p-i, x11(i+1,i), 1,
632 $ conjg(tauq1(i)), x21(i+1,i+1), ldx21, work )
633 END IF
634 CALL clarf( 'L', m-q-i+1, p-i, x12(i,i), 1, conjg(tauq2(i)),
635 $ x12(i,i+1), ldx12, work )
636
637 IF ( m-p .GT. i ) THEN
638 CALL clarf( 'L', m-q-i+1, m-p-i, x12(i,i), 1,
639 $ conjg(tauq2(i)), x22(i,i+1), ldx22, work )
640 END IF
641 END DO
642*
643* Reduce columns Q + 1, ..., P of X12, X22
644*
645 DO i = q + 1, p
646*
647 CALL cscal( m-q-i+1, cmplx( -z1*z4, 0.0e0 ), x12(i,i), 1 )
648 CALL clarfgp( 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 clarf( 'L', m-q-i+1, p-i, x12(i,i), 1,
653 $ conjg(tauq2(i)), x12(i,i+1), ldx12, work )
654 END IF
655 IF( m-p-q .GE. 1 )
656 $ CALL clarf( 'L', m-q-i+1, m-p-q, x12(i,i), 1,
657 $ conjg(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 cscal( m-p-q-i+1, cmplx( z2*z4, 0.0e0 ),
666 $ x22(p+i,q+i), 1 )
667 CALL clarfgp( 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 IF ( m-p-q .NE. i ) THEN
671 CALL clarf( 'L', m-p-q-i+1, m-p-q-i, x22(p+i,q+i), 1,
672 $ conjg(tauq2(p+i)), x22(p+i,q+i+1), ldx22,
673 $ work )
674 END IF
675 END DO
676*
677 END IF
678*
679 RETURN
680*
681* End of CUNBDB
682*
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine caxpy(N, CA, CX, INCX, CY, INCY)
CAXPY
Definition: caxpy.f:88
subroutine cscal(N, CA, CX, INCX)
CSCAL
Definition: cscal.f:78
subroutine clarfgp(N, ALPHA, X, INCX, TAU)
CLARFGP generates an elementary reflector (Householder matrix) with non-negative beta.
Definition: clarfgp.f:104
subroutine clacgv(N, X, INCX)
CLACGV conjugates a complex vector.
Definition: clacgv.f:74
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
Definition: clarf.f:128
real(wp) function scnrm2(n, x, incx)
SCNRM2
Definition: scnrm2.f90:90
Here is the call graph for this function:
Here is the caller graph for this function: