Commit dfed4896 authored by rtoy's avatar rtoy
Browse files

packages/lapack.system:

o Add DGEEVX and dependencies

packages/lapack/lapack-tests.lisp:
o Add test for DGEEVX
o Add comments

packages/lapack/dgeevx.f:
packages/lapack/dlacon.f:
packages/lapack/dlaexc.f:
packages/lapack/dlaqtr.f:
packages/lapack/dlasy2.f:
packages/lapack/dtrexc.f:
packages/lapack/dtrsna.f:
o New files for DGEEVX and dependencies.
parent ae180aa2
......@@ -116,6 +116,17 @@
(:file "dorg2r")
(:file "dlarft")
(:file "dladiv")
(:file "dgeevx"
:depends-on ("dtrsna"))
(:file "dtrsna"
:depends-on ("dtrexc" "dlacon" "dlaqtr"))
(:file "dtrexc"
:depends-on ("dlaexc"))
(:file "dlacon")
(:file "dlaqtr")
(:file "dlaexc"
:depends-on ("dlasy2"))
(:file "dlasy2")
;;(:file "zgeev")
))))
......
This diff is collapsed.
SUBROUTINE DLACON( N, V, X, ISGN, EST, KASE )
*
* -- LAPACK auxiliary routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* February 29, 1992
*
* .. Scalar Arguments ..
INTEGER KASE, N
DOUBLE PRECISION EST
* ..
* .. Array Arguments ..
INTEGER ISGN( * )
DOUBLE PRECISION V( * ), X( * )
* ..
*
* Purpose
* =======
*
* DLACON estimates the 1-norm of a square, real matrix A.
* Reverse communication is used for evaluating matrix-vector products.
*
* Arguments
* =========
*
* N (input) INTEGER
* The order of the matrix. N >= 1.
*
* V (workspace) DOUBLE PRECISION array, dimension (N)
* On the final return, V = A*W, where EST = norm(V)/norm(W)
* (W is not returned).
*
* X (input/output) DOUBLE PRECISION array, dimension (N)
* On an intermediate return, X should be overwritten by
* A * X, if KASE=1,
* A' * X, if KASE=2,
* and DLACON must be re-called with all the other parameters
* unchanged.
*
* ISGN (workspace) INTEGER array, dimension (N)
*
* EST (output) DOUBLE PRECISION
* An estimate (a lower bound) for norm(A).
*
* KASE (input/output) INTEGER
* On the initial call to DLACON, KASE should be 0.
* On an intermediate return, KASE will be 1 or 2, indicating
* whether X should be overwritten by A * X or A' * X.
* On the final return from DLACON, KASE will again be 0.
*
* Further Details
* ======= =======
*
* Contributed by Nick Higham, University of Manchester.
* Originally named SONEST, dated March 16, 1988.
*
* Reference: N.J. Higham, "FORTRAN codes for estimating the one-norm of
* a real or complex matrix, with applications to condition estimation",
* ACM Trans. Math. Soft., vol. 14, no. 4, pp. 381-396, December 1988.
*
* =====================================================================
*
* .. Parameters ..
INTEGER ITMAX
PARAMETER ( ITMAX = 5 )
DOUBLE PRECISION ZERO, ONE, TWO
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0 )
* ..
* .. Local Scalars ..
INTEGER I, ITER, J, JLAST, JUMP
DOUBLE PRECISION ALTSGN, ESTOLD, TEMP
* ..
* .. External Functions ..
INTEGER IDAMAX
DOUBLE PRECISION DASUM
EXTERNAL IDAMAX, DASUM
* ..
* .. External Subroutines ..
EXTERNAL DCOPY
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, DBLE, NINT, SIGN
* ..
* .. Save statement ..
SAVE
* ..
* .. Executable Statements ..
*
IF( KASE.EQ.0 ) THEN
DO 10 I = 1, N
X( I ) = ONE / DBLE( N )
10 CONTINUE
KASE = 1
JUMP = 1
RETURN
END IF
*
GO TO ( 20, 40, 70, 110, 140 )JUMP
*
* ................ ENTRY (JUMP = 1)
* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY A*X.
*
20 CONTINUE
IF( N.EQ.1 ) THEN
V( 1 ) = X( 1 )
EST = ABS( V( 1 ) )
* ... QUIT
GO TO 150
END IF
EST = DASUM( N, X, 1 )
*
DO 30 I = 1, N
X( I ) = SIGN( ONE, X( I ) )
ISGN( I ) = NINT( X( I ) )
30 CONTINUE
KASE = 2
JUMP = 2
RETURN
*
* ................ ENTRY (JUMP = 2)
* FIRST ITERATION. X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
*
40 CONTINUE
J = IDAMAX( N, X, 1 )
ITER = 2
*
* MAIN LOOP - ITERATIONS 2,3,...,ITMAX.
*
50 CONTINUE
DO 60 I = 1, N
X( I ) = ZERO
60 CONTINUE
X( J ) = ONE
KASE = 1
JUMP = 3
RETURN
*
* ................ ENTRY (JUMP = 3)
* X HAS BEEN OVERWRITTEN BY A*X.
*
70 CONTINUE
CALL DCOPY( N, X, 1, V, 1 )
ESTOLD = EST
EST = DASUM( N, V, 1 )
DO 80 I = 1, N
IF( NINT( SIGN( ONE, X( I ) ) ).NE.ISGN( I ) )
$ GO TO 90
80 CONTINUE
* REPEATED SIGN VECTOR DETECTED, HENCE ALGORITHM HAS CONVERGED.
GO TO 120
*
90 CONTINUE
* TEST FOR CYCLING.
IF( EST.LE.ESTOLD )
$ GO TO 120
*
DO 100 I = 1, N
X( I ) = SIGN( ONE, X( I ) )
ISGN( I ) = NINT( X( I ) )
100 CONTINUE
KASE = 2
JUMP = 4
RETURN
*
* ................ ENTRY (JUMP = 4)
* X HAS BEEN OVERWRITTEN BY TRANDPOSE(A)*X.
*
110 CONTINUE
JLAST = J
J = IDAMAX( N, X, 1 )
IF( ( X( JLAST ).NE.ABS( X( J ) ) ) .AND. ( ITER.LT.ITMAX ) ) THEN
ITER = ITER + 1
GO TO 50
END IF
*
* ITERATION COMPLETE. FINAL STAGE.
*
120 CONTINUE
ALTSGN = ONE
DO 130 I = 1, N
X( I ) = ALTSGN*( ONE+DBLE( I-1 ) / DBLE( N-1 ) )
ALTSGN = -ALTSGN
130 CONTINUE
KASE = 1
JUMP = 5
RETURN
*
* ................ ENTRY (JUMP = 5)
* X HAS BEEN OVERWRITTEN BY A*X.
*
140 CONTINUE
TEMP = TWO*( DASUM( N, X, 1 ) / DBLE( 3*N ) )
IF( TEMP.GT.EST ) THEN
CALL DCOPY( N, X, 1, V, 1 )
EST = TEMP
END IF
*
150 CONTINUE
KASE = 0
RETURN
*
* End of DLACON
*
END
SUBROUTINE DLAEXC( WANTQ, N, T, LDT, Q, LDQ, J1, N1, N2, WORK,
$ INFO )
*
* -- LAPACK auxiliary routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* February 29, 1992
*
* .. Scalar Arguments ..
LOGICAL WANTQ
INTEGER INFO, J1, LDQ, LDT, N, N1, N2
* ..
* .. Array Arguments ..
DOUBLE PRECISION Q( LDQ, * ), T( LDT, * ), WORK( * )
* ..
*
* Purpose
* =======
*
* DLAEXC swaps adjacent diagonal blocks T11 and T22 of order 1 or 2 in
* an upper quasi-triangular matrix T by an orthogonal similarity
* transformation.
*
* T must be in Schur canonical form, that is, block upper triangular
* with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2 diagonal block
* has its diagonal elemnts equal and its off-diagonal elements of
* opposite sign.
*
* Arguments
* =========
*
* WANTQ (input) LOGICAL
* = .TRUE. : accumulate the transformation in the matrix Q;
* = .FALSE.: do not accumulate the transformation.
*
* N (input) INTEGER
* The order of the matrix T. N >= 0.
*
* T (input/output) DOUBLE PRECISION array, dimension (LDT,N)
* On entry, the upper quasi-triangular matrix T, in Schur
* canonical form.
* On exit, the updated matrix T, again in Schur canonical form.
*
* LDT (input) INTEGER
* The leading dimension of the array T. LDT >= max(1,N).
*
* Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N)
* On entry, if WANTQ is .TRUE., the orthogonal matrix Q.
* On exit, if WANTQ is .TRUE., the updated matrix Q.
* If WANTQ is .FALSE., Q is not referenced.
*
* LDQ (input) INTEGER
* The leading dimension of the array Q.
* LDQ >= 1; and if WANTQ is .TRUE., LDQ >= N.
*
* J1 (input) INTEGER
* The index of the first row of the first block T11.
*
* N1 (input) INTEGER
* The order of the first block T11. N1 = 0, 1 or 2.
*
* N2 (input) INTEGER
* The order of the second block T22. N2 = 0, 1 or 2.
*
* WORK (workspace) DOUBLE PRECISION array, dimension (N)
*
* INFO (output) INTEGER
* = 0: successful exit
* = 1: the transformed matrix T would be too far from Schur
* form; the blocks are not swapped and T and Q are
* unchanged.
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
DOUBLE PRECISION TEN
PARAMETER ( TEN = 1.0D+1 )
INTEGER LDD, LDX
PARAMETER ( LDD = 4, LDX = 2 )
* ..
* .. Local Scalars ..
INTEGER IERR, J2, J3, J4, K, ND
DOUBLE PRECISION CS, DNORM, EPS, SCALE, SMLNUM, SN, T11, T22,
$ T33, TAU, TAU1, TAU2, TEMP, THRESH, WI1, WI2,
$ WR1, WR2, XNORM
* ..
* .. Local Arrays ..
DOUBLE PRECISION D( LDD, 4 ), U( 3 ), U1( 3 ), U2( 3 ),
$ X( LDX, 2 )
* ..
* .. External Functions ..
DOUBLE PRECISION DLAMCH, DLANGE
EXTERNAL DLAMCH, DLANGE
* ..
* .. External Subroutines ..
EXTERNAL DLACPY, DLANV2, DLARFG, DLARFX, DLARTG, DLASY2,
$ DROT
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX
* ..
* .. Executable Statements ..
*
INFO = 0
*
* Quick return if possible
*
IF( N.EQ.0 .OR. N1.EQ.0 .OR. N2.EQ.0 )
$ RETURN
IF( J1+N1.GT.N )
$ RETURN
*
J2 = J1 + 1
J3 = J1 + 2
J4 = J1 + 3
*
IF( N1.EQ.1 .AND. N2.EQ.1 ) THEN
*
* Swap two 1-by-1 blocks.
*
T11 = T( J1, J1 )
T22 = T( J2, J2 )
*
* Determine the transformation to perform the interchange.
*
CALL DLARTG( T( J1, J2 ), T22-T11, CS, SN, TEMP )
*
* Apply transformation to the matrix T.
*
IF( J3.LE.N )
$ CALL DROT( N-J1-1, T( J1, J3 ), LDT, T( J2, J3 ), LDT, CS,
$ SN )
CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
*
T( J1, J1 ) = T22
T( J2, J2 ) = T11
*
IF( WANTQ ) THEN
*
* Accumulate transformation in the matrix Q.
*
CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
END IF
*
ELSE
*
* Swapping involves at least one 2-by-2 block.
*
* Copy the diagonal block of order N1+N2 to the local array D
* and compute its norm.
*
ND = N1 + N2
CALL DLACPY( 'Full', ND, ND, T( J1, J1 ), LDT, D, LDD )
DNORM = DLANGE( 'Max', ND, ND, D, LDD, WORK )
*
* Compute machine-dependent threshold for test for accepting
* swap.
*
EPS = DLAMCH( 'P' )
SMLNUM = DLAMCH( 'S' ) / EPS
THRESH = MAX( TEN*EPS*DNORM, SMLNUM )
*
* Solve T11*X - X*T22 = scale*T12 for X.
*
CALL DLASY2( .FALSE., .FALSE., -1, N1, N2, D, LDD,
$ D( N1+1, N1+1 ), LDD, D( 1, N1+1 ), LDD, SCALE, X,
$ LDX, XNORM, IERR )
*
* Swap the adjacent diagonal blocks.
*
K = N1 + N1 + N2 - 3
GO TO ( 10, 20, 30 )K
*
10 CONTINUE
*
* N1 = 1, N2 = 2: generate elementary reflector H so that:
*
* ( scale, X11, X12 ) H = ( 0, 0, * )
*
U( 1 ) = SCALE
U( 2 ) = X( 1, 1 )
U( 3 ) = X( 1, 2 )
CALL DLARFG( 3, U( 3 ), U, 1, TAU )
U( 3 ) = ONE
T11 = T( J1, J1 )
*
* Perform swap provisionally on diagonal block in D.
*
CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
*
* Test whether to reject swap.
*
IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 3,
$ 3 )-T11 ) ).GT.THRESH )GO TO 50
*
* Accept swap: apply transformation to the entire matrix T.
*
CALL DLARFX( 'L', 3, N-J1+1, U, TAU, T( J1, J1 ), LDT, WORK )
CALL DLARFX( 'R', J2, 3, U, TAU, T( 1, J1 ), LDT, WORK )
*
T( J3, J1 ) = ZERO
T( J3, J2 ) = ZERO
T( J3, J3 ) = T11
*
IF( WANTQ ) THEN
*
* Accumulate transformation in the matrix Q.
*
CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
END IF
GO TO 40
*
20 CONTINUE
*
* N1 = 2, N2 = 1: generate elementary reflector H so that:
*
* H ( -X11 ) = ( * )
* ( -X21 ) = ( 0 )
* ( scale ) = ( 0 )
*
U( 1 ) = -X( 1, 1 )
U( 2 ) = -X( 2, 1 )
U( 3 ) = SCALE
CALL DLARFG( 3, U( 1 ), U( 2 ), 1, TAU )
U( 1 ) = ONE
T33 = T( J3, J3 )
*
* Perform swap provisionally on diagonal block in D.
*
CALL DLARFX( 'L', 3, 3, U, TAU, D, LDD, WORK )
CALL DLARFX( 'R', 3, 3, U, TAU, D, LDD, WORK )
*
* Test whether to reject swap.
*
IF( MAX( ABS( D( 2, 1 ) ), ABS( D( 3, 1 ) ), ABS( D( 1,
$ 1 )-T33 ) ).GT.THRESH )GO TO 50
*
* Accept swap: apply transformation to the entire matrix T.
*
CALL DLARFX( 'R', J3, 3, U, TAU, T( 1, J1 ), LDT, WORK )
CALL DLARFX( 'L', 3, N-J1, U, TAU, T( J1, J2 ), LDT, WORK )
*
T( J1, J1 ) = T33
T( J2, J1 ) = ZERO
T( J3, J1 ) = ZERO
*
IF( WANTQ ) THEN
*
* Accumulate transformation in the matrix Q.
*
CALL DLARFX( 'R', N, 3, U, TAU, Q( 1, J1 ), LDQ, WORK )
END IF
GO TO 40
*
30 CONTINUE
*
* N1 = 2, N2 = 2: generate elementary reflectors H(1) and H(2) so
* that:
*
* H(2) H(1) ( -X11 -X12 ) = ( * * )
* ( -X21 -X22 ) ( 0 * )
* ( scale 0 ) ( 0 0 )
* ( 0 scale ) ( 0 0 )
*
U1( 1 ) = -X( 1, 1 )
U1( 2 ) = -X( 2, 1 )
U1( 3 ) = SCALE
CALL DLARFG( 3, U1( 1 ), U1( 2 ), 1, TAU1 )
U1( 1 ) = ONE
*
TEMP = -TAU1*( X( 1, 2 )+U1( 2 )*X( 2, 2 ) )
U2( 1 ) = -TEMP*U1( 2 ) - X( 2, 2 )
U2( 2 ) = -TEMP*U1( 3 )
U2( 3 ) = SCALE
CALL DLARFG( 3, U2( 1 ), U2( 2 ), 1, TAU2 )
U2( 1 ) = ONE
*
* Perform swap provisionally on diagonal block in D.
*
CALL DLARFX( 'L', 3, 4, U1, TAU1, D, LDD, WORK )
CALL DLARFX( 'R', 4, 3, U1, TAU1, D, LDD, WORK )
CALL DLARFX( 'L', 3, 4, U2, TAU2, D( 2, 1 ), LDD, WORK )
CALL DLARFX( 'R', 4, 3, U2, TAU2, D( 1, 2 ), LDD, WORK )
*
* Test whether to reject swap.
*
IF( MAX( ABS( D( 3, 1 ) ), ABS( D( 3, 2 ) ), ABS( D( 4, 1 ) ),
$ ABS( D( 4, 2 ) ) ).GT.THRESH )GO TO 50
*
* Accept swap: apply transformation to the entire matrix T.
*
CALL DLARFX( 'L', 3, N-J1+1, U1, TAU1, T( J1, J1 ), LDT, WORK )
CALL DLARFX( 'R', J4, 3, U1, TAU1, T( 1, J1 ), LDT, WORK )
CALL DLARFX( 'L', 3, N-J1+1, U2, TAU2, T( J2, J1 ), LDT, WORK )
CALL DLARFX( 'R', J4, 3, U2, TAU2, T( 1, J2 ), LDT, WORK )
*
T( J3, J1 ) = ZERO
T( J3, J2 ) = ZERO
T( J4, J1 ) = ZERO
T( J4, J2 ) = ZERO
*
IF( WANTQ ) THEN
*
* Accumulate transformation in the matrix Q.
*
CALL DLARFX( 'R', N, 3, U1, TAU1, Q( 1, J1 ), LDQ, WORK )
CALL DLARFX( 'R', N, 3, U2, TAU2, Q( 1, J2 ), LDQ, WORK )
END IF
*
40 CONTINUE
*
IF( N2.EQ.2 ) THEN
*
* Standardize new 2-by-2 block T11
*
CALL DLANV2( T( J1, J1 ), T( J1, J2 ), T( J2, J1 ),
$ T( J2, J2 ), WR1, WI1, WR2, WI2, CS, SN )
CALL DROT( N-J1-1, T( J1, J1+2 ), LDT, T( J2, J1+2 ), LDT,
$ CS, SN )
CALL DROT( J1-1, T( 1, J1 ), 1, T( 1, J2 ), 1, CS, SN )
IF( WANTQ )
$ CALL DROT( N, Q( 1, J1 ), 1, Q( 1, J2 ), 1, CS, SN )
END IF
*
IF( N1.EQ.2 ) THEN
*
* Standardize new 2-by-2 block T22
*
J3 = J1 + N2
J4 = J3 + 1
CALL DLANV2( T( J3, J3 ), T( J3, J4 ), T( J4, J3 ),