diff options
Diffstat (limited to 'tests/examplefiles/ahcon.f')
| -rw-r--r-- | tests/examplefiles/ahcon.f | 340 |
1 files changed, 0 insertions, 340 deletions
diff --git a/tests/examplefiles/ahcon.f b/tests/examplefiles/ahcon.f deleted file mode 100644 index 48ae920b..00000000 --- a/tests/examplefiles/ahcon.f +++ /dev/null @@ -1,340 +0,0 @@ - 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 |