266 SUBROUTINE claqr2( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
267 $ IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT,
268 $ NV, WV, LDWV, WORK, LWORK )
275 INTEGER IHIZ, ILOZ, KBOT, KTOP, LDH, LDT, LDV, LDWV,
276 $ LDZ, LWORK, N, ND, NH, NS, NV, NW
280 COMPLEX H( LDH, * ), SH( * ), T( LDT, * ), V( LDV, * ),
281 $ WORK( * ), WV( LDWV, * ), Z( LDZ, * )
288 PARAMETER ( ZERO = ( 0.0e0, 0.0e0 ),
289 $ one = ( 1.0e0, 0.0e0 ) )
291 PARAMETER ( RZERO = 0.0e0, rone = 1.0e0 )
294 COMPLEX BETA, CDUM, S, TAU
295 REAL FOO, SAFMAX, SAFMIN, SMLNUM, ULP
296 INTEGER I, IFST, ILST, INFO, INFQR, J, JW, KCOL, KLN,
297 $ knt, krow, kwtop, ltop, lwk1, lwk2, lwkopt
308 INTRINSIC abs, aimag, cmplx, conjg, int, max, min, real
314 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
320 jw = min( nw, kbot-ktop+1 )
327 CALL cgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
328 lwk1 = int( work( 1 ) )
332 CALL cunmhr(
'R',
'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
334 lwk2 = int( work( 1 ) )
338 lwkopt = jw + max( lwk1, lwk2 )
343 IF( lwork.EQ.-1 )
THEN
344 work( 1 ) = cmplx( lwkopt, 0 )
361 safmin = slamch(
'SAFE MINIMUM' )
362 safmax = rone / safmin
363 CALL slabad( safmin, safmax )
364 ulp = slamch(
'PRECISION' )
365 smlnum = safmin*( real( n ) / ulp )
369 jw = min( nw, kbot-ktop+1 )
370 kwtop = kbot - jw + 1
371 IF( kwtop.EQ.ktop )
THEN
374 s = h( kwtop, kwtop-1 )
377 IF( kbot.EQ.kwtop )
THEN
381 sh( kwtop ) = h( kwtop, kwtop )
384 IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
389 $ h( kwtop, kwtop-1 ) = zero
401 CALL clacpy(
'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
402 CALL ccopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
404 CALL claset(
'A', jw, jw, zero, one, v, ldv )
405 CALL clahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
406 $ jw, v, ldv, infqr )
412 DO 10 knt = infqr + 1, jw
416 foo = cabs1( t( ns, ns ) )
419 IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
431 CALL ctrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
446 DO 30 i = infqr + 1, ns
449 IF( cabs1( t( j, j ) ).GT.cabs1( t( ifst, ifst ) ) )
454 $
CALL ctrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
460 DO 40 i = infqr + 1, jw
461 sh( kwtop+i-1 ) = t( i, i )
465 IF( ns.LT.jw .OR. s.EQ.zero )
THEN
466 IF( ns.GT.1 .AND. s.NE.zero )
THEN
470 CALL ccopy( ns, v, ldv, work, 1 )
472 work( i ) = conjg( work( i ) )
475 CALL clarfg( ns, beta, work( 2 ), 1, tau )
478 CALL claset(
'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
480 CALL clarf(
'L', ns, jw, work, 1, conjg( tau ), t, ldt,
482 CALL clarf(
'R', ns, ns, work, 1, tau, t, ldt,
484 CALL clarf(
'R', jw, ns, work, 1, tau, v, ldv,
487 CALL cgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
494 $ h( kwtop, kwtop-1 ) = s*conjg( v( 1, 1 ) )
495 CALL clacpy(
'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
496 CALL ccopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
502 IF( ns.GT.1 .AND. s.NE.zero )
503 $
CALL cunmhr(
'R',
'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
504 $ work( jw+1 ), lwork-jw, info )
513 DO 60 krow = ltop, kwtop - 1, nv
514 kln = min( nv, kwtop-krow )
515 CALL cgemm(
'N',
'N', kln, jw, jw, one, h( krow, kwtop ),
516 $ ldh, v, ldv, zero, wv, ldwv )
517 CALL clacpy(
'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
523 DO 70 kcol = kbot + 1, n, nh
524 kln = min( nh, n-kcol+1 )
525 CALL cgemm(
'C',
'N', jw, kln, jw, one, v, ldv,
526 $ h( kwtop, kcol ), ldh, zero, t, ldt )
527 CALL clacpy(
'A', jw, kln, t, ldt, h( kwtop, kcol ),
535 DO 80 krow = iloz, ihiz, nv
536 kln = min( nv, ihiz-krow+1 )
537 CALL cgemm(
'N',
'N', kln, jw, jw, one, z( krow, kwtop ),
538 $ ldz, v, ldv, zero, wv, ldwv )
539 CALL clacpy(
'A', kln, jw, wv, ldwv, z( krow, kwtop ),
559 work( 1 ) = cmplx( lwkopt, 0 )
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.
subroutine claqr2(WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ, IHIZ, Z, LDZ, NS, ND, SH, V, LDV, NH, T, LDT, NV, WV, LDWV, WORK, LWORK)
CLAQR2 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
subroutine clahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
CLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMHR
subroutine ctrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
CTREXC