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