SUBROUTINE AHCON (SIZE,N,M,A,B,OLEVR,OLEVI,CLEVR,CLEVI,
& SCR1,SCR2,IPVT,JPVT,CON,WORK,ISEED,IERR)
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