265 SUBROUTINE claqr3( WANTT, WANTZ, N, KTOP, KBOT, NW, H, LDH, ILOZ,
266 $ ihiz, z, ldz, ns, nd, sh, v, ldv, nh, t, ldt,
267 $ 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, lwk3,
310 INTRINSIC abs, aimag, cmplx, conjg, int, max, min, real
316 cabs1( cdum ) = abs(
REAL( CDUM ) ) + abs( aimag( cdum ) )
322 jw = min( nw, kbot-ktop+1 )
329 CALL
cgehrd( jw, 1, jw-1, t, ldt, work, work, -1, info )
330 lwk1 = int( work( 1 ) )
334 CALL
cunmhr(
'R',
'N', jw, jw, 1, jw-1, t, ldt, work, v, ldv,
336 lwk2 = int( work( 1 ) )
340 CALL
claqr4( .true., .true., jw, 1, jw, t, ldt, sh, 1, jw, v,
341 $ ldv, work, -1, infqr )
342 lwk3 = int( work( 1 ) )
346 lwkopt = max( jw+max( lwk1, lwk2 ), lwk3 )
351 IF( lwork.EQ.-1 )
THEN
352 work( 1 ) = cmplx( lwkopt, 0 )
369 safmin =
slamch(
'SAFE MINIMUM' )
370 safmax = rone / safmin
371 CALL
slabad( safmin, safmax )
372 ulp =
slamch(
'PRECISION' )
373 smlnum = safmin*(
REAL( N ) / ulp )
377 jw = min( nw, kbot-ktop+1 )
378 kwtop = kbot - jw + 1
379 IF( kwtop.EQ.ktop )
THEN
382 s = h( kwtop, kwtop-1 )
385 IF( kbot.EQ.kwtop )
THEN
389 sh( kwtop ) = h( kwtop, kwtop )
392 IF( cabs1( s ).LE.max( smlnum, ulp*cabs1( h( kwtop,
397 $ h( kwtop, kwtop-1 ) = zero
409 CALL
clacpy(
'U', jw, jw, h( kwtop, kwtop ), ldh, t, ldt )
410 CALL
ccopy( jw-1, h( kwtop+1, kwtop ), ldh+1, t( 2, 1 ), ldt+1 )
412 CALL
claset(
'A', jw, jw, zero, one, v, ldv )
413 nmin =
ilaenv( 12,
'CLAQR3',
'SV', jw, 1, jw, lwork )
414 IF( jw.GT.nmin )
THEN
415 CALL
claqr4( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
416 $ jw, v, ldv, work, lwork, infqr )
418 CALL
clahqr( .true., .true., jw, 1, jw, t, ldt, sh( kwtop ), 1,
419 $ jw, v, ldv, infqr )
426 DO 10 knt = infqr + 1, jw
430 foo = cabs1( t( ns, ns ) )
433 IF( cabs1( s )*cabs1( v( 1, ns ) ).LE.max( smlnum, ulp*foo ) )
445 CALL
ctrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
460 DO 30 i = infqr + 1, ns
463 IF( cabs1( t(
j,
j ) ).GT.cabs1( t( ifst, ifst ) ) )
468 $ CALL
ctrexc(
'V', jw, t, ldt, v, ldv, ifst, ilst, info )
474 DO 40 i = infqr + 1, jw
475 sh( kwtop+i-1 ) = t( i, i )
479 IF( ns.LT.jw .OR. s.EQ.zero )
THEN
480 IF( ns.GT.1 .AND. s.NE.zero )
THEN
484 CALL
ccopy( ns, v, ldv, work, 1 )
486 work( i ) = conjg( work( i ) )
489 CALL
clarfg( ns, beta, work( 2 ), 1, tau )
492 CALL
claset(
'L', jw-2, jw-2, zero, zero, t( 3, 1 ), ldt )
494 CALL
clarf(
'L', ns, jw, work, 1, conjg( tau ), t, ldt,
496 CALL
clarf(
'R', ns, ns, work, 1, tau, t, ldt,
498 CALL
clarf(
'R', jw, ns, work, 1, tau, v, ldv,
501 CALL
cgehrd( jw, 1, ns, t, ldt, work, work( jw+1 ),
508 $ h( kwtop, kwtop-1 ) = s*conjg( v( 1, 1 ) )
509 CALL
clacpy(
'U', jw, jw, t, ldt, h( kwtop, kwtop ), ldh )
510 CALL
ccopy( jw-1, t( 2, 1 ), ldt+1, h( kwtop+1, kwtop ),
516 IF( ns.GT.1 .AND. s.NE.zero )
517 $ CALL
cunmhr(
'R',
'N', jw, ns, 1, ns, t, ldt, work, v, ldv,
518 $ work( jw+1 ), lwork-jw, info )
527 DO 60 krow = ltop, kwtop - 1, nv
528 kln = min( nv, kwtop-krow )
529 CALL
cgemm(
'N',
'N', kln, jw, jw, one, h( krow, kwtop ),
530 $ ldh, v, ldv, zero, wv, ldwv )
531 CALL
clacpy(
'A', kln, jw, wv, ldwv, h( krow, kwtop ), ldh )
537 DO 70 kcol = kbot + 1, n, nh
538 kln = min( nh, n-kcol+1 )
539 CALL
cgemm(
'C',
'N', jw, kln, jw, one, v, ldv,
540 $ h( kwtop, kcol ), ldh, zero, t, ldt )
541 CALL
clacpy(
'A', jw, kln, t, ldt, h( kwtop, kcol ),
549 DO 80 krow = iloz, ihiz, nv
550 kln = min( nv, ihiz-krow+1 )
551 CALL
cgemm(
'N',
'N', kln, jw, jw, one, z( krow, kwtop ),
552 $ ldz, v, ldv, zero, wv, ldwv )
553 CALL
clacpy(
'A', kln, jw, wv, ldwv, z( krow, kwtop ),
573 work( 1 ) = cmplx( lwkopt, 0 )
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 ctrexc(COMPQ, N, T, LDT, Q, LDQ, IFST, ILST, INFO)
CTREXC
subroutine claqr4(WANTT, WANTZ, N, ILO, IHI, H, LDH, W, ILOZ, IHIZ, Z, LDZ, WORK, LWORK, INFO)
CLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
subroutine cunmhr(SIDE, TRANS, M, N, ILO, IHI, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMHR
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(CMACH)
SLAMCH
subroutine claqr3(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)
CLAQR3 performs the unitary similarity transformation of a Hessenberg matrix to detect and deflate fu...
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine cgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
CGEHRD
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine ccopy(N, CX, INCX, CY, INCY)
CCOPY
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
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, using the double-shift/single-shift QR algorithm.
subroutine clarf(SIDE, M, N, V, INCV, TAU, C, LDC, WORK)
CLARF applies an elementary reflector to a general rectangular matrix.