341 lines
11 KiB
Fortran
341 lines
11 KiB
Fortran
SUBROUTINE AHCON (SIZE,N,M,A,B,OLEVR,OLEVI,CLEVR,CLEVI, TRUNCATED
|
|
& SCR1,SCR2,IPVT,JPVT,CON,WORK,ISEED,IERR) !Test inline comment
|
|
C
|
|
C FUNCTION:
|
|
CF
|
|
CF Determines whether the pair (A,B) is controllable and flags
|
|
CF the eigenvalues corresponding to uncontrollable modes.
|
|
CF this ad-hoc controllability calculation uses a random matrix F
|
|
CF and computes whether eigenvalues move from A to the controlled
|
|
CF system A+B*F.
|
|
CF
|
|
C USAGE:
|
|
CU
|
|
CU CALL AHCON (SIZE,N,M,A,B,OLEVR,OLEVI,CLEVR,CLEVI,SCR1,SCR2,IPVT,
|
|
CU JPVT,CON,WORK,ISEED,IERR)
|
|
CU
|
|
CU since AHCON generates different random F matrices for each
|
|
CU call, as long as iseed is not re-initialized by the main
|
|
CU program, and since this code has the potential to be fooled
|
|
CU by extremely ill-conditioned problems, the cautious user
|
|
CU may wish to call it multiple times and rely, perhaps, on
|
|
CU a 2-of-3 vote. We believe, but have not proved, that any
|
|
CU errors this routine may produce are conservative--i.e., that
|
|
CU it may flag a controllable mode as uncontrollable, but
|
|
CU not vice-versa.
|
|
CU
|
|
C INPUTS:
|
|
CI
|
|
CI SIZE integer - first dimension of all 2-d arrays.
|
|
CI
|
|
CI N integer - number of states.
|
|
CI
|
|
CI M integer - number of inputs.
|
|
CI
|
|
CI A double precision - SIZE by N array containing the
|
|
CI N by N system dynamics matrix A.
|
|
CI
|
|
CI B double precision - SIZE by M array containing the
|
|
CI N by M system input matrix B.
|
|
CI
|
|
CI ISEED initial seed for random number generator; if ISEED=0,
|
|
CI then AHCON will set ISEED to a legal value.
|
|
CI
|
|
C OUTPUTS:
|
|
CO
|
|
CO OLEVR double precision - N dimensional vector containing the
|
|
CO real parts of the eigenvalues of A.
|
|
CO
|
|
CO OLEVI double precision - N dimensional vector containing the
|
|
CO imaginary parts of the eigenvalues of A.
|
|
CO
|
|
CO CLEVR double precision - N dimensional vector work space
|
|
CO containing the real parts of the eigenvalues of A+B*F,
|
|
CO where F is the random matrix.
|
|
CO
|
|
CO CLEVI double precision - N dimensional vector work space
|
|
CO containing the imaginary parts of the eigenvalues of
|
|
CO A+B*F, where F is the random matrix.
|
|
CO
|
|
CO SCR1 double precision - N dimensional vector containing the
|
|
CO magnitudes of the corresponding eigenvalues of A.
|
|
CO
|
|
CO SCR2 double precision - N dimensional vector containing the
|
|
CO damping factors of the corresponding eigenvalues of A.
|
|
CO
|
|
CO IPVT integer - N dimensional vector; contains the row pivots
|
|
CO used in finding the nearest neighbor eigenvalues between
|
|
CO those of A and of A+B*F. The IPVT(1)th eigenvalue of
|
|
CO A and the JPVT(1)th eigenvalue of A+B*F are the closest
|
|
CO pair.
|
|
CO
|
|
CO JPVT integer - N dimensional vector; contains the column
|
|
CO pivots used in finding the nearest neighbor eigenvalues;
|
|
CO see IPVT.
|
|
CO
|
|
CO CON logical - N dimensional vector; flagging the uncontrollable
|
|
CO modes of the system. CON(I)=.TRUE. implies the
|
|
CO eigenvalue of A given by DCMPLX(OLEVR(IPVT(I)),OLEVI(IPVT(i)))
|
|
CO corresponds to a controllable mode; CON(I)=.FALSE.
|
|
CO implies an uncontrollable mode for that eigenvalue.
|
|
CO
|
|
CO WORK double precision - SIZE by N dimensional array containing
|
|
CO an N by N matrix. WORK(I,J) is the distance between
|
|
CO the open loop eigenvalue given by DCMPLX(OLEVR(I),OLEVI(I))
|
|
CO and the closed loop eigenvalue of A+B*F given by
|
|
CO DCMPLX(CLEVR(J),CLEVI(J)).
|
|
CO
|
|
CO IERR integer - IERR=0 indicates normal return; a non-zero
|
|
CO value indicates trouble in the eigenvalue calculation.
|
|
CO see the EISPACK and EIGEN documentation for details.
|
|
CO
|
|
C ALGORITHM:
|
|
CA
|
|
CA Calculate eigenvalues of A and of A+B*F for a randomly
|
|
CA generated F, and see which ones change. Use a full pivot
|
|
CA search through a matrix of euclidean distance measures
|
|
CA between each pair of eigenvalues from (A,A+BF) to
|
|
CA determine the closest pairs.
|
|
CA
|
|
C MACHINE DEPENDENCIES:
|
|
CM
|
|
CM NONE
|
|
CM
|
|
C HISTORY:
|
|
CH
|
|
CH written by: Birdwell & Laub
|
|
CH date: May 18, 1985
|
|
CH current version: 1.0
|
|
CH modifications: made machine independent and modified for
|
|
CH f77:bb:8-86.
|
|
CH changed cmplx -> dcmplx: 7/27/88 jdb
|
|
CH
|
|
C ROUTINES CALLED:
|
|
CC
|
|
CC EIGEN,RAND
|
|
CC
|
|
C COMMON MEMORY USED:
|
|
CM
|
|
CM none
|
|
CM
|
|
C----------------------------------------------------------------------
|
|
C written for: The CASCADE Project
|
|
C Oak Ridge National Laboratory
|
|
C U.S. Department of Energy
|
|
C contract number DE-AC05-840R21400
|
|
C subcontract number 37B-7685 S13
|
|
C organization: The University of Tennessee
|
|
C----------------------------------------------------------------------
|
|
C THIS SOFTWARE IS IN THE PUBLIC DOMAIN
|
|
C NO RESTRICTIONS ON ITS USE ARE IMPLIED
|
|
C----------------------------------------------------------------------
|
|
C
|
|
C--global variables:
|
|
C
|
|
INTEGER SIZE
|
|
INTEGER N
|
|
INTEGER M
|
|
INTEGER IPVT(1)
|
|
INTEGER JPVT(1)
|
|
INTEGER IERR
|
|
C
|
|
DOUBLE PRECISION A(SIZE,N)
|
|
DOUBLE PRECISION B(SIZE,M)
|
|
DOUBLE PRECISION WORK(SIZE,N)
|
|
DOUBLE PRECISION CLEVR(N)
|
|
DOUBLE PRECISION CLEVI(N)
|
|
DOUBLE PRECISION OLEVR(N)
|
|
DOUBLE PRECISION OLEVI(N)
|
|
DOUBLE PRECISION SCR1(N)
|
|
DOUBLE PRECISION SCR2(N)
|
|
C
|
|
LOGICAL CON(N)
|
|
C
|
|
C--local variables:
|
|
C
|
|
INTEGER ISEED
|
|
INTEGER ITEMP
|
|
INTEGER K1
|
|
INTEGER K2
|
|
INTEGER I
|
|
INTEGER J
|
|
INTEGER K
|
|
INTEGER IMAX
|
|
INTEGER JMAX
|
|
C
|
|
DOUBLE PRECISION VALUE
|
|
DOUBLE PRECISION EPS
|
|
DOUBLE PRECISION EPS1
|
|
DOUBLE PRECISION TEMP
|
|
DOUBLE PRECISION CURR
|
|
DOUBLE PRECISION ANORM
|
|
DOUBLE PRECISION BNORM
|
|
DOUBLE PRECISION COLNRM
|
|
DOUBLE PRECISION RNDMNO
|
|
C
|
|
DOUBLE COMPLEX DCMPLX
|
|
C
|
|
C--compute machine epsilon
|
|
C
|
|
EPS = 1.D0
|
|
100 CONTINUE
|
|
EPS = EPS / 2.D0
|
|
EPS1 = 1.D0 + EPS
|
|
IF (EPS1 .NE. 1.D0) GO TO 100
|
|
EPS = EPS * 2.D0
|
|
C
|
|
C--compute the l-1 norm of a
|
|
C
|
|
ANORM = 0.0D0
|
|
DO 120 J = 1, N
|
|
COLNRM = 0.D0
|
|
DO 110 I = 1, N
|
|
COLNRM = COLNRM + ABS(A(I,J))
|
|
110 CONTINUE
|
|
IF (COLNRM .GT. ANORM) ANORM = COLNRM
|
|
120 CONTINUE
|
|
C
|
|
C--compute the l-1 norm of b
|
|
C
|
|
BNORM = 0.0D0
|
|
DO 140 J = 1, M
|
|
COLNRM = 0.D0
|
|
DO 130 I = 1, N
|
|
COLNRM = COLNRM + ABS(B(I,J))
|
|
130 CONTINUE
|
|
IF (COLNRM .GT. BNORM) BNORM = COLNRM
|
|
140 CONTINUE
|
|
C
|
|
C--compute a + b * f
|
|
C
|
|
DO 160 J = 1, N
|
|
DO 150 I = 1, N
|
|
WORK(I,J) = A(I,J)
|
|
150 CONTINUE
|
|
160 CONTINUE
|
|
C
|
|
C--the elements of f are random with uniform distribution
|
|
C--from -anorm/bnorm to +anorm/bnorm
|
|
C--note that f is not explicitly stored as a matrix
|
|
C--pathalogical floating point notes: the if (bnorm .gt. 0.d0)
|
|
C--test should actually be if (bnorm .gt. dsmall), where dsmall
|
|
C--is the smallest representable number whose reciprocal does
|
|
C--not generate an overflow or loss of precision.
|
|
C
|
|
IF (ISEED .EQ. 0) ISEED = 86345823
|
|
IF (ANORM .EQ. 0.D0) ANORM = 1.D0
|
|
IF (BNORM .GT. 0.D0) THEN
|
|
TEMP = 2.D0 * ANORM / BNORM
|
|
ELSE
|
|
TEMP = 2.D0
|
|
END IF
|
|
DO 190 K = 1, M
|
|
DO 180 J = 1, N
|
|
CALL RAND(ISEED,ISEED,RNDMNO)
|
|
VALUE = (RNDMNO - 0.5D0) * TEMP
|
|
DO 170 I = 1, N
|
|
WORK(I,J) = WORK(I,J) + B(I,K)*VALUE
|
|
170 CONTINUE
|
|
180 CONTINUE
|
|
190 CONTINUE
|
|
C
|
|
C--compute the eigenvalues of a + b*f, and several other things
|
|
C
|
|
CALL EIGEN (0,SIZE,N,WORK,CLEVR,CLEVI,WORK,SCR1,SCR2,IERR)
|
|
IF (IERR .NE. 0) RETURN
|
|
C
|
|
C--copy a so it is not destroyed
|
|
C
|
|
DO 210 J = 1, N
|
|
DO 200 I = 1, N
|
|
WORK(I,J) = A(I,J)
|
|
200 CONTINUE
|
|
210 CONTINUE
|
|
C
|
|
C--compute the eigenvalues of a, and several other things
|
|
C
|
|
CALL EIGEN (0,SIZE,N,WORK,OLEVR,OLEVI,WORK,SCR1,SCR2,IERR)
|
|
IF (IERR .NE. 0) RETURN
|
|
C
|
|
C--form the matrix of distances between eigenvalues of a and
|
|
C--EIGENVALUES OF A+B*F
|
|
C
|
|
DO 230 J = 1, N
|
|
DO 220 I = 1, N
|
|
WORK(I,J) =
|
|
& ABS(DCMPLX(OLEVR(I),OLEVI(I))-DCMPLX(CLEVR(J),CLEVI(J)))
|
|
220 CONTINUE
|
|
230 CONTINUE
|
|
C
|
|
C--initialize row and column pivots
|
|
C
|
|
DO 240 I = 1, N
|
|
IPVT(I) = I
|
|
JPVT(I) = I
|
|
240 CONTINUE
|
|
C
|
|
C--a little bit messy to avoid swapping columns and
|
|
C--rows of work
|
|
C
|
|
DO 270 I = 1, N-1
|
|
C
|
|
C--find the minimum element of each lower right square
|
|
C--submatrix of work, for submatrices of size n x n
|
|
C--through 2 x 2
|
|
C
|
|
CURR = WORK(IPVT(I),JPVT(I))
|
|
IMAX = I
|
|
JMAX = I
|
|
TEMP = CURR
|
|
C
|
|
C--find the minimum element
|
|
C
|
|
DO 260 K1 = I, N
|
|
DO 250 K2 = I, N
|
|
IF (WORK(IPVT(K1),JPVT(K2)) .LT. TEMP) THEN
|
|
TEMP = WORK(IPVT(K1),JPVT(K2))
|
|
IMAX = K1
|
|
JMAX = K2
|
|
END IF
|
|
250 CONTINUE
|
|
260 CONTINUE
|
|
C
|
|
C--update row and column pivots for indirect addressing of work
|
|
C
|
|
ITEMP = IPVT(I)
|
|
IPVT(I) = IPVT(IMAX)
|
|
IPVT(IMAX) = ITEMP
|
|
C
|
|
ITEMP = JPVT(I)
|
|
JPVT(I) = JPVT(JMAX)
|
|
JPVT(JMAX) = ITEMP
|
|
C
|
|
C--do next submatrix
|
|
C
|
|
270 CONTINUE
|
|
C
|
|
C--this threshold for determining when an eigenvalue has
|
|
C--not moved, and is therefore uncontrollable, is critical,
|
|
C--and may require future changes with more experience.
|
|
C
|
|
EPS1 = SQRT(EPS)
|
|
C
|
|
C--for each eigenvalue pair, decide if it is controllable
|
|
C
|
|
DO 280 I = 1, N
|
|
C
|
|
C--note that we are working with the "pivoted" work matrix
|
|
C--and are looking at its diagonal elements
|
|
C
|
|
IF (WORK(IPVT(I),JPVT(I))/ANORM .LE. EPS1) THEN
|
|
CON(I) = .FALSE.
|
|
ELSE
|
|
CON(I) = .TRUE.
|
|
END IF
|
|
280 CONTINUE
|
|
C
|
|
C--finally!
|
|
C
|
|
RETURN
|
|
END
|