266 SUBROUTINE zlaqr3( 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 )
276 INTEGER ihiz, iloz, kbot, ktop, ldh, ldt, ldv, ldwv,
277 $ ldz, lwork, n, nd, nh, ns, nv, nw
281 COMPLEX*16 h( ldh, * ), sh( * ), t( ldt, * ), v( ldv, * ),
282 $ work( * ), wv( ldwv, * ), z( ldz, * )
289 parameter( zero = ( 0.0d0, 0.0d0 ),
290 $ one = ( 1.0d0, 0.0d0 ) )
291 DOUBLE PRECISION rzero, rone
292 parameter( rzero = 0.0d0, rone = 1.0d0 )
295 COMPLEX*16 beta, cdum, s, tau
296 DOUBLE PRECISION foo, safmax, safmin, smlnum, ulp
297 INTEGER i, ifst, ilst, info, infqr,
j, jw, kcol, kln,
298 $ knt, krow, kwtop, ltop, lwk1, lwk2, lwk3,
311 INTRINSIC abs, dble, dcmplx, dconjg, dimag, int, max, min
314 DOUBLE PRECISION cabs1
317 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
323 jw = min( nw, kbot-ktop+1 )
330 CALL
zgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
331 lwk1 = int( work( 1 ) )
335 CALL
zunmhr(
'R',
'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
337 lwk2 = int( work( 1 ) )
341 CALL
zlaqr4( .true., .true., jw, 1, jw, t, ldt, sh, 1, jw, v,
342 $ ldv, work, -1, infqr )
343 lwk3 = int( work( 1 ) )
347 lwkopt = max( jw+max( lwk1, lwk2 ), lwk3 )
352 IF( lwork.EQ.-1 )
THEN
353 work( 1 ) = dcmplx( lwkopt, 0 )
370 safmin =
dlamch(
'SAFE MINIMUM' )
371 safmax = rone / safmin
372 CALL
dlabad( safmin, safmax )
373 ulp =
dlamch(
'PRECISION' )
374 smlnum = safmin*( dble( n ) / ulp )
378 jw = min( nw, kbot-ktop+1 )
379 kwtop = kbot - jw + 1
380 IF( kwtop.EQ.ktop )
THEN
383 s = h( kwtop, kwtop-1 )
386 IF( kbot.EQ.kwtop )
THEN
390 sh( kwtop ) = h( kwtop, kwtop )
393 IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
398 $ h( kwtop, kwtop-1 ) = zero
410 CALL
zlacpy(
'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
411 CALL
zcopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
413 CALL
zlaset(
'A', jw, jw, zero, one, v, ldv )
414 nmin =
ilaenv( 12,
'ZLAQR3',
'SV', jw, 1, jw, lwork )
415 IF( jw.GT.nmin )
THEN
416 CALL
zlaqr4( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
417 $ jw, v, ldv, work, lwork, infqr )
419 CALL
zlahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
420 $ jw, v, ldv, infqr )
427 DO 10 knt = infqr + 1, jw
431 foo = cabs1( t( ns, ns ) )
434 IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
446 CALL
ztrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
461 DO 30 i = infqr + 1, ns
464 IF( cabs1( t(
j,
j ) ).GT.cabs1( t( ifst, ifst ) ) )
469 $ CALL
ztrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
475 DO 40 i = infqr + 1, jw
476 sh( kwtop+i-1 ) = t( i, i )
480 IF( ns.LT.jw .OR. s.EQ.zero )
THEN
481 IF( ns.GT.1 .AND. s.NE.zero )
THEN
485 CALL
zcopy( ns, v, ldv, work, 1 )
487 work( i ) = dconjg( work( i ) )
490 CALL
zlarfg( ns, beta, work( 2 ), 1, tau )
493 CALL
zlaset(
'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
495 CALL
zlarf(
'L', ns, jw, work, 1, dconjg( tau ), t, ldt,
497 CALL
zlarf(
'R', ns, ns, work, 1, tau, t, ldt,
499 CALL
zlarf(
'R', jw, ns, work, 1, tau, v, ldv,
502 CALL
zgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
509 $ h( kwtop, kwtop-1 ) = s*dconjg( v( 1, 1 ) )
510 CALL
zlacpy(
'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
511 CALL
zcopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
517 IF( ns.GT.1 .AND. s.NE.zero )
518 $ CALL
zunmhr(
'R',
'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
519 $ work( jw+1 ), lwork-jw, info )
528 DO 60 krow = ltop, kwtop - 1, nv
529 kln = min( nv, kwtop-krow )
530 CALL
zgemm(
'N',
'N', kln, jw, jw, one, h( krow, kwtop ),
531 $ ldh, v, ldv, zero, wv, ldwv )
532 CALL
zlacpy(
'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
538 DO 70 kcol = kbot + 1, n, nh
539 kln = min( nh, n-kcol+1 )
540 CALL
zgemm(
'C',
'N', jw, kln, jw, one, v, ldv,
541 $ h( kwtop, kcol ), ldh, zero, t, ldt )
542 CALL
zlacpy(
'A', jw, kln, t, ldt, h( kwtop, kcol ),
550 DO 80 krow = iloz, ihiz, nv
551 kln = min( nv, ihiz-krow+1 )
552 CALL
zgemm(
'N',
'N', kln, jw, jw, one, z( krow, kwtop ),
553 $ ldz, v, ldv, zero, wv, ldwv )
554 CALL
zlacpy(
'A', kln, jw, wv, ldwv, z( krow, kwtop ),
574 work( 1 ) = dcmplx( lwkopt, 0 )
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
ZLARF applies an elementary reflector to a general rectangular matrix.
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine zcopy(N, ZX, INCX, ZY, INCY)
ZCOPY
subroutine zlahqr(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, INFO)
ZLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix, using the double-shift/single-shift QR algorithm.
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine ztrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
ZTREXC
subroutine zunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMHR
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine zlarfg(N, ALPHA, X, INCX, TAU)
ZLARFG generates an elementary reflector (Householder matrix).
subroutine zlaqr3(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)
ZLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
subroutine zlaqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
ZLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...