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