Mercurial > forge
changeset 9469:fccb049af652 octave-forge
control: reorganize package makefile (3)
| author | paramaniac |
|---|---|
| date | Wed, 22 Feb 2012 18:34:17 +0000 |
| parents | 08aac5184699 |
| children | 19411ae007a3 |
| files | main/control/devel/makefile_control.m main/control/src/Makefile main/control/src/TG04BX.f main/control/src/TG04BX.fortran |
| diffstat | 4 files changed, 264 insertions(+), 300 deletions(-) [+] |
line wrap: on
line diff
--- a/main/control/devel/makefile_control.m Wed Feb 22 18:24:10 2012 +0000 +++ b/main/control/devel/makefile_control.m Wed Feb 22 18:34:17 2012 +0000 @@ -30,7 +30,7 @@ makefile_zero %} -% system ("make clean"); +system ("make clean"); system ("make -j4 all"); system ("rm *.o");
--- a/main/control/src/Makefile Wed Feb 22 18:24:10 2012 +0000 +++ b/main/control/src/Makefile Wed Feb 22 18:34:17 2012 +0000 @@ -27,6 +27,7 @@ tar -xzf slicot.tar.gz mv slicot/src/*.f . mv slicot/src_aux/*.f . + cp TG04BX.fortran TG04BX.f mkoctfile *.f \ ${LAPACK_LIBS} ${BLAS_LIBS} ${FLIBS} ar -r slicotlibrary.a *.o @@ -35,83 +36,46 @@ # transmission zeros of state-space models slab08nd.oct: slab08nd.cc slicotlibrary.a mkoctfile slab08nd.cc slicotlibrary.a \ - AB08ND.f AB08NX.f TB01ID.f MB03OY.f MB03PY.f \ ${LAPACK_LIBS} ${BLAS_LIBS} # L-infinity norm slab13dd.oct: slab13dd.cc slicotlibrary.a mkoctfile slab13dd.cc slicotlibrary.a \ - AB13DD.f MA02AD.f MB01SD.f MB03XD.f TB01ID.f \ - TG01AD.f TG01BD.f AB13DX.f MA01AD.f MA02ID.f \ - MB03XP.f MB04DD.f MB04QB.f MB04TB.f MB03XU.f \ - MB04TS.f UE01MD.f MB02RD.f MB02SD.f MB04QC.f \ - MB04QF.f MB03YA.f MB03YD.f MB02RZ.f MB04QU.f \ - MB02SZ.f MB03YT.f \ ${LAPACK_LIBS} ${BLAS_LIBS} ${FLIBS} # H-2 controller synthesis - continuous-time slsb10hd.oct: slsb10hd.cc slicotlibrary.a mkoctfile slsb10hd.cc slicotlibrary.a \ - SB10HD.f SB10UD.f SB10VD.f SB10WD.f SB02RD.f \ - MB01RU.f SB02SD.f MA02ED.f SB02RU.f SB02MR.f \ - MB01SD.f SB02MS.f SB02MV.f SB02MW.f MA02AD.f \ - SB02QD.f MB02PD.f SB03QX.f SB03QY.f MB01RX.f \ - MB01RY.f SB03SX.f SB03SY.f select.f SB03MX.f \ - SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f \ ${LAPACK_LIBS} ${BLAS_LIBS} # H-2 controller synthesis - discrete-time slsb10ed.oct: slsb10ed.cc slicotlibrary.a mkoctfile slsb10ed.cc slicotlibrary.a \ - SB10ED.f SB10SD.f SB10TD.f SB10PD.f MB01RX.f \ - SB02SD.f SB02OD.f MB01RU.f SB02OU.f SB02OV.f \ - SB02OW.f MB01RY.f SB02OY.f SB03SX.f SB03SY.f \ - MA02ED.f select.f SB03MX.f SB02MR.f SB02MV.f \ - MB01UD.f SB03MV.f SB04PX.f \ ${LAPACK_LIBS} ${BLAS_LIBS} # H-2 norm slab13bd.oct: slab13bd.cc slicotlibrary.a mkoctfile slab13bd.cc slicotlibrary.a \ - AB13BD.f SB08DD.f SB03OU.f SB01FY.f TB01LD.f \ - SB03OT.f MB04ND.f MB04OD.f MB03QX.f select.f \ - SB03OR.f MB04OX.f MB03QD.f SB03OY.f MA02AD.f \ - MB03QY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \ ${LAPACK_LIBS} ${BLAS_LIBS} # Pole assignment slsb01bd.oct: slsb01bd.cc slicotlibrary.a mkoctfile slsb01bd.cc slicotlibrary.a \ - SB01BD.f MB03QD.f MB03QY.f SB01BX.f SB01BY.f \ - select.f \ ${LAPACK_LIBS} ${BLAS_LIBS} # H-inf controller synthesis - continuous-time slsb10fd.oct: slsb10fd.cc slicotlibrary.a mkoctfile slsb10fd.cc slicotlibrary.a \ - SB10FD.f SB10PD.f SB10QD.f SB10RD.f SB02RD.f \ - MB01RU.f MB01RX.f MA02AD.f SB02SD.f MA02ED.f \ - SB02RU.f SB02MR.f MB01SD.f SB02MS.f SB02MV.f \ - SB02MW.f SB02QD.f MB02PD.f SB03QX.f SB03QY.f \ - MB01RY.f SB03SX.f SB03SY.f select.f SB03MX.f \ - SB03MY.f MB01UD.f SB03MV.f SB03MW.f SB04PX.f \ ${LAPACK_LIBS} ${BLAS_LIBS} # H-inf controller synthesis - discrete-time slsb10dd.oct: slsb10dd.cc slicotlibrary.a mkoctfile slsb10dd.cc slicotlibrary.a \ - SB10DD.f MB01RU.f MB01RX.f SB02SD.f SB02OD.f \ - MA02AD.f SB02OU.f SB02OV.f SB02OW.f MB01RY.f \ - SB02OY.f SB03SX.f SB03SY.f MA02ED.f select.f \ - SB03MX.f SB02MR.f SB02MV.f MB01UD.f SB03MV.f \ - SB04PX.f \ ${LAPACK_LIBS} ${BLAS_LIBS} # Lyapunov equations slsb03md.oct: slsb03md.cc slicotlibrary.a mkoctfile slsb03md.cc slicotlibrary.a \ - SB03MD.f select.f SB03MX.f SB03MY.f MB01RD.f \ - SB03MV.f SB03MW.f SB04PX.f \ ${LAPACK_LIBS} ${BLAS_LIBS} # Sylvester equations - continuous-time @@ -256,4 +220,4 @@ rm -rf *.o core octave-core *.oct *~ *.f slicot realclean: clean - rm -rf *.a \ No newline at end of file + rm -rf *.a
--- a/main/control/src/TG04BX.f Wed Feb 22 18:24:10 2012 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,261 +0,0 @@ - SUBROUTINE TG04BX( IP, IZ, A, LDA, E, B, C, D, PR, PI, ZR, ZI, - $ GAIN, IWORK ) -C -C WARNING -C -C This routine is a modified version of TB04BX. It is intended -C for the Octave Control Systems Package and supports Descriptor -C State-Space models. TG04BX is *NOT* part of SLICOT and the -C authors from NICONET e.V. are *NOT* responsible for it. -C See file DESCRIPTION for the current maintainer of the Octave -C control package. -C -C SLICOT RELEASE 5.0. -C -C Copyright (c) 2002-2009 NICONET e.V. -C -C This program is free software: you can redistribute it and/or -C modify it under the terms of the GNU General Public License as -C published by the Free Software Foundation, either version 2 of -C the License, or (at your option) any later version. -C -C This program is distributed in the hope that it will be useful, -C but WITHOUT ANY WARRANTY; without even the implied warranty of -C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -C GNU General Public License for more details. -C -C You should have received a copy of the GNU General Public License -C along with this program. If not, see -C <http://www.gnu.org/licenses/>. -C -C PURPOSE -C -C To compute the gain of a single-input single-output linear system, -C given its state-space representation (A,b,c,d), and its poles and -C zeros. The matrix A is assumed to be in an upper Hessenberg form. -C The gain is computed using the formula -C -C -1 IP IZ -C g = (c*( S0*E - A ) *b + d)*Prod( S0 - Pi )/Prod( S0 - Zi ) , -C i=1 i=1 (1) -C -C where Pi, i = 1 : IP, and Zj, j = 1 : IZ, are the poles and zeros, -C respectively, and S0 is a real scalar different from all poles and -C zeros. -C -C ARGUMENTS -C -C Input/Output Parameters -C -C IP (input) INTEGER -C The number of the system poles. IP >= 0. -C -C IZ (input) INTEGER -C The number of the system zeros. IZ >= 0. -C -C A (input/output) DOUBLE PRECISION array, dimension (LDA,IP) -C On entry, the leading IP-by-IP part of this array must -C contain the state dynamics matrix A in an upper Hessenberg -C form. The elements below the second diagonal are not -C referenced. -C On exit, the leading IP-by-IP upper Hessenberg part of -C this array contains the LU factorization of the matrix -C A - S0*I, as computed by SLICOT Library routine MB02SD. -C -C LDA INTEGER -C The leading dimension of array A. LDA >= max(1,IP). -C -C B (input/output) DOUBLE PRECISION array, dimension (IP) -C On entry, this array must contain the system input -C vector b. -C On exit, this array contains the solution of the linear -C system ( A - S0*I )x = b . -C -C C (input) DOUBLE PRECISION array, dimension (IP) -C This array must contain the system output vector c. -C -C D (input) DOUBLE PRECISION -C The variable must contain the system feedthrough scalar d. -C -C PR (input) DOUBLE PRECISION array, dimension (IP) -C This array must contain the real parts of the system -C poles. Pairs of complex conjugate poles must be stored in -C consecutive memory locations. -C -C PI (input) DOUBLE PRECISION array, dimension (IP) -C This array must contain the imaginary parts of the system -C poles. -C -C ZR (input) DOUBLE PRECISION array, dimension (IZ) -C This array must contain the real parts of the system -C zeros. Pairs of complex conjugate zeros must be stored in -C consecutive memory locations. -C -C ZI (input) DOUBLE PRECISION array, dimension (IZ) -C This array must contain the imaginary parts of the system -C zeros. -C -C GAIN (output) DOUBLE PRECISION -C The gain of the linear system (A,b,c,d), given by (1). -C -C Workspace -C -C IWORK INTEGER array, dimension (IP) -C On exit, it contains the pivot indices; for 1 <= i <= IP, -C row i of the matrix A - S0*I was interchanged with -C row IWORK(i). -C -C METHOD -C -C The routine implements the method presented in [1]. A suitable -C value of S0 is chosen based on the system poles and zeros. -C Then, the LU factorization of the upper Hessenberg, nonsingular -C matrix A - S0*I is computed and used to solve the linear system -C in (1). -C -C REFERENCES -C -C [1] Varga, A. and Sima, V. -C Numerically Stable Algorithm for Transfer Function Matrix -C Evaluation. -C Int. J. Control, vol. 33, nr. 6, pp. 1123-1133, 1981. -C -C NUMERICAL ASPECTS -C -C The algorithm is numerically stable in practice and requires -C O(IP*IP) floating point operations. -C -C CONTRIBUTORS -C -C V. Sima, Research Institute for Informatics, Bucharest, May 2002. -C Partly based on the BIMASC Library routine GAIN by A. Varga. -C -C REVISIONS -C -C 2011-02-08 (Lukas Reichlin) Modifications for Descriptor Systems. -C -C KEYWORDS -C -C Eigenvalue, state-space representation, transfer function, zeros. -C -C ****************************************************************** -C -C .. Parameters .. - DOUBLE PRECISION ZERO, ONE, TWO, P1, ONEP1 - PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, - $ P1 = 0.1D0, ONEP1 = 1.1D0 ) -C .. Scalar Arguments .. - DOUBLE PRECISION D, GAIN - INTEGER IP, IZ, LDA -C .. Array Arguments .. - DOUBLE PRECISION A(LDA,*), E(LDA,*), B(*), C(*), PI(*), PR(*), - $ ZI(*), ZR(*) - INTEGER IWORK(*) -C .. Local Scalars .. - INTEGER I, J, INFO - DOUBLE PRECISION S0, S -C .. External Functions .. - DOUBLE PRECISION DDOT - EXTERNAL DDOT -C .. External Subroutines .. - EXTERNAL MB02RD, MB02SD -C .. Intrinsic Functions .. - INTRINSIC ABS, MAX -C .. -C .. Executable Statements .. -C -C For efficiency, the input scalar parameters are not checked. -C -C Quick return if possible. -C -C IF( IP.EQ.0 ) THEN -C GAIN = ZERO -C RETURN -C END IF -C -C Compute a suitable value for S0 . -C - S0 = ZERO -C - DO 10 I = 1, IP - S = ABS( PR(I) ) - IF ( PI(I).NE.ZERO ) - $ S = S + ABS( PI(I) ) - S0 = MAX( S0, S ) - 10 CONTINUE -C - DO 20 I = 1, IZ - S = ABS( ZR(I) ) - IF ( ZI(I).NE.ZERO ) - $ S = S + ABS( ZI(I) ) - S0 = MAX( S0, S ) - 20 CONTINUE -C - S0 = TWO*S0 + P1 - IF ( S0.LE.ONE ) - $ S0 = ONEP1 -C -C Form A - S0*E . -C - DO 30 J = 1, LDA - DO 25 I = 1, LDA - A(I,J) = A(I,J) - S0*E(I,J) - 25 CONTINUE - 30 CONTINUE -C -C Compute the LU factorization of the matrix A - S0*E -C (guaranteed to be nonsingular). -C -C CALL MB02SD( IP, A, LDA, IWORK, INFO ) - CALL MB02SD( LDA, A, LDA, IWORK, INFO ) -C -C Solve the linear system (A - S0*E)*x = b . -C -C CALL MB02RD( 'No Transpose', IP, 1, A, LDA, IWORK, B, IP, INFO ) -C CALL MB02RD( 'No Transpose', IP, 1, A, LDA, IWORK, B, LDA, INFO ) - CALL MB02RD( 'No Transpose', LDA, 1, A, LDA, IWORK, B, LDA, INFO ) -C -1 -C Compute c*(S0*E - A) *b + d . -C -C GAIN = D - DDOT( IP, C, 1, B, 1 ) - GAIN = D - DDOT( LDA, C, 1, B, 1 ) -C -C Multiply by the products in terms of poles and zeros in (1). -C - I = 1 -C -C WHILE ( I <= IP ) DO -C - 40 IF ( I.LE.IP ) THEN - IF ( PI(I).EQ.ZERO ) THEN - GAIN = GAIN*( S0 - PR(I) ) - I = I + 1 - ELSE - GAIN = GAIN*( S0*( S0 - TWO*PR(I) ) + PR(I)**2 + PI(I)**2 ) - I = I + 2 - END IF - GO TO 40 - END IF -C -C END WHILE 40 -C - I = 1 -C -C WHILE ( I <= IZ ) DO -C - 50 IF ( I.LE.IZ ) THEN - IF ( ZI(I).EQ.ZERO ) THEN - GAIN = GAIN/( S0 - ZR(I) ) - I = I + 1 - ELSE - GAIN = GAIN/( S0*( S0 - TWO*ZR(I) ) + ZR(I)**2 + ZI(I)**2 ) - I = I + 2 - END IF - GO TO 50 - END IF -C -C END WHILE 50 -C - RETURN -C *** Last line of TG04BX *** - END
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/main/control/src/TG04BX.fortran Wed Feb 22 18:34:17 2012 +0000 @@ -0,0 +1,261 @@ + SUBROUTINE TG04BX( IP, IZ, A, LDA, E, B, C, D, PR, PI, ZR, ZI, + $ GAIN, IWORK ) +C +C WARNING +C +C This routine is a modified version of TB04BX. It is intended +C for the Octave Control Systems Package and supports Descriptor +C State-Space models. TG04BX is *NOT* part of SLICOT and the +C authors from NICONET e.V. are *NOT* responsible for it. +C See file DESCRIPTION for the current maintainer of the Octave +C control package. +C +C SLICOT RELEASE 5.0. +C +C Copyright (c) 2002-2009 NICONET e.V. +C +C This program is free software: you can redistribute it and/or +C modify it under the terms of the GNU General Public License as +C published by the Free Software Foundation, either version 2 of +C the License, or (at your option) any later version. +C +C This program is distributed in the hope that it will be useful, +C but WITHOUT ANY WARRANTY; without even the implied warranty of +C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +C GNU General Public License for more details. +C +C You should have received a copy of the GNU General Public License +C along with this program. If not, see +C <http://www.gnu.org/licenses/>. +C +C PURPOSE +C +C To compute the gain of a single-input single-output linear system, +C given its state-space representation (A,b,c,d), and its poles and +C zeros. The matrix A is assumed to be in an upper Hessenberg form. +C The gain is computed using the formula +C +C -1 IP IZ +C g = (c*( S0*E - A ) *b + d)*Prod( S0 - Pi )/Prod( S0 - Zi ) , +C i=1 i=1 (1) +C +C where Pi, i = 1 : IP, and Zj, j = 1 : IZ, are the poles and zeros, +C respectively, and S0 is a real scalar different from all poles and +C zeros. +C +C ARGUMENTS +C +C Input/Output Parameters +C +C IP (input) INTEGER +C The number of the system poles. IP >= 0. +C +C IZ (input) INTEGER +C The number of the system zeros. IZ >= 0. +C +C A (input/output) DOUBLE PRECISION array, dimension (LDA,IP) +C On entry, the leading IP-by-IP part of this array must +C contain the state dynamics matrix A in an upper Hessenberg +C form. The elements below the second diagonal are not +C referenced. +C On exit, the leading IP-by-IP upper Hessenberg part of +C this array contains the LU factorization of the matrix +C A - S0*I, as computed by SLICOT Library routine MB02SD. +C +C LDA INTEGER +C The leading dimension of array A. LDA >= max(1,IP). +C +C B (input/output) DOUBLE PRECISION array, dimension (IP) +C On entry, this array must contain the system input +C vector b. +C On exit, this array contains the solution of the linear +C system ( A - S0*I )x = b . +C +C C (input) DOUBLE PRECISION array, dimension (IP) +C This array must contain the system output vector c. +C +C D (input) DOUBLE PRECISION +C The variable must contain the system feedthrough scalar d. +C +C PR (input) DOUBLE PRECISION array, dimension (IP) +C This array must contain the real parts of the system +C poles. Pairs of complex conjugate poles must be stored in +C consecutive memory locations. +C +C PI (input) DOUBLE PRECISION array, dimension (IP) +C This array must contain the imaginary parts of the system +C poles. +C +C ZR (input) DOUBLE PRECISION array, dimension (IZ) +C This array must contain the real parts of the system +C zeros. Pairs of complex conjugate zeros must be stored in +C consecutive memory locations. +C +C ZI (input) DOUBLE PRECISION array, dimension (IZ) +C This array must contain the imaginary parts of the system +C zeros. +C +C GAIN (output) DOUBLE PRECISION +C The gain of the linear system (A,b,c,d), given by (1). +C +C Workspace +C +C IWORK INTEGER array, dimension (IP) +C On exit, it contains the pivot indices; for 1 <= i <= IP, +C row i of the matrix A - S0*I was interchanged with +C row IWORK(i). +C +C METHOD +C +C The routine implements the method presented in [1]. A suitable +C value of S0 is chosen based on the system poles and zeros. +C Then, the LU factorization of the upper Hessenberg, nonsingular +C matrix A - S0*I is computed and used to solve the linear system +C in (1). +C +C REFERENCES +C +C [1] Varga, A. and Sima, V. +C Numerically Stable Algorithm for Transfer Function Matrix +C Evaluation. +C Int. J. Control, vol. 33, nr. 6, pp. 1123-1133, 1981. +C +C NUMERICAL ASPECTS +C +C The algorithm is numerically stable in practice and requires +C O(IP*IP) floating point operations. +C +C CONTRIBUTORS +C +C V. Sima, Research Institute for Informatics, Bucharest, May 2002. +C Partly based on the BIMASC Library routine GAIN by A. Varga. +C +C REVISIONS +C +C 2011-02-08 (Lukas Reichlin) Modifications for Descriptor Systems. +C +C KEYWORDS +C +C Eigenvalue, state-space representation, transfer function, zeros. +C +C ****************************************************************** +C +C .. Parameters .. + DOUBLE PRECISION ZERO, ONE, TWO, P1, ONEP1 + PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, + $ P1 = 0.1D0, ONEP1 = 1.1D0 ) +C .. Scalar Arguments .. + DOUBLE PRECISION D, GAIN + INTEGER IP, IZ, LDA +C .. Array Arguments .. + DOUBLE PRECISION A(LDA,*), E(LDA,*), B(*), C(*), PI(*), PR(*), + $ ZI(*), ZR(*) + INTEGER IWORK(*) +C .. Local Scalars .. + INTEGER I, J, INFO + DOUBLE PRECISION S0, S +C .. External Functions .. + DOUBLE PRECISION DDOT + EXTERNAL DDOT +C .. External Subroutines .. + EXTERNAL MB02RD, MB02SD +C .. Intrinsic Functions .. + INTRINSIC ABS, MAX +C .. +C .. Executable Statements .. +C +C For efficiency, the input scalar parameters are not checked. +C +C Quick return if possible. +C +C IF( IP.EQ.0 ) THEN +C GAIN = ZERO +C RETURN +C END IF +C +C Compute a suitable value for S0 . +C + S0 = ZERO +C + DO 10 I = 1, IP + S = ABS( PR(I) ) + IF ( PI(I).NE.ZERO ) + $ S = S + ABS( PI(I) ) + S0 = MAX( S0, S ) + 10 CONTINUE +C + DO 20 I = 1, IZ + S = ABS( ZR(I) ) + IF ( ZI(I).NE.ZERO ) + $ S = S + ABS( ZI(I) ) + S0 = MAX( S0, S ) + 20 CONTINUE +C + S0 = TWO*S0 + P1 + IF ( S0.LE.ONE ) + $ S0 = ONEP1 +C +C Form A - S0*E . +C + DO 30 J = 1, LDA + DO 25 I = 1, LDA + A(I,J) = A(I,J) - S0*E(I,J) + 25 CONTINUE + 30 CONTINUE +C +C Compute the LU factorization of the matrix A - S0*E +C (guaranteed to be nonsingular). +C +C CALL MB02SD( IP, A, LDA, IWORK, INFO ) + CALL MB02SD( LDA, A, LDA, IWORK, INFO ) +C +C Solve the linear system (A - S0*E)*x = b . +C +C CALL MB02RD( 'No Transpose', IP, 1, A, LDA, IWORK, B, IP, INFO ) +C CALL MB02RD( 'No Transpose', IP, 1, A, LDA, IWORK, B, LDA, INFO ) + CALL MB02RD( 'No Transpose', LDA, 1, A, LDA, IWORK, B, LDA, INFO ) +C -1 +C Compute c*(S0*E - A) *b + d . +C +C GAIN = D - DDOT( IP, C, 1, B, 1 ) + GAIN = D - DDOT( LDA, C, 1, B, 1 ) +C +C Multiply by the products in terms of poles and zeros in (1). +C + I = 1 +C +C WHILE ( I <= IP ) DO +C + 40 IF ( I.LE.IP ) THEN + IF ( PI(I).EQ.ZERO ) THEN + GAIN = GAIN*( S0 - PR(I) ) + I = I + 1 + ELSE + GAIN = GAIN*( S0*( S0 - TWO*PR(I) ) + PR(I)**2 + PI(I)**2 ) + I = I + 2 + END IF + GO TO 40 + END IF +C +C END WHILE 40 +C + I = 1 +C +C WHILE ( I <= IZ ) DO +C + 50 IF ( I.LE.IZ ) THEN + IF ( ZI(I).EQ.ZERO ) THEN + GAIN = GAIN/( S0 - ZR(I) ) + I = I + 1 + ELSE + GAIN = GAIN/( S0*( S0 - TWO*ZR(I) ) + ZR(I)**2 + ZI(I)**2 ) + I = I + 2 + END IF + GO TO 50 + END IF +C +C END WHILE 50 +C + RETURN +C *** Last line of TG04BX *** + END