Commit a413d7e6 authored by rtoy's avatar rtoy
Browse files

Initial revision for LAPACK routines.

parent a0ce62d4
;; -*- Mode: lisp -*-
(mk:define-language :f2cl
:compiler #'f2cl:f2cl-compile
:source-extension "f")
(mk:defsystem lapack
:source-pathname (directory-namestring *load-truename*)
:depends-on ("blas")
:compiler-options (:keep-lisp-file t
:relaxed-array-decls t
:array-type :array
:array-slicing t
)
:components
(
(:module lapack
:binary-pathname ".."
:source-extension "f"
:language :f2cl
:components
(
(:file "dgeev"
:depends-on ("dgebak"
"dgebal"
"dgehrd"
"dhseqr"
"dlabad"
"dlacpy"
"dlamch"
"dlange"
"dlapy2"
"dlartg"
"dlascl"
"dorghr"
"dtrevc"
))
(:file "dgebak"
:depends-on ("ilaenv"
"dlamch"
"dlabad"
"dlange"
"dlascl"
"dgebal"
"dgehrd"
"dlacpy"
"dorghr"
"dhseqr"
"dtrevc"
"dlapy2"
"dlartg"
))
(:file "dgebal")
(:file "dgehrd"
:depends-on ("dlahrd"
"dlarfb"
"dgehd2"))
(:file "dhseqr"
:depends-on ("dlaset"
"dlahqr"
"dlanhs"
"dlapy2"
"dlarfg"
"dlarfx"))
(:file "dlabad")
(:file "dlacpy")
(:file "dlamch")
(:file "dlange"
:depends-on ("dlassq"))
(:file "dlapy2")
(:file "dlartg")
(:file "dlascl")
(:file "dorghr"
:depends-on ("dorgqr"))
(:file "dtrevc"
:depends-on ("dlaln2"))
(:file "ilaenv"
:depends-on ("ieeeck"))
(:file "dlahrd"
:depends-on ("dlarfg"))
(:file "dlarfb")
(:file "dgehd2"
:depends-on ("dlarfg"
"dlarf"))
(:file "dlaset")
(:file "dlahqr"
:depends-on ("dlanhs"
"dlarfg"
"dlanv2"))
(:file "dlanhs")
(:file "dlarfg")
(:file "dlarfx")
(:file "dlassq")
(:file "dorgqr"
:depends-on ("dorg2r"
"dlarft"))
(:file "dlaln2"
:depends-on ("dladiv"))
(:file "ieeeck")
(:file "dlarf")
(:file "dlanv2")
(:file "dorg2r")
(:file "dlarft")
(:file "dladiv")
;;(:file "zgeev")
))))
SUBROUTINE DGEBAK( JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV,
$ INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* September 30, 1994
*
* .. Scalar Arguments ..
CHARACTER JOB, SIDE
INTEGER IHI, ILO, INFO, LDV, M, N
* ..
* .. Array Arguments ..
DOUBLE PRECISION SCALE( * ), V( LDV, * )
* ..
*
* Purpose
* =======
*
* DGEBAK forms the right or left eigenvectors of a real general matrix
* by backward transformation on the computed eigenvectors of the
* balanced matrix output by DGEBAL.
*
* Arguments
* =========
*
* JOB (input) CHARACTER*1
* Specifies the type of backward transformation required:
* = 'N', do nothing, return immediately;
* = 'P', do backward transformation for permutation only;
* = 'S', do backward transformation for scaling only;
* = 'B', do backward transformations for both permutation and
* scaling.
* JOB must be the same as the argument JOB supplied to DGEBAL.
*
* SIDE (input) CHARACTER*1
* = 'R': V contains right eigenvectors;
* = 'L': V contains left eigenvectors.
*
* N (input) INTEGER
* The number of rows of the matrix V. N >= 0.
*
* ILO (input) INTEGER
* IHI (input) INTEGER
* The integers ILO and IHI determined by DGEBAL.
* 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
*
* SCALE (input) DOUBLE PRECISION array, dimension (N)
* Details of the permutation and scaling factors, as returned
* by DGEBAL.
*
* M (input) INTEGER
* The number of columns of the matrix V. M >= 0.
*
* V (input/output) DOUBLE PRECISION array, dimension (LDV,M)
* On entry, the matrix of right or left eigenvectors to be
* transformed, as returned by DHSEIN or DTREVC.
* On exit, V is overwritten by the transformed eigenvectors.
*
* LDV (input) INTEGER
* The leading dimension of the array V. LDV >= max(1,N).
*
* INFO (output) INTEGER
* = 0: successful exit
* < 0: if INFO = -i, the i-th argument had an illegal value.
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ONE
PARAMETER ( ONE = 1.0D+0 )
* ..
* .. Local Scalars ..
LOGICAL LEFTV, RIGHTV
INTEGER I, II, K
DOUBLE PRECISION S
* ..
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* ..
* .. External Subroutines ..
EXTERNAL DSCAL, DSWAP, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC MAX, MIN
* ..
* .. Executable Statements ..
*
* Decode and Test the input parameters
*
RIGHTV = LSAME( SIDE, 'R' )
LEFTV = LSAME( SIDE, 'L' )
*
INFO = 0
IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
$ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
INFO = -1
ELSE IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
INFO = -2
ELSE IF( N.LT.0 ) THEN
INFO = -3
ELSE IF( ILO.LT.1 .OR. ILO.GT.MAX( 1, N ) ) THEN
INFO = -4
ELSE IF( IHI.LT.MIN( ILO, N ) .OR. IHI.GT.N ) THEN
INFO = -5
ELSE IF( M.LT.0 ) THEN
INFO = -7
ELSE IF( LDV.LT.MAX( 1, N ) ) THEN
INFO = -9
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGEBAK', -INFO )
RETURN
END IF
*
* Quick return if possible
*
IF( N.EQ.0 )
$ RETURN
IF( M.EQ.0 )
$ RETURN
IF( LSAME( JOB, 'N' ) )
$ RETURN
*
IF( ILO.EQ.IHI )
$ GO TO 30
*
* Backward balance
*
IF( LSAME( JOB, 'S' ) .OR. LSAME( JOB, 'B' ) ) THEN
*
IF( RIGHTV ) THEN
DO 10 I = ILO, IHI
S = SCALE( I )
CALL DSCAL( M, S, V( I, 1 ), LDV )
10 CONTINUE
END IF
*
IF( LEFTV ) THEN
DO 20 I = ILO, IHI
S = ONE / SCALE( I )
CALL DSCAL( M, S, V( I, 1 ), LDV )
20 CONTINUE
END IF
*
END IF
*
* Backward permutation
*
* For I = ILO-1 step -1 until 1,
* IHI+1 step 1 until N do --
*
30 CONTINUE
IF( LSAME( JOB, 'P' ) .OR. LSAME( JOB, 'B' ) ) THEN
IF( RIGHTV ) THEN
DO 40 II = 1, N
I = II
IF( I.GE.ILO .AND. I.LE.IHI )
$ GO TO 40
IF( I.LT.ILO )
$ I = ILO - II
K = SCALE( I )
IF( K.EQ.I )
$ GO TO 40
CALL DSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV )
40 CONTINUE
END IF
*
IF( LEFTV ) THEN
DO 50 II = 1, N
I = II
IF( I.GE.ILO .AND. I.LE.IHI )
$ GO TO 50
IF( I.LT.ILO )
$ I = ILO - II
K = SCALE( I )
IF( K.EQ.I )
$ GO TO 50
CALL DSWAP( M, V( I, 1 ), LDV, V( K, 1 ), LDV )
50 CONTINUE
END IF
END IF
*
RETURN
*
* End of DGEBAK
*
END
SUBROUTINE DGEBAL( JOB, N, A, LDA, ILO, IHI, SCALE, INFO )
*
* -- LAPACK routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* June 30, 1999
*
* .. Scalar Arguments ..
CHARACTER JOB
INTEGER IHI, ILO, INFO, LDA, N
* ..
* .. Array Arguments ..
DOUBLE PRECISION A( LDA, * ), SCALE( * )
* ..
*
* Purpose
* =======
*
* DGEBAL balances a general real matrix A. This involves, first,
* permuting A by a similarity transformation to isolate eigenvalues
* in the first 1 to ILO-1 and last IHI+1 to N elements on the
* diagonal; and second, applying a diagonal similarity transformation
* to rows and columns ILO to IHI to make the rows and columns as
* close in norm as possible. Both steps are optional.
*
* Balancing may reduce the 1-norm of the matrix, and improve the
* accuracy of the computed eigenvalues and/or eigenvectors.
*
* Arguments
* =========
*
* JOB (input) CHARACTER*1
* Specifies the operations to be performed on A:
* = 'N': none: simply set ILO = 1, IHI = N, SCALE(I) = 1.0
* for i = 1,...,N;
* = 'P': permute only;
* = 'S': scale only;
* = 'B': both permute and scale.
*
* N (input) INTEGER
* The order of the matrix A. N >= 0.
*
* A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
* On entry, the input matrix A.
* On exit, A is overwritten by the balanced matrix.
* If JOB = 'N', A is not referenced.
* See Further Details.
*
* LDA (input) INTEGER
* The leading dimension of the array A. LDA >= max(1,N).
*
* ILO (output) INTEGER
* IHI (output) INTEGER
* ILO and IHI are set to integers such that on exit
* A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.
* If JOB = 'N' or 'S', ILO = 1 and IHI = N.
*
* SCALE (output) DOUBLE PRECISION array, dimension (N)
* Details of the permutations and scaling factors applied to
* A. If P(j) is the index of the row and column interchanged
* with row and column j and D(j) is the scaling factor
* applied to row and column j, then
* SCALE(j) = P(j) for j = 1,...,ILO-1
* = D(j) for j = ILO,...,IHI
* = P(j) for j = IHI+1,...,N.
* The order in which the interchanges are made is N to IHI+1,
* then 1 to ILO-1.
*
* INFO (output) INTEGER
* = 0: successful exit.
* < 0: if INFO = -i, the i-th argument had an illegal value.
*
* Further Details
* ===============
*
* The permutations consist of row and column interchanges which put
* the matrix in the form
*
* ( T1 X Y )
* P A P = ( 0 B Z )
* ( 0 0 T2 )
*
* where T1 and T2 are upper triangular matrices whose eigenvalues lie
* along the diagonal. The column indices ILO and IHI mark the starting
* and ending columns of the submatrix B. Balancing consists of applying
* a diagonal similarity transformation inv(D) * B * D to make the
* 1-norms of each row of B and its corresponding column nearly equal.
* The output matrix is
*
* ( T1 X*D Y )
* ( 0 inv(D)*B*D inv(D)*Z ).
* ( 0 0 T2 )
*
* Information about the permutations P and the diagonal matrix D is
* returned in the vector SCALE.
*
* This subroutine is based on the EISPACK routine BALANC.
*
* Modified by Tzu-Yi Chen, Computer Science Division, University of
* California at Berkeley, USA
*
* =====================================================================
*
* .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
DOUBLE PRECISION SCLFAC
PARAMETER ( SCLFAC = 0.8D+1 )
DOUBLE PRECISION FACTOR
PARAMETER ( FACTOR = 0.95D+0 )
* ..
* .. Local Scalars ..
LOGICAL NOCONV
INTEGER I, ICA, IEXC, IRA, J, K, L, M
DOUBLE PRECISION C, CA, F, G, R, RA, S, SFMAX1, SFMAX2, SFMIN1,
$ SFMIN2
* ..
* .. External Functions ..
LOGICAL LSAME
INTEGER IDAMAX
DOUBLE PRECISION DLAMCH
EXTERNAL LSAME, IDAMAX, DLAMCH
* ..
* .. External Subroutines ..
EXTERNAL DSCAL, DSWAP, XERBLA
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN
* ..
* .. Executable Statements ..
*
* Test the input parameters
*
INFO = 0
IF( .NOT.LSAME( JOB, 'N' ) .AND. .NOT.LSAME( JOB, 'P' ) .AND.
$ .NOT.LSAME( JOB, 'S' ) .AND. .NOT.LSAME( JOB, 'B' ) ) THEN
INFO = -1
ELSE IF( N.LT.0 ) THEN
INFO = -2
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -4
END IF
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'DGEBAL', -INFO )
RETURN
END IF
*
K = 1
L = N
*
IF( N.EQ.0 )
$ GO TO 210
*
IF( LSAME( JOB, 'N' ) ) THEN
DO 10 I = 1, N
SCALE( I ) = ONE
10 CONTINUE
GO TO 210
END IF
*
IF( LSAME( JOB, 'S' ) )
$ GO TO 120
*
* Permutation to isolate eigenvalues if possible
*
GO TO 50
*
* Row and column exchange.
*
20 CONTINUE
SCALE( M ) = J
IF( J.EQ.M )
$ GO TO 30
*
CALL DSWAP( L, A( 1, J ), 1, A( 1, M ), 1 )
CALL DSWAP( N-K+1, A( J, K ), LDA, A( M, K ), LDA )
*
30 CONTINUE
GO TO ( 40, 80 )IEXC
*
* Search for rows isolating an eigenvalue and push them down.
*
40 CONTINUE
IF( L.EQ.1 )
$ GO TO 210
L = L - 1
*
50 CONTINUE
DO 70 J = L, 1, -1
*
DO 60 I = 1, L
IF( I.EQ.J )
$ GO TO 60
IF( A( J, I ).NE.ZERO )
$ GO TO 70
60 CONTINUE
*
M = L
IEXC = 1
GO TO 20
70 CONTINUE
*
GO TO 90
*
* Search for columns isolating an eigenvalue and push them left.
*
80 CONTINUE
K = K + 1
*
90 CONTINUE
DO 110 J = K, L
*
DO 100 I = K, L
IF( I.EQ.J )
$ GO TO 100
IF( A( I, J ).NE.ZERO )
$ GO TO 110
100 CONTINUE
*
M = K
IEXC = 2
GO TO 20
110 CONTINUE
*
120 CONTINUE
DO 130 I = K, L
SCALE( I ) = ONE
130 CONTINUE
*
IF( LSAME( JOB, 'P' ) )
$ GO TO 210
*
* Balance the submatrix in rows K to L.
*
* Iterative loop for norm reduction
*
SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
SFMAX1 = ONE / SFMIN1
SFMIN2 = SFMIN1*SCLFAC
SFMAX2 = ONE / SFMIN2
140 CONTINUE
NOCONV = .FALSE.
*
DO 200 I = K, L
C = ZERO
R = ZERO
*
DO 150 J = K, L
IF( J.EQ.I )
$ GO TO 150
C = C + ABS( A( J, I ) )
R = R + ABS( A( I, J ) )
150 CONTINUE
ICA = IDAMAX( L, A( 1, I ), 1 )
CA = ABS( A( ICA, I ) )
IRA = IDAMAX( N-K+1, A( I, K ), LDA )
RA = ABS( A( I, IRA+K-1 ) )
*
* Guard against zero C or R due to underflow.
*
IF( C.EQ.ZERO .OR. R.EQ.ZERO )
$ GO TO 200
G = R / SCLFAC
F = ONE
S = C + R
160 CONTINUE
IF( C.GE.G .OR. MAX( F, C, CA ).GE.SFMAX2 .OR.
$ MIN( R, G, RA ).LE.SFMIN2 )GO TO 170
F = F*SCLFAC
C = C*SCLFAC
CA = CA*SCLFAC
R = R / SCLFAC
G = G / SCLFAC
RA = RA / SCLFAC
GO TO 160
*
170 CONTINUE