SUBROUTINE ASYMTX( AA, EVEC, EVAL, M, IA, IEVEC, IER, WKD, AAD, 1,3
&                   EVECD, EVALD )

c    =======  D O U B L E    P R E C I S I O N    V E R S I O N  ======

c       Solves eigenfunction problem for real asymmetric matrix
c       for which it is known a priori that the eigenvalues are real.

c       This is an adaptation of a subroutine EIGRF in the IMSL
c       library to use real instead of complex arithmetic, accounting
c       for the known fact that the eigenvalues and eigenvectors in
c       the discrete ordinate solution are real.  Other changes include
c       putting all the called subroutines in-line, deleting the
c       performance index calculation, updating many DO-loops
c       to Fortran77, and in calculating the machine precision
c       TOL instead of specifying it in a data statement.

c       EIGRF is based primarily on EISPACK routines.  The matrix is
c       first balanced using the Parlett-Reinsch algorithm.  Then
c       the Martin-Wilkinson algorithm is applied.

c       There is a statement 'J  = WKD( I )' that converts a double
c       precision variable to an integer variable, that seems dangerous
c       to us in principle, but seems to work fine in practice.

c       References:
c          Dongarra, J. and C. Moler, EISPACK -- A Package for Solving
c             Matrix Eigenvalue Problems, in Cowell, ed., 1984:
c             Sources and Development of Mathematical Software,
c             Prentice-Hall, Englewood Cliffs, NJ
c         Parlett and Reinsch, 1969: Balancing a Matrix for Calculation
c             of Eigenvalues and Eigenvectors, Num. Math. 13, 293-304
c         Wilkinson, J., 1965: The Algebraic Eigenvalue Problem,
c             Clarendon Press, Oxford
c
c
c   I N P U T    V A R I A B L E S:
c
c       AA    :  input asymmetric matrix, destroyed after solved
c
c        M    :  order of  AA
c
c       IA    :  first dimension of  AA
c
c    IEVEC    :  first dimension of  EVEC
c
c
c   O U T P U T    V A R I A B L E S:
c
c       EVEC  :  (unnormalized) eigenvectors of  AA
c                ( column J corresponds to EVAL(J) )
c
c       EVAL  :  (unordered) eigenvalues of AA ( dimension at least M )
c
c       IER   :  if .NE. 0, signals that EVAL(IER) failed to converge;
c                in that case eigenvalues IER+1,IER+2,...,M  are
c                correct but eigenvalues 1,...,IER are set to zero.
c
c
c   S C R A T C H   V A R I A B L E S:
c
c       WKD   :  work area ( dimension at least 2*M )
c       AAD   :  double precision stand-in for AA
c       EVECD :  double precision stand-in for EVEC
c       EVALD :  double precision stand-in for EVAL
c
c   Called by- SOLEIG
c   Calls- D1MACH, ERRMSG
c +-------------------------------------------------------------------+

c     .. Scalar Arguments ..

INTEGER   IA, IER, IEVEC, M
c     ..
c     .. Array Arguments ..

REAL      AA( IA, M ), EVAL( M ), EVEC( IEVEC, M )
DOUBLE PRECISION AAD( IA, M ), EVALD( M ), EVECD( IA, M ),
&                 WKD( * )
c     ..
c     .. Local Scalars ..

LOGICAL   NOCONV, NOTLAS
INTEGER   I, II, IN, J, K, KA, KKK, L, LB, LLL, N, N1, N2
DOUBLE PRECISION C1, C2, C3, C4, C5, C6, COL, DISCRI, F, G, H,
&                 ONE, P, Q, R, REPL, RNORM, ROW, S, SCALE, SGN, T,
&                 TOL, UU, VV, W, X, Y, Z, ZERO
c     ..
c     .. External Functions ..

DOUBLE PRECISION D1MACH
EXTERNAL  D1MACH
c     ..
c     .. External Subroutines ..

EXTERNAL  ERRMSG
c     ..
c     .. Intrinsic Functions ..

INTRINSIC ABS, MIN, SIGN, SQRT
c     ..
DATA      C1 / 0.4375D0 / , C2 / 0.5D0 / , C3 / 0.75D0 / ,
&          C4 / 0.95D0 / , C5 / 16.D0 / , C6 / 256.D0 / ,
&          ZERO / 0.D0 / , ONE / 1.D0 /

IER  = 0
TOL  = D1MACH( 4 )

IF( M.LT.1 .OR. IA.LT.M .OR. IEVEC.LT.M )
&    CALL ERRMSG( 'ASYMTX--bad input variable(s)', .TRUE. )

c                           ** Handle 1x1 and 2x2 special cases
IF( M.EQ.1 ) THEN

EVAL( 1 )   = AA( 1,1 )
EVEC( 1,1 ) = 1.0
RETURN

ELSE IF( M.EQ.2 ) THEN

DISCRI = ( AA( 1,1 ) - AA( 2,2 ) )**2 + 4.*AA( 1,2 )*AA( 2,1 )

IF( DISCRI .LT. 0.0 )
&       CALL ERRMSG( 'ASYMTX--complex evals in 2x2 case',.TRUE. )

SGN  = ONE

IF( AA( 1,1 ) .LT. AA( 2,2 ) ) SGN  = - ONE

EVAL( 1 ) = 0.5*( AA( 1,1 ) + AA( 2,2 ) + SGN*SQRT( DISCRI ) )
EVAL( 2 ) = 0.5*( AA( 1,1 ) + AA( 2,2 ) - SGN*SQRT( DISCRI ) )
EVEC( 1,1 ) = 1.0
EVEC( 2,2 ) = 1.0

IF( AA( 1,1 ) .EQ. AA( 2,2 ) .AND.
&       ( AA( 2,1 ).EQ.0.0 .OR. AA( 1,2 ).EQ.0.0 ) ) THEN

RNORM = ABS( AA( 1,1 ) ) + ABS( AA( 1,2 ) ) +
&              ABS( AA( 2,1 ) ) + ABS( AA( 2,2 ) )
W     = TOL * RNORM
EVEC( 2,1 ) =   AA( 2,1 ) / W
EVEC( 1,2 ) = - AA( 1,2 ) / W

ELSE

EVEC( 2,1 ) = AA( 2,1 ) / ( EVAL( 1 ) - AA( 2,2 ) )
EVEC( 1,2 ) = AA( 1,2 ) / ( EVAL( 2 ) - AA( 1,1 ) )

END IF

RETURN

END IF

c                               ** Convert single-prec. matrix to double
DO 20 J = 1, M

DO 10 K = 1, M
AAD( J,K ) = AA( J,K )
10    CONTINUE

20 CONTINUE

c                                ** Initialize output variables
IER  = 0

DO 40 I = 1, M
EVALD( I ) = ZERO

DO 30 J = 1, M
EVECD( I, J ) = ZERO
30    CONTINUE

EVECD( I, I ) = ONE
40 CONTINUE

c                  ** Balance the input matrix and reduce its norm by
c                  ** diagonal similarity transformation stored in WK;
c                  ** then search for rows isolating an eigenvalue
c                  ** and push them down
RNORM  = ZERO
L  = 1
K  = M

50 CONTINUE
KKK  = K

DO 90 J = KKK, 1, -1

ROW  = ZERO

DO 60 I = 1, K

IF( I.NE.J ) ROW  = ROW + ABS( AAD( J,I ) )

60    CONTINUE

IF( ROW.EQ.ZERO ) THEN

WKD( K ) = J

IF( J.NE.K ) THEN

DO 70 I = 1, K
REPL        = AAD( I, J )
AAD( I, K ) = REPL
70          CONTINUE

DO 80 I = L, M
REPL        = AAD( J, I )
AAD( K, I ) = REPL
80          CONTINUE

END IF

K  = K - 1
GO TO  50

END IF

90 CONTINUE
c                                ** Search for columns isolating an
c                                ** eigenvalue and push them left
100 CONTINUE
LLL  = L

DO 140 J = LLL, K

COL  = ZERO

DO 110 I = L, K

IF( I.NE.J ) COL  = COL + ABS( AAD( I,J ) )

110    CONTINUE

IF( COL.EQ.ZERO ) THEN

WKD( L ) = J

IF( J.NE.L ) THEN

DO 120 I = 1, K
REPL        = AAD( I, J )
AAD( I, L ) = REPL
120          CONTINUE

DO 130 I = L, M
REPL        = AAD( J, I )
AAD( L, I ) = REPL
130          CONTINUE

END IF

L  = L + 1
GO TO  100

END IF

140 CONTINUE

c                           ** Balance the submatrix in rows L through K
DO 150 I = L, K
WKD( I ) = ONE
150 CONTINUE

160 CONTINUE
NOCONV = .FALSE.

DO 220 I = L, K

COL  = ZERO
ROW  = ZERO

DO 170 J = L, K

IF( J.NE.I ) THEN

COL  = COL + ABS( AAD( J,I ) )
ROW  = ROW + ABS( AAD( I,J ) )

END IF

170    CONTINUE

F  = ONE
G  = ROW / C5
H  = COL + ROW

180    CONTINUE
IF( COL.LT.G ) THEN

F    = F*C5
COL  = COL*C6
GO TO  180

END IF

G  = ROW*C5

190    CONTINUE
IF( COL.GE.G ) THEN

F    = F / C5
COL  = COL / C6
GO TO  190

END IF
c                                                ** Now balance
IF(
WKD( I ) = WKD( I )*F
NOCONV = .TRUE.

DO 200 J = L, M
200       CONTINUE

DO 210 J = 1, K
210       CONTINUE

END IF

220 CONTINUE

IF( NOCONV ) GO TO  160
c                                   ** Is A already in Hessenberg form?
IF( K-1 .LT. L+1 ) GO TO  370

c                                   ** Transfer A to a Hessenberg form
DO 310 N = L + 1, K - 1

H  = ZERO
WKD( N + M ) = ZERO
SCALE  = ZERO
c                                                 ** Scale column
DO 230 I = N, K
SCALE  = SCALE + ABS( AAD( I,N - 1 ) )
230    CONTINUE

IF( SCALE.NE.ZERO ) THEN

DO 240 I = K, N, -1
WKD( I + M ) = AAD( I, N - 1 ) / SCALE
H  = H + WKD( I + M )**2
240       CONTINUE

G    = - SIGN( SQRT( H ), WKD( N + M ) )
H    = H - WKD( N + M )*G
WKD( N + M ) = WKD( N + M ) - G
c                                            ** Form (I-(U*UT)/H)*A
DO 270 J = N, M

F  = ZERO

DO 250 I = K, N, -1
F  = F + WKD( I + M )*AAD( I, J )
250          CONTINUE

DO 260 I = N, K
AAD( I, J ) = AAD( I, J ) - WKD( I + M )*F / H
260          CONTINUE

270       CONTINUE
c                                    ** Form (I-(U*UT)/H)*A*(I-(U*UT)/H)
DO 300 I = 1, K

F  = ZERO

DO 280 J = K, N, -1
F  = F + WKD( J + M )*AAD( I, J )
280          CONTINUE

DO 290 J = N, K
AAD( I, J ) = AAD( I, J ) - WKD( J + M )*F / H
290          CONTINUE

300       CONTINUE

WKD( N + M ) = SCALE*WKD( N + M )
AAD( N, N - 1 ) = SCALE*G

END IF

310 CONTINUE

DO 360 N = K - 2, L, -1

N1   = N + 1
N2   = N + 2
F  = AAD( N + 1, N )

IF( F.NE.ZERO ) THEN

F  = F*WKD( N + 1 + M )

DO 320 I = N + 2, K
WKD( I + M ) = AAD( I, N )
320       CONTINUE

IF( N + 1.LE.K ) THEN

DO 350 J = 1, M

G  = ZERO

DO 330 I = N + 1, K
G  = G + WKD( I + M )*EVECD( I, J )
330             CONTINUE

G  = G / F

DO 340 I = N + 1, K
EVECD( I, J ) = EVECD( I, J ) + G*WKD( I + M )
340             CONTINUE

350          CONTINUE

END IF

END IF

360 CONTINUE

370 CONTINUE

N  = 1

DO 390 I = 1, M

DO 380 J = N, M
RNORM  = RNORM + ABS( AAD( I,J ) )
380    CONTINUE

N  = I

IF( I.LT.L .OR. I.GT.K ) EVALD( I ) = AAD( I, I )

390 CONTINUE

N  = K
T  = ZERO

c                                      ** Search for next eigenvalues
400 CONTINUE
IF( N.LT.L ) GO TO  550

IN  = 0
N1  = N - 1
N2  = N - 2
c                          ** Look for single small sub-diagonal element
410 CONTINUE

DO 420 I = L, N
LB  = N + L - I

IF( LB.EQ.L ) GO TO  430

S  = ABS( AAD( LB - 1,LB - 1 ) ) + ABS( AAD( LB,LB ) )

IF( S.EQ.ZERO ) S  = RNORM

IF( ABS( AAD( LB, LB-1 ) ).LE. TOL*S ) GO TO  430

420 CONTINUE

430 CONTINUE
X  = AAD( N, N )

IF( LB.EQ.N ) THEN
c                                        ** One eigenvalue found
AAD( N, N ) = X + T
EVALD( N ) = AAD( N, N )
N  = N1
GO TO  400

END IF

Y  = AAD( N1, N1 )

IF( LB.EQ.N1 ) THEN
c                                        ** Two eigenvalues found
P  = ( Y - X )*C2
Q  = P**2 + W
Z  = SQRT( ABS( Q ) )
AAD( N, N ) = X + T
X  = AAD( N, N )
AAD( N1, N1 ) = Y + T
c                                        ** Real pair
Z  = P + SIGN( Z, P )
EVALD( N1 ) = X + Z
EVALD( N ) = EVALD( N1 )

IF( Z.NE.ZERO ) EVALD( N ) = X - W / Z

X  = AAD( N, N1 )
c                                  ** Employ scale factor in case
c                                  ** X and Z are very small
R  = SQRT( X*X + Z*Z )
P  = X / R
Q  = Z / R
c                                             ** Row modification
DO 440 J = N1, M
Z  = AAD( N1, J )
440    CONTINUE
c                                             ** Column modification
DO 450 I = 1, N
Z  = AAD( I, N1 )
450    CONTINUE
c                                          ** Accumulate transformations
DO 460 I = L, K
Z  = EVECD( I, N1 )
EVECD( I, N1 ) = Q*Z + P*EVECD( I, N )
EVECD( I, N ) = Q*EVECD( I, N ) - P*Z
460    CONTINUE

N  = N2
GO TO  400

END IF

IF( IN.EQ.30 ) THEN

c                    ** No convergence after 30 iterations; set error
c                    ** indicator to the index of the current eigenvalue
IER  = N
GO TO  700

END IF
c                                                  ** Form shift
IF( IN.EQ.10 .OR. IN.EQ.20 ) THEN

T  = T + X

DO 470 I = L, N
470    CONTINUE

S  = ABS( AAD( N,N1 ) ) + ABS( AAD( N1,N2 ) )
X  = C3*S
Y  = X
W  = -C1*S**2

END IF

IN  = IN + 1

c                ** Look for two consecutive small sub-diagonal elements

DO 480 J = LB, N2
I  = N2 + LB - J
Z  = AAD( I, I )
R  = X - Z
S  = Y - Z
P  = ( R*S - W ) / AAD( I + 1, I ) + AAD( I, I + 1 )
Q  = AAD( I + 1, I + 1 ) - Z - R - S
R  = AAD( I + 2, I + 1 )
S  = ABS( P ) + ABS( Q ) + ABS( R )
P  = P / S
Q  = Q / S
R  = R / S

IF( I.EQ.LB ) GO TO  490

UU   = ABS( AAD( I, I-1 ) )*( ABS( Q ) + ABS( R ) )
VV   = ABS( P ) * ( ABS( AAD( I-1, I-1 ) ) + ABS( Z ) +
&                       ABS( AAD( I+1, I+1 ) ) )

IF( UU .LE. TOL*VV ) GO TO  490

480 CONTINUE

490 CONTINUE
AAD( I+2, I ) = ZERO

DO 500 J = I + 3, N
AAD( J, J - 2 ) = ZERO
AAD( J, J - 3 ) = ZERO
500 CONTINUE

c             ** Double QR step involving rows K to N and columns M to N

DO 540 KA = I, N1

NOTLAS = KA.NE.N1

IF( KA.EQ.I ) THEN

S  = SIGN( SQRT( P*P + Q*Q + R*R ), P )

IF( LB.NE.I ) AAD( KA, KA - 1 ) = -AAD( KA, KA - 1 )

ELSE

P  = AAD( KA, KA - 1 )
Q  = AAD( KA + 1, KA - 1 )
R  = ZERO

IF( NOTLAS ) R  = AAD( KA + 2, KA - 1 )

X  = ABS( P ) + ABS( Q ) + ABS( R )

IF( X.EQ.ZERO ) GO TO  540

P  = P / X
Q  = Q / X
R  = R / X
S  = SIGN( SQRT( P*P + Q*Q + R*R ), P )
AAD( KA, KA - 1 ) = -S*X

END IF

P  = P + S
X  = P / S
Y  = Q / S
Z  = R / S
Q  = Q / P
R  = R / P
c                                              ** Row modification
DO 510 J = KA, M

P  = AAD( KA, J ) + Q*AAD( KA + 1, J )

IF( NOTLAS ) THEN

P  = P + R*AAD( KA + 2, J )
AAD( KA + 2, J ) = AAD( KA + 2, J ) - P*Z

END IF

AAD( KA + 1, J ) = AAD( KA + 1, J ) - P*Y
510    CONTINUE
c                                                 ** Column modification
DO 520 II = 1, MIN( N, KA + 3 )

P  = X*AAD( II, KA ) + Y*AAD( II, KA + 1 )

IF( NOTLAS ) THEN

P  = P + Z*AAD( II, KA + 2 )
AAD( II, KA + 2 ) = AAD( II, KA + 2 ) - P*R

END IF

AAD( II, KA + 1 ) = AAD( II, KA + 1 ) - P*Q
520    CONTINUE
c                                          ** Accumulate transformations
DO 530 II = L, K

P  = X*EVECD( II, KA ) + Y*EVECD( II, KA + 1 )

IF( NOTLAS ) THEN

P  = P + Z*EVECD( II, KA + 2 )
EVECD( II, KA + 2 ) = EVECD( II, KA + 2 ) - P*R

END IF

EVECD( II, KA + 1 ) = EVECD( II, KA + 1 ) - P*Q
EVECD( II, KA ) = EVECD( II, KA ) - P
530    CONTINUE

540 CONTINUE

GO TO  410
c                     ** All evals found, now backsubstitute real vector
550 CONTINUE

IF( RNORM.NE.ZERO ) THEN

DO 580 N = M, 1, -1
N2   = N
AAD( N, N ) = ONE

DO 570 I = N - 1, 1, -1
W  = AAD( I, I ) - EVALD( N )

IF( W.EQ.ZERO ) W  = TOL*RNORM

R  = AAD( I, N )

DO 560 J = N2, N - 1
560          CONTINUE

AAD( I, N ) = -R / W
N2   = I
570       CONTINUE

580    CONTINUE
c                      ** End backsubstitution vectors of isolated evals
DO 600 I = 1, M

IF( I.LT.L .OR. I.GT.K ) THEN

DO 590 J = I, M
EVECD( I, J ) = AAD( I, J )
590          CONTINUE

END IF

600    CONTINUE
c                                   ** Multiply by transformation matrix
IF( K.NE.0 ) THEN

DO 630 J = M, L, -1

DO 620 I = L, K
Z  = ZERO

DO 610 N = L, MIN( J, K )
Z  = Z + EVECD( I, N )*AAD( N, J )
610             CONTINUE

EVECD( I, J ) = Z
620          CONTINUE

630       CONTINUE

END IF

END IF

DO 650 I = L, K

DO 640 J = 1, M
EVECD( I, J ) = EVECD( I, J ) * WKD( I )
640    CONTINUE
650 CONTINUE

c                           ** Interchange rows if permutations occurred
DO 670 I = L-1, 1, -1

J  = WKD( I )

IF( I.NE.J ) THEN

DO 660 N = 1, M
REPL   = EVECD( I, N )
EVECD( I, N ) = EVECD( J, N )
EVECD( J, N ) = REPL
660       CONTINUE

END IF

670 CONTINUE

DO 690 I = K + 1, M

J  = WKD( I )

IF( I.NE.J ) THEN

DO 680 N = 1, M
REPL   = EVECD( I, N )
EVECD( I, N ) = EVECD( J, N )
EVECD( J, N ) = REPL
680       CONTINUE

END IF

690 CONTINUE

c                         ** Put results into output arrays
700 CONTINUE

DO 720 J = 1, M

EVAL( J ) = EVALD( J )

DO 710 K = 1, M
EVEC( J, K ) = EVECD( J, K )
710    CONTINUE

720 CONTINUE

RETURN
END