Revision: 9649
http://octave.svn.sourceforge.net/octave/?rev=9649&view=rev
Author: paramaniac
Date: 2012-02-22 16:37:14 +0000 (Wed, 22 Feb 2012)
Log Message:
-----------
control-devel: reorganize makefile
Modified Paths:
--------------
trunk/octave-forge/extra/control-devel/devel/makefile_modred.m
trunk/octave-forge/extra/control-devel/src/Makefile
Added Paths:
-----------
trunk/octave-forge/extra/control-devel/src/slicot.tar.gz
Removed Paths:
-------------
trunk/octave-forge/extra/control-devel/src/AB04MD.f
trunk/octave-forge/extra/control-devel/src/AB05PD.f
trunk/octave-forge/extra/control-devel/src/AB05QD.f
trunk/octave-forge/extra/control-devel/src/AB07MD.f
trunk/octave-forge/extra/control-devel/src/AB07ND.f
trunk/octave-forge/extra/control-devel/src/AB08MD.f
trunk/octave-forge/extra/control-devel/src/AB08NX.f
trunk/octave-forge/extra/control-devel/src/AB09AD.f
trunk/octave-forge/extra/control-devel/src/AB09AX.f
trunk/octave-forge/extra/control-devel/src/AB09BD.f
trunk/octave-forge/extra/control-devel/src/AB09BX.f
trunk/octave-forge/extra/control-devel/src/AB09CX.f
trunk/octave-forge/extra/control-devel/src/AB09DD.f
trunk/octave-forge/extra/control-devel/src/AB09HD.f
trunk/octave-forge/extra/control-devel/src/AB09HY.f
trunk/octave-forge/extra/control-devel/src/AB09ID.f
trunk/octave-forge/extra/control-devel/src/AB09IX.f
trunk/octave-forge/extra/control-devel/src/AB09IY.f
trunk/octave-forge/extra/control-devel/src/AB09JD.f
trunk/octave-forge/extra/control-devel/src/AB09JV.f
trunk/octave-forge/extra/control-devel/src/AB09JW.f
trunk/octave-forge/extra/control-devel/src/AB09JX.f
trunk/octave-forge/extra/control-devel/src/AG07BD.f
trunk/octave-forge/extra/control-devel/src/DG01MD.f
trunk/octave-forge/extra/control-devel/src/IB01AD.f
trunk/octave-forge/extra/control-devel/src/IB01BD.f
trunk/octave-forge/extra/control-devel/src/IB01CD.f
trunk/octave-forge/extra/control-devel/src/IB01MD.f
trunk/octave-forge/extra/control-devel/src/IB01MY.f
trunk/octave-forge/extra/control-devel/src/IB01ND.f
trunk/octave-forge/extra/control-devel/src/IB01OD.f
trunk/octave-forge/extra/control-devel/src/IB01OY.f
trunk/octave-forge/extra/control-devel/src/IB01PD.f
trunk/octave-forge/extra/control-devel/src/IB01PX.f
trunk/octave-forge/extra/control-devel/src/IB01PY.f
trunk/octave-forge/extra/control-devel/src/IB01QD.f
trunk/octave-forge/extra/control-devel/src/IB01RD.f
trunk/octave-forge/extra/control-devel/src/MA02AD.f
trunk/octave-forge/extra/control-devel/src/MA02BD.f
trunk/octave-forge/extra/control-devel/src/MA02DD.f
trunk/octave-forge/extra/control-devel/src/MA02ED.f
trunk/octave-forge/extra/control-devel/src/MA02FD.f
trunk/octave-forge/extra/control-devel/src/MA02GD.f
trunk/octave-forge/extra/control-devel/src/MB01PD.f
trunk/octave-forge/extra/control-devel/src/MB01QD.f
trunk/octave-forge/extra/control-devel/src/MB01RU.f
trunk/octave-forge/extra/control-devel/src/MB01RX.f
trunk/octave-forge/extra/control-devel/src/MB01RY.f
trunk/octave-forge/extra/control-devel/src/MB01SD.f
trunk/octave-forge/extra/control-devel/src/MB01TD.f
trunk/octave-forge/extra/control-devel/src/MB01UD.f
trunk/octave-forge/extra/control-devel/src/MB01VD.f
trunk/octave-forge/extra/control-devel/src/MB01WD.f
trunk/octave-forge/extra/control-devel/src/MB01YD.f
trunk/octave-forge/extra/control-devel/src/MB01ZD.f
trunk/octave-forge/extra/control-devel/src/MB02PD.f
trunk/octave-forge/extra/control-devel/src/MB02QY.f
trunk/octave-forge/extra/control-devel/src/MB02UD.f
trunk/octave-forge/extra/control-devel/src/MB03OD.f
trunk/octave-forge/extra/control-devel/src/MB03OY.f
trunk/octave-forge/extra/control-devel/src/MB03PY.f
trunk/octave-forge/extra/control-devel/src/MB03QD.f
trunk/octave-forge/extra/control-devel/src/MB03QX.f
trunk/octave-forge/extra/control-devel/src/MB03QY.f
trunk/octave-forge/extra/control-devel/src/MB03UD.f
trunk/octave-forge/extra/control-devel/src/MB04ID.f
trunk/octave-forge/extra/control-devel/src/MB04IY.f
trunk/octave-forge/extra/control-devel/src/MB04KD.f
trunk/octave-forge/extra/control-devel/src/MB04ND.f
trunk/octave-forge/extra/control-devel/src/MB04NY.f
trunk/octave-forge/extra/control-devel/src/MB04OD.f
trunk/octave-forge/extra/control-devel/src/MB04OX.f
trunk/octave-forge/extra/control-devel/src/MB04OY.f
trunk/octave-forge/extra/control-devel/src/MC01PD.f
trunk/octave-forge/extra/control-devel/src/SB01FY.f
trunk/octave-forge/extra/control-devel/src/SB02MD.f
trunk/octave-forge/extra/control-devel/src/SB02MR.f
trunk/octave-forge/extra/control-devel/src/SB02MS.f
trunk/octave-forge/extra/control-devel/src/SB02MT.f
trunk/octave-forge/extra/control-devel/src/SB02MU.f
trunk/octave-forge/extra/control-devel/src/SB02MV.f
trunk/octave-forge/extra/control-devel/src/SB02MW.f
trunk/octave-forge/extra/control-devel/src/SB02ND.f
trunk/octave-forge/extra/control-devel/src/SB02QD.f
trunk/octave-forge/extra/control-devel/src/SB02RD.f
trunk/octave-forge/extra/control-devel/src/SB02RU.f
trunk/octave-forge/extra/control-devel/src/SB02SD.f
trunk/octave-forge/extra/control-devel/src/SB03MV.f
trunk/octave-forge/extra/control-devel/src/SB03MW.f
trunk/octave-forge/extra/control-devel/src/SB03MX.f
trunk/octave-forge/extra/control-devel/src/SB03MY.f
trunk/octave-forge/extra/control-devel/src/SB03OD.f
trunk/octave-forge/extra/control-devel/src/SB03OR.f
trunk/octave-forge/extra/control-devel/src/SB03OT.f
trunk/octave-forge/extra/control-devel/src/SB03OU.f
trunk/octave-forge/extra/control-devel/src/SB03OV.f
trunk/octave-forge/extra/control-devel/src/SB03OY.f
trunk/octave-forge/extra/control-devel/src/SB03QX.f
trunk/octave-forge/extra/control-devel/src/SB03QY.f
trunk/octave-forge/extra/control-devel/src/SB03SX.f
trunk/octave-forge/extra/control-devel/src/SB03SY.f
trunk/octave-forge/extra/control-devel/src/SB04PX.f
trunk/octave-forge/extra/control-devel/src/SB04PY.f
trunk/octave-forge/extra/control-devel/src/SB08CD.f
trunk/octave-forge/extra/control-devel/src/SB08DD.f
trunk/octave-forge/extra/control-devel/src/SB08GD.f
trunk/octave-forge/extra/control-devel/src/SB08HD.f
trunk/octave-forge/extra/control-devel/src/SB10YD.f
trunk/octave-forge/extra/control-devel/src/SB10ZP.f
trunk/octave-forge/extra/control-devel/src/SB16AD.f
trunk/octave-forge/extra/control-devel/src/SB16AY.f
trunk/octave-forge/extra/control-devel/src/SB16BD.f
trunk/octave-forge/extra/control-devel/src/SB16CD.f
trunk/octave-forge/extra/control-devel/src/SB16CY.f
trunk/octave-forge/extra/control-devel/src/TB01ID.f
trunk/octave-forge/extra/control-devel/src/TB01KD.f
trunk/octave-forge/extra/control-devel/src/TB01LD.f
trunk/octave-forge/extra/control-devel/src/TB01PD.f
trunk/octave-forge/extra/control-devel/src/TB01UD.f
trunk/octave-forge/extra/control-devel/src/TB01WD.f
trunk/octave-forge/extra/control-devel/src/TB01XD.f
trunk/octave-forge/extra/control-devel/src/TD03AY.f
trunk/octave-forge/extra/control-devel/src/TD04AD.f
trunk/octave-forge/extra/control-devel/src/delctg.f
trunk/octave-forge/extra/control-devel/src/select.f
Modified: trunk/octave-forge/extra/control-devel/devel/makefile_modred.m
===================================================================
--- trunk/octave-forge/extra/control-devel/devel/makefile_modred.m 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/devel/makefile_modred.m 2012-02-22 16:37:14 UTC (rev 9649)
@@ -13,36 +13,17 @@
cd (srcdir);
mkoctfile slab09hd.cc \
- AB09HD.f TB01ID.f AB04MD.f TB01KD.f AB09HY.f \
- AB09IX.f MB03UD.f SB02MD.f AB09DD.f TB01LD.f \
- SB03OU.f MA02AD.f MB03QX.f select.f SB03OT.f \
- SB02MR.f SB02MS.f MB03QD.f SB02MU.f SB02MV.f \
- SB02MW.f MB04ND.f MB04OD.f MB03QY.f SB03OR.f \
- SB03OY.f SB04PX.f MB04NY.f MB04OY.f SB03OV.f \
+ slicotlibrary.a \
"$(mkoctfile -p LAPACK_LIBS)" \
"$(mkoctfile -p BLAS_LIBS)"
mkoctfile slab09id.cc \
- AB09ID.f TB01PD.f SB08DD.f TB01ID.f TB01KD.f \
- AB09IX.f AB09IY.f SB08CD.f MB04ND.f TB01XD.f \
- MB04OD.f MB01WD.f MB03UD.f AB07MD.f SB01FY.f \
- AB09DD.f TB01LD.f SB03OU.f TB01UD.f MA02AD.f \
- MA02BD.f MB03OY.f MB03QX.f MB01PD.f select.f \
- MB01YD.f MB04NY.f MB01ZD.f SB03OT.f MB04OX.f \
- MB04OY.f MB03QD.f SB03OY.f MB03QY.f MB01QD.f \
- SB03OR.f SB03OV.f SB04PX.f \
+ slicotlibrary.a \
"$(mkoctfile -p LAPACK_LIBS)" \
"$(mkoctfile -p BLAS_LIBS)"
mkoctfile slab09jd.cc \
- AB09JD.f TB01ID.f TB01KD.f AB07ND.f AB09JV.f \
- AB09JW.f AB09CX.f AG07BD.f AB08MD.f AB04MD.f \
- TB01LD.f delctg.f SB04PY.f AB09AX.f AB08NX.f \
- MB01SD.f AB09JX.f MA02AD.f TB01WD.f MB03OY.f \
- MB03PY.f MA02DD.f MB03UD.f MB03QX.f select.f \
- SB04PX.f SB03OU.f MB03QD.f MB03QY.f SB03OT.f \
- MB04ND.f MB04OD.f SB03OR.f SB03OY.f MB04NY.f \
- MB04OY.f SB03OV.f \
+ slicotlibrary.a \
"$(mkoctfile -p LAPACK_LIBS)" \
"$(mkoctfile -p BLAS_LIBS)"
Deleted: trunk/octave-forge/extra/control-devel/src/AB04MD.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB04MD.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB04MD.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,345 +0,0 @@
- SUBROUTINE AB04MD( TYPE, N, M, P, ALPHA, BETA, A, LDA, B, LDB, C,
- $ LDC, D, LDD, IWORK, DWORK, LDWORK, INFO )
-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 perform a transformation on the parameters (A,B,C,D) of a
-C system, which is equivalent to a bilinear transformation of the
-C corresponding transfer function matrix.
-C
-C ARGUMENTS
-C
-C Mode Parameters
-C
-C TYPE CHARACTER*1
-C Indicates the type of the original system and the
-C transformation to be performed as follows:
-C = 'D': discrete-time -> continuous-time;
-C = 'C': continuous-time -> discrete-time.
-C
-C Input/Output Parameters
-C
-C N (input) INTEGER
-C The order of the state matrix A. N >= 0.
-C
-C M (input) INTEGER
-C The number of system inputs. M >= 0.
-C
-C P (input) INTEGER
-C The number of system outputs. P >= 0.
-C
-C ALPHA, (input) DOUBLE PRECISION
-C BETA Parameters specifying the bilinear transformation.
-C Recommended values for stable systems: ALPHA = 1,
-C BETA = 1. ALPHA <> 0, BETA <> 0.
-C
-C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C On entry, the leading N-by-N part of this array must
-C contain the state matrix A of the original system.
-C On exit, the leading N-by-N part of this array contains
-C _
-C the state matrix A of the transformed system.
-C
-C LDA INTEGER
-C The leading dimension of array A. LDA >= MAX(1,N).
-C
-C B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
-C On entry, the leading N-by-M part of this array must
-C contain the input matrix B of the original system.
-C On exit, the leading N-by-M part of this array contains
-C _
-C the input matrix B of the transformed system.
-C
-C LDB INTEGER
-C The leading dimension of array B. LDB >= MAX(1,N).
-C
-C C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C On entry, the leading P-by-N part of this array must
-C contain the output matrix C of the original system.
-C On exit, the leading P-by-N part of this array contains
-C _
-C the output matrix C of the transformed system.
-C
-C LDC INTEGER
-C The leading dimension of array C. LDC >= MAX(1,P).
-C
-C D (input/output) DOUBLE PRECISION array, dimension (LDD,M)
-C On entry, the leading P-by-M part of this array must
-C contain the input/output matrix D for the original system.
-C On exit, the leading P-by-M part of this array contains
-C _
-C the input/output matrix D of the transformed system.
-C
-C LDD INTEGER
-C The leading dimension of array D. LDD >= MAX(1,P).
-C
-C Workspace
-C
-C IWORK INTEGER array, dimension (N)
-C
-C DWORK DOUBLE PRECISION array, dimension (LDWORK)
-C On exit, if INFO = 0, DWORK(1) returns the optimal value
-C of LDWORK.
-C
-C LDWORK INTEGER
-C The length of the array DWORK. LDWORK >= MAX(1,N).
-C For optimum performance LDWORK >= MAX(1,N*NB), where NB
-C is the optimal blocksize.
-C
-C Error Indicator
-C
-C INFO INTEGER
-C = 0: successful exit;
-C < 0: if INFO = -i, the i-th argument had an illegal
-C value;
-C = 1: if the matrix (ALPHA*I + A) is exactly singular;
-C = 2: if the matrix (BETA*I - A) is exactly singular.
-C
-C METHOD
-C
-C The parameters of the discrete-time system are transformed into
-C the parameters of the continuous-time system (TYPE = 'D'), or
-C vice-versa (TYPE = 'C') by the transformation:
-C
-C 1. Discrete -> continuous
-C _ -1
-C A = beta*(alpha*I + A) * (A - alpha*I)
-C _ -1
-C B = sqrt(2*alpha*beta) * (alpha*I + A) * B
-C _ -1
-C C = sqrt(2*alpha*beta) * C * (alpha*I + A)
-C _ -1
-C D = D - C * (alpha*I + A) * B
-C
-C which is equivalent to the bilinear transformation
-C
-C z - alpha
-C z -> s = beta --------- .
-C z + alpha
-C
-C of one transfer matrix onto the other.
-C
-C 2. Continuous -> discrete
-C _ -1
-C A = alpha*(beta*I - A) * (beta*I + A)
-C _ -1
-C B = sqrt(2*alpha*beta) * (beta*I - A) * B
-C _ -1
-C C = sqrt(2*alpha*beta) * C * (beta*I - A)
-C _ -1
-C D = D + C * (beta*I - A) * B
-C
-C which is equivalent to the bilinear transformation
-C
-C beta + s
-C s -> z = alpha -------- .
-C beta - s
-C
-C of one transfer matrix onto the other.
-C
-C REFERENCES
-C
-C [1] Al-Saggaf, U.M. and Franklin, G.F.
-C Model reduction via balanced realizations: a extension and
-C frequency weighting techniques.
-C IEEE Trans. Autom. Contr., AC-33, pp. 687-692, 1988.
-C
-C NUMERICAL ASPECTS
-C 3
-C The time taken is approximately proportional to N .
-C The accuracy depends mainly on the condition number of the matrix
-C to be inverted.
-C
-C CONTRIBUTORS
-C
-C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, and
-C A. Varga, German Aerospace Research Establishment,
-C Oberpfaffenhofen, Germany, Nov. 1996.
-C Supersedes Release 2.0 routine AB04AD by W. van der Linden, and
-C A.J. Geurts, Technische Hogeschool Eindhoven, Holland.
-C
-C REVISIONS
-C
-C -
-C
-C KEYWORDS
-C
-C Bilinear transformation, continuous-time system, discrete-time
-C system, state-space model.
-C
-C ******************************************************************
-C
-C .. Parameters ..
- DOUBLE PRECISION ZERO, ONE, TWO
- PARAMETER ( ZERO=0.0D0, ONE=1.0D0, TWO=2.0D0 )
-C .. Scalar Arguments ..
- CHARACTER TYPE
- INTEGER INFO, LDA, LDB, LDC, LDD, LDWORK, M, N, P
- DOUBLE PRECISION ALPHA, BETA
-C .. Array Arguments ..
- INTEGER IWORK(*)
- DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), DWORK(*)
-C .. Local Scalars ..
- LOGICAL LTYPE
- INTEGER I, IP
- DOUBLE PRECISION AB2, PALPHA, PBETA, SQRAB2
-C .. External Functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-C .. External Subroutines ..
- EXTERNAL DGEMM, DGETRF, DGETRS, DGETRI, DLASCL, DSCAL,
- $ DSWAP, XERBLA
-C .. Intrinsic Functions ..
- INTRINSIC ABS, MAX, SIGN, SQRT
-C .. Executable Statements ..
-C
- INFO = 0
- LTYPE = LSAME( TYPE, 'D' )
-C
-C Test the input scalar arguments.
-C
- IF( .NOT.LTYPE .AND. .NOT.LSAME( TYPE, 'C' ) ) THEN
- INFO = -1
- ELSE IF( N.LT.0 ) THEN
- INFO = -2
- ELSE IF( M.LT.0 ) THEN
- INFO = -3
- ELSE IF( P.LT.0 ) THEN
- INFO = -4
- ELSE IF( ALPHA.EQ.ZERO ) THEN
- INFO = -5
- ELSE IF( BETA.EQ.ZERO ) THEN
- INFO = -6
- ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
- INFO = -8
- ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
- INFO = -10
- ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
- INFO = -12
- ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
- INFO = -14
- ELSE IF( LDWORK.LT.MAX( 1, N ) ) THEN
- INFO = -17
- END IF
-C
- IF ( INFO.NE.0 ) THEN
-C
-C Error return.
-C
- CALL XERBLA( 'AB04MD', -INFO )
- RETURN
- END IF
-C
-C Quick return if possible.
-C
- IF ( MAX( N, M, P ).EQ.0 )
- $ RETURN
-C
-C (Note: Comments in the code beginning "Workspace:" describe the
-C minimal amount of real workspace needed at that point in the
-C code, as well as the preferred amount for good performance.
-C NB refers to the optimal block size for the immediately
-C following subroutine, as returned by ILAENV.)
-C
- IF (LTYPE) THEN
-C
-C Discrete-time to continuous-time with (ALPHA, BETA).
-C
- PALPHA = ALPHA
- PBETA = BETA
- ELSE
-C
-C Continuous-time to discrete-time with (ALPHA, BETA) is
-C equivalent with discrete-time to continuous-time with
-C (-BETA, -ALPHA), if B and C change the sign.
-C
- PALPHA = -BETA
- PBETA = -ALPHA
- END IF
-C
- AB2 = PALPHA*PBETA*TWO
- SQRAB2 = SIGN( SQRT( ABS( AB2 ) ), PALPHA )
-C -1
-C Compute (alpha*I + A) .
-C
- DO 10 I = 1, N
- A(I,I) = A(I,I) + PALPHA
- 10 CONTINUE
-C
- CALL DGETRF( N, N, A, LDA, IWORK, INFO )
-C
- IF (INFO.NE.0) THEN
-C
-C Error return.
-C
- IF (LTYPE) THEN
- INFO = 1
- ELSE
- INFO = 2
- END IF
- RETURN
- END IF
-C -1
-C Compute (alpha*I+A) *B.
-C
- CALL DGETRS( 'No transpose', N, M, A, LDA, IWORK, B, LDB, INFO )
-C -1
-C Compute D - C*(alpha*I+A) *B.
-C
- CALL DGEMM( 'No transpose', 'No transpose', P, M, N, -ONE, C,
- $ LDC, B, LDB, ONE, D, LDD )
-C
-C Scale B by sqrt(2*alpha*beta).
-C
- CALL DLASCL( 'General', 0, 0, ONE, SQRAB2, N, M, B, LDB, INFO )
-C -1
-C Compute sqrt(2*alpha*beta)*C*(alpha*I + A) .
-C
- CALL DTRSM( 'Right', 'Upper', 'No transpose', 'Non-unit', P, N,
- $ SQRAB2, A, LDA, C, LDC )
-C
- CALL DTRSM( 'Right', 'Lower', 'No transpose', 'Unit', P, N, ONE,
- $ A, LDA, C, LDC )
-C
-C Apply column interchanges to the solution matrix.
-C
- DO 20 I = N-1, 1, -1
- IP = IWORK(I)
- IF ( IP.NE.I )
- $ CALL DSWAP( P, C(1,I), 1, C(1,IP), 1 )
- 20 CONTINUE
-C -1
-C Compute beta*(alpha*I + A) *(A - alpha*I) as
-C -1
-C beta*I - 2*alpha*beta*(alpha*I + A) .
-C
-C Workspace: need N; prefer N*NB.
-C
- CALL DGETRI( N, A, LDA, IWORK, DWORK, LDWORK, INFO )
-C
- DO 30 I = 1, N
- CALL DSCAL(N, -AB2, A(1,I), 1)
- A(I,I) = A(I,I) + PBETA
- 30 CONTINUE
-C
- RETURN
-C *** Last line of AB04MD ***
- END
Deleted: trunk/octave-forge/extra/control-devel/src/AB05PD.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB05PD.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB05PD.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,385 +0,0 @@
- SUBROUTINE AB05PD( OVER, N1, M, P, N2, ALPHA, A1, LDA1, B1, LDB1,
- $ C1, LDC1, D1, LDD1, A2, LDA2, B2, LDB2, C2,
- $ LDC2, D2, LDD2, N, A, LDA, B, LDB, C, LDC, D,
- $ LDD, INFO)
-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 state-space model G = (A,B,C,D) corresponding to
-C the sum G = G1 + alpha*G2, where G1 = (A1,B1,C1,D1) and
-C G2 = (A2,B2,C2,D2). G, G1, and G2 are the transfer-function
-C matrices of the corresponding state-space models.
-C
-C ARGUMENTS
-C
-C Mode Parameters
-C
-C OVER CHARACTER*1
-C Indicates whether the user wishes to overlap pairs of
-C arrays, as follows:
-C = 'N': Do not overlap;
-C = 'O': Overlap pairs of arrays: A1 and A, B1 and B,
-C C1 and C, and D1 and D, i.e. the same name is
-C effectively used for each pair (for all pairs)
-C in the routine call. In this case, setting
-C LDA1 = LDA, LDB1 = LDB, LDC1 = LDC, and LDD1 = LDD
-C will give maximum efficiency.
-C
-C Input/Output Parameters
-C
-C N1 (input) INTEGER
-C The number of state variables in the first system, i.e.
-C the order of the matrix A1, the number of rows of B1 and
-C the number of columns of C1. N1 >= 0.
-C
-C M (input) INTEGER
-C The number of input variables of the two systems, i.e. the
-C number of columns of matrices B1, D1, B2 and D2. M >= 0.
-C
-C P (input) INTEGER
-C The number of output variables of the two systems, i.e.
-C the number of rows of matrices C1, D1, C2 and D2. P >= 0.
-C
-C N2 (input) INTEGER
-C The number of state variables in the second system, i.e.
-C the order of the matrix A2, the number of rows of B2 and
-C the number of columns of C2. N2 >= 0.
-C
-C ALPHA (input) DOUBLE PRECISION
-C The coefficient multiplying G2.
-C
-C A1 (input) DOUBLE PRECISION array, dimension (LDA1,N1)
-C The leading N1-by-N1 part of this array must contain the
-C state transition matrix A1 for the first system.
-C
-C LDA1 INTEGER
-C The leading dimension of array A1. LDA1 >= MAX(1,N1).
-C
-C B1 (input) DOUBLE PRECISION array, dimension (LDB1,M)
-C The leading N1-by-M part of this array must contain the
-C input/state matrix B1 for the first system.
-C
-C LDB1 INTEGER
-C The leading dimension of array B1. LDB1 >= MAX(1,N1).
-C
-C C1 (input) DOUBLE PRECISION array, dimension (LDC1,N1)
-C The leading P-by-N1 part of this array must contain the
-C state/output matrix C1 for the first system.
-C
-C LDC1 INTEGER
-C The leading dimension of array C1.
-C LDC1 >= MAX(1,P) if N1 > 0.
-C LDC1 >= 1 if N1 = 0.
-C
-C D1 (input) DOUBLE PRECISION array, dimension (LDD1,M)
-C The leading P-by-M part of this array must contain the
-C input/output matrix D1 for the first system.
-C
-C LDD1 INTEGER
-C The leading dimension of array D1. LDD1 >= MAX(1,P).
-C
-C A2 (input) DOUBLE PRECISION array, dimension (LDA2,N2)
-C The leading N2-by-N2 part of this array must contain the
-C state transition matrix A2 for the second system.
-C
-C LDA2 INTEGER
-C The leading dimension of array A2. LDA2 >= MAX(1,N2).
-C
-C B2 (input) DOUBLE PRECISION array, dimension (LDB2,M)
-C The leading N2-by-M part of this array must contain the
-C input/state matrix B2 for the second system.
-C
-C LDB2 INTEGER
-C The leading dimension of array B2. LDB2 >= MAX(1,N2).
-C
-C C2 (input) DOUBLE PRECISION array, dimension (LDC2,N2)
-C The leading P-by-N2 part of this array must contain the
-C state/output matrix C2 for the second system.
-C
-C LDC2 INTEGER
-C The leading dimension of array C2.
-C LDC2 >= MAX(1,P) if N2 > 0.
-C LDC2 >= 1 if N2 = 0.
-C
-C D2 (input) DOUBLE PRECISION array, dimension (LDD2,M)
-C The leading P-by-M part of this array must contain the
-C input/output matrix D2 for the second system.
-C
-C LDD2 INTEGER
-C The leading dimension of array D2. LDD2 >= MAX(1,P).
-C
-C N (output) INTEGER
-C The number of state variables (N1 + N2) in the resulting
-C system, i.e. the order of the matrix A, the number of rows
-C of B and the number of columns of C.
-C
-C A (output) DOUBLE PRECISION array, dimension (LDA,N1+N2)
-C The leading N-by-N part of this array contains the state
-C transition matrix A for the resulting system.
-C The array A can overlap A1 if OVER = 'O'.
-C
-C LDA INTEGER
-C The leading dimension of array A. LDA >= MAX(1,N1+N2).
-C
-C B (output) DOUBLE PRECISION array, dimension (LDB,M)
-C The leading N-by-M part of this array contains the
-C input/state matrix B for the resulting system.
-C The array B can overlap B1 if OVER = 'O'.
-C
-C LDB INTEGER
-C The leading dimension of array B. LDB >= MAX(1,N1+N2).
-C
-C C (output) DOUBLE PRECISION array, dimension (LDC,N1+N2)
-C The leading P-by-N part of this array contains the
-C state/output matrix C for the resulting system.
-C The array C can overlap C1 if OVER = 'O'.
-C
-C LDC INTEGER
-C The leading dimension of array C.
-C LDC >= MAX(1,P) if N1+N2 > 0.
-C LDC >= 1 if N1+N2 = 0.
-C
-C D (output) DOUBLE PRECISION array, dimension (LDD,M)
-C The leading P-by-M part of this array contains the
-C input/output matrix D for the resulting system.
-C The array D can overlap D1 if OVER = 'O'.
-C
-C LDD INTEGER
-C The leading dimension of array D. LDD >= MAX(1,P).
-C
-C Error Indicator
-C
-C INFO INTEGER
-C = 0: successful exit;
-C < 0: if INFO = -i, the i-th argument had an illegal
-C value.
-C
-C METHOD
-C
-C The matrices of the resulting systems are determined as:
-C
-C ( A1 0 ) ( B1 )
-C A = ( ) , B = ( ) ,
-C ( 0 A2 ) ( B2 )
-C
-C C = ( C1 alpha*C2 ) , D = D1 + alpha*D2 .
-C
-C REFERENCES
-C
-C None
-C
-C NUMERICAL ASPECTS
-C
-C None
-C
-C CONTRIBUTORS
-C
-C A. Varga, German Aerospace Research Establishment,
-C Oberpfaffenhofen, Germany, and V. Sima, Katholieke Univ. Leuven,
-C Belgium, Nov. 1996.
-C
-C REVISIONS
-C
-C V. Sima, Research Institute for Informatics, Bucharest, July 2003,
-C Feb. 2004.
-C
-C KEYWORDS
-C
-C Multivariable system, state-space model, state-space
-C representation.
-C
-C ******************************************************************
-C
-C .. Parameters ..
- DOUBLE PRECISION ZERO, ONE
- PARAMETER ( ZERO=0.0D0, ONE = 1.0D0 )
-C .. Scalar Arguments ..
- CHARACTER OVER
- INTEGER INFO, LDA, LDA1, LDA2, LDB, LDB1, LDB2, LDC,
- $ LDC1, LDC2, LDD, LDD1, LDD2, M, N, N1, N2, P
- DOUBLE PRECISION ALPHA
-C .. Array Arguments ..
- DOUBLE PRECISION A(LDA,*), A1(LDA1,*), A2(LDA2,*), B(LDB,*),
- $ B1(LDB1,*), B2(LDB2,*), C(LDC,*), C1(LDC1,*),
- $ C2(LDC2,*), D(LDD,*), D1(LDD1,*), D2(LDD2,*)
-C .. Local Scalars ..
- LOGICAL LOVER
- INTEGER I, J, N1P1
-C .. External Functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-C .. External Subroutines ..
- EXTERNAL DAXPY, DLACPY, DLASCL, DLASET, XERBLA
-C .. Intrinsic Functions ..
- INTRINSIC MAX, MIN
-C .. Executable Statements ..
-C
- LOVER = LSAME( OVER, 'O' )
- N = N1 + N2
- INFO = 0
-C
-C Test the input scalar arguments.
-C
- IF( .NOT.LOVER .AND. .NOT.LSAME( OVER, 'N' ) ) THEN
- INFO = -1
- ELSE IF( N1.LT.0 ) THEN
- INFO = -2
- ELSE IF( M.LT.0 ) THEN
- INFO = -3
- ELSE IF( P.LT.0 ) THEN
- INFO = -4
- ELSE IF( N2.LT.0 ) THEN
- INFO = -5
- ELSE IF( LDA1.LT.MAX( 1, N1 ) ) THEN
- INFO = -8
- ELSE IF( LDB1.LT.MAX( 1, N1 ) ) THEN
- INFO = -10
- ELSE IF( ( N1.GT.0 .AND. LDC1.LT.MAX( 1, P ) ) .OR.
- $ ( N1.EQ.0 .AND. LDC1.LT.1 ) ) THEN
- INFO = -12
- ELSE IF( LDD1.LT.MAX( 1, P ) ) THEN
- INFO = -14
- ELSE IF( LDA2.LT.MAX( 1, N2 ) ) THEN
- INFO = -16
- ELSE IF( LDB2.LT.MAX( 1, N2 ) ) THEN
- INFO = -18
- ELSE IF( ( N2.GT.0 .AND. LDC2.LT.MAX( 1, P ) ) .OR.
- $ ( N2.EQ.0 .AND. LDC2.LT.1 ) ) THEN
- INFO = -20
- ELSE IF( LDD2.LT.MAX( 1, P ) ) THEN
- INFO = -22
- ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
- INFO = -25
- ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
- INFO = -27
- ELSE IF( ( N.GT.0 .AND. LDC.LT.MAX( 1, P ) ) .OR.
- $ ( N.EQ.0 .AND. LDC.LT.1 ) ) THEN
- INFO = -29
- ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
- INFO = -31
- END IF
-C
- IF ( INFO.NE.0 ) THEN
-C
-C Error return.
-C
- CALL XERBLA( 'AB05PD', -INFO )
- RETURN
- END IF
-C
-C Quick return if possible.
-C
- IF ( MAX( N, MIN( M, P ) ).EQ.0 )
- $ RETURN
-C
- N1P1 = N1 + 1
-C
-C ( A1 0 )
-C Construct A = ( ) .
-C ( 0 A2 )
-C
- IF ( LOVER .AND. LDA1.LE.LDA ) THEN
- IF ( LDA1.LT.LDA ) THEN
-C
- DO 20 J = N1, 1, -1
- DO 10 I = N1, 1, -1
- A(I,J) = A1(I,J)
- 10 CONTINUE
- 20 CONTINUE
-C
- END IF
- ELSE
- CALL DLACPY( 'F', N1, N1, A1, LDA1, A, LDA )
- END IF
-C
- IF ( N2.GT.0 ) THEN
- CALL DLASET( 'F', N1, N2, ZERO, ZERO, A(1,N1P1), LDA )
- CALL DLASET( 'F', N2, N1, ZERO, ZERO, A(N1P1,1), LDA )
- CALL DLACPY( 'F', N2, N2, A2, LDA2, A(N1P1,N1P1), LDA )
- END IF
-C
-C ( B1 )
-C Construct B = ( ) .
-C ( B2 )
-C
- IF ( LOVER .AND. LDB1.LE.LDB ) THEN
- IF ( LDB1.LT.LDB ) THEN
-C
- DO 40 J = M, 1, -1
- DO 30 I = N1, 1, -1
- B(I,J) = B1(I,J)
- 30 CONTINUE
- 40 CONTINUE
-C
- END IF
- ELSE
- CALL DLACPY( 'F', N1, M, B1, LDB1, B, LDB )
- END IF
-C
- IF ( N2.GT.0 )
- $ CALL DLACPY( 'F', N2, M, B2, LDB2, B(N1P1,1), LDB )
-C
-C Construct C = ( C1 alpha*C2 ) .
-C
- IF ( LOVER .AND. LDC1.LE.LDC ) THEN
- IF ( LDC1.LT.LDC ) THEN
-C
- DO 60 J = N1, 1, -1
- DO 50 I = P, 1, -1
- C(I,J) = C1(I,J)
- 50 CONTINUE
- 60 CONTINUE
-C
- END IF
- ELSE
- CALL DLACPY( 'F', P, N1, C1, LDC1, C, LDC )
- END IF
-C
- IF ( N2.GT.0 ) THEN
- CALL DLACPY( 'F', P, N2, C2, LDC2, C(1,N1P1), LDC )
- IF ( ALPHA.NE.ONE )
- $ CALL DLASCL( 'G', 0, 0, ONE, ALPHA, P, N2, C(1,N1P1), LDC,
- $ INFO )
- END IF
-C
-C Construct D = D1 + alpha*D2 .
-C
- IF ( LOVER .AND. LDD1.LE.LDD ) THEN
- IF ( LDD1.LT.LDD ) THEN
-C
- DO 80 J = M, 1, -1
- DO 70 I = P, 1, -1
- D(I,J) = D1(I,J)
- 70 CONTINUE
- 80 CONTINUE
-C
- END IF
- ELSE
- CALL DLACPY( 'F', P, M, D1, LDD1, D, LDD )
- END IF
-C
- DO 90 J = 1, M
- CALL DAXPY( P, ALPHA, D2(1,J), 1, D(1,J), 1 )
- 90 CONTINUE
-C
- RETURN
-C *** Last line of AB05PD ***
- END
Deleted: trunk/octave-forge/extra/control-devel/src/AB05QD.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB05QD.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB05QD.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,419 +0,0 @@
- SUBROUTINE AB05QD( OVER, N1, M1, P1, N2, M2, P2, A1, LDA1, B1,
- $ LDB1, C1, LDC1, D1, LDD1, A2, LDA2, B2, LDB2,
- $ C2, LDC2, D2, LDD2, N, M, P, A, LDA, B, LDB,
- $ C, LDC, D, LDD, INFO )
-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 append two systems G1 and G2 in state-space form together.
-C If G1 = (A1,B1,C1,D1) and G2 = (A2,B2,C2,D2) are the state-space
-C models of the given two systems having the transfer-function
-C matrices G1 and G2, respectively, this subroutine constructs the
-C state-space model G = (A,B,C,D) which corresponds to the
-C transfer-function matrix
-C
-C ( G1 0 )
-C G = ( )
-C ( 0 G2 )
-C
-C ARGUMENTS
-C
-C Mode Parameters
-C
-C OVER CHARACTER*1
-C Indicates whether the user wishes to overlap pairs of
-C arrays, as follows:
-C = 'N': Do not overlap;
-C = 'O': Overlap pairs of arrays: A1 and A, B1 and B,
-C C1 and C, and D1 and D, i.e. the same name is
-C effectively used for each pair (for all pairs)
-C in the routine call. In this case, setting
-C LDA1 = LDA, LDB1 = LDB, LDC1 = LDC, and LDD1 = LDD
-C will give maximum efficiency.
-C
-C Input/Output Parameters
-C
-C N1 (input) INTEGER
-C The number of state variables in the first system, i.e.
-C the order of the matrix A1, the number of rows of B1 and
-C the number of columns of C1. N1 >= 0.
-C
-C M1 (input) INTEGER
-C The number of input variables in the first system, i.e.
-C the number of columns of matrices B1 and D1. M1 >= 0.
-C
-C P1 (input) INTEGER
-C The number of output variables in the first system, i.e.
-C the number of rows of matrices C1 and D1. P1 >= 0.
-C
-C N2 (input) INTEGER
-C The number of state variables in the second system, i.e.
-C the order of the matrix A2, the number of rows of B2 and
-C the number of columns of C2. N2 >= 0.
-C
-C M2 (input) INTEGER
-C The number of input variables in the second system, i.e.
-C the number of columns of matrices B2 and D2. M2 >= 0.
-C
-C P2 (input) INTEGER
-C The number of output variables in the second system, i.e.
-C the number of rows of matrices C2 and D2. P2 >= 0.
-C
-C A1 (input) DOUBLE PRECISION array, dimension (LDA1,N1)
-C The leading N1-by-N1 part of this array must contain the
-C state transition matrix A1 for the first system.
-C
-C LDA1 INTEGER
-C The leading dimension of array A1. LDA1 >= MAX(1,N1).
-C
-C B1 (input) DOUBLE PRECISION array, dimension (LDB1,M1)
-C The leading N1-by-M1 part of this array must contain the
-C input/state matrix B1 for the first system.
-C
-C LDB1 INTEGER
-C The leading dimension of array B1. LDB1 >= MAX(1,N1).
-C
-C C1 (input) DOUBLE PRECISION array, dimension (LDC1,N1)
-C The leading P1-by-N1 part of this array must contain the
-C state/output matrix C1 for the first system.
-C
-C LDC1 INTEGER
-C The leading dimension of array C1.
-C LDC1 >= MAX(1,P1) if N1 > 0.
-C LDC1 >= 1 if N1 = 0.
-C
-C D1 (input) DOUBLE PRECISION array, dimension (LDD1,M1)
-C The leading P1-by-M1 part of this array must contain the
-C input/output matrix D1 for the first system.
-C
-C LDD1 INTEGER
-C The leading dimension of array D1. LDD1 >= MAX(1,P1).
-C
-C A2 (input) DOUBLE PRECISION array, dimension (LDA2,N2)
-C The leading N2-by-N2 part of this array must contain the
-C state transition matrix A2 for the second system.
-C
-C LDA2 INTEGER
-C The leading dimension of array A2. LDA2 >= MAX(1,N2).
-C
-C B2 (input) DOUBLE PRECISION array, dimension (LDB2,M2)
-C The leading N2-by-M2 part of this array must contain the
-C input/state matrix B2 for the second system.
-C
-C LDB2 INTEGER
-C The leading dimension of array B2. LDB2 >= MAX(1,N2).
-C
-C C2 (input) DOUBLE PRECISION array, dimension (LDC2,N2)
-C The leading P2-by-N2 part of this array must contain the
-C state/output matrix C2 for the second system.
-C
-C LDC2 INTEGER
-C The leading dimension of array C2.
-C LDC2 >= MAX(1,P2) if N2 > 0.
-C LDC2 >= 1 if N2 = 0.
-C
-C D2 (input) DOUBLE PRECISION array, dimension (LDD2,M2)
-C The leading P2-by-M2 part of this array must contain the
-C input/output matrix D2 for the second system.
-C
-C LDD2 INTEGER
-C The leading dimension of array D2. LDD2 >= MAX(1,P2).
-C
-C N (output) INTEGER
-C The number of state variables (N1 + N2) in the resulting
-C system, i.e. the order of the matrix A, the number of rows
-C of B and the number of columns of C.
-C
-C M (output) INTEGER
-C The number of input variables (M1 + M2) in the resulting
-C system, i.e. the number of columns of B and D.
-C
-C P (output) INTEGER
-C The number of output variables (P1 + P2) of the resulting
-C system, i.e. the number of rows of C and D.
-C
-C A (output) DOUBLE PRECISION array, dimension (LDA,N1+N2)
-C The leading N-by-N part of this array contains the state
-C transition matrix A for the resulting system.
-C The array A can overlap A1 if OVER = 'O'.
-C
-C LDA INTEGER
-C The leading dimension of array A. LDA >= MAX(1,N1+N2).
-C
-C B (output) DOUBLE PRECISION array, dimension (LDB,M1+M2)
-C The leading N-by-M part of this array contains the
-C input/state matrix B for the resulting system.
-C The array B can overlap B1 if OVER = 'O'.
-C
-C LDB INTEGER
-C The leading dimension of array B. LDB >= MAX(1,N1+N2).
-C
-C C (output) DOUBLE PRECISION array, dimension (LDC,N1+N2)
-C The leading P-by-N part of this array contains the
-C state/output matrix C for the resulting system.
-C The array C can overlap C1 if OVER = 'O'.
-C
-C LDC INTEGER
-C The leading dimension of array C.
-C LDC >= MAX(1,P1+P2) if N1+N2 > 0.
-C LDC >= 1 if N1+N2 = 0.
-C
-C D (output) DOUBLE PRECISION array, dimension (LDD,M1+M2)
-C The leading P-by-M part of this array contains the
-C input/output matrix D for the resulting system.
-C The array D can overlap D1 if OVER = 'O'.
-C
-C LDD INTEGER
-C The leading dimension of array D. LDD >= MAX(1,P1+P2).
-C
-C Error Indicator
-C
-C INFO INTEGER
-C = 0: successful exit;
-C < 0: if INFO = -i, the i-th argument had an illegal
-C value.
-C
-C METHOD
-C
-C The matrices of the resulting systems are determined as:
-C
-C ( A1 0 ) ( B1 0 )
-C A = ( ) , B = ( ) ,
-C ( 0 A2 ) ( 0 B2 )
-C
-C ( C1 0 ) ( D1 0 )
-C C = ( ) , D = ( ) .
-C ( 0 C2 ) ( 0 D2 )
-C
-C REFERENCES
-C
-C None
-C
-C CONTRIBUTORS
-C
-C A. Varga, German Aerospace Research Establishment,
-C Oberpfaffenhofen, Germany, and V. Sima, Katholieke Univ. Leuven,
-C Belgium, Nov. 1996.
-C
-C REVISIONS
-C
-C V. Sima, Research Institute for Informatics, Bucharest, Feb. 2004.
-C
-C KEYWORDS
-C
-C Multivariable system, state-space model, state-space
-C representation.
-C
-C ******************************************************************
-C
-C .. Parameters ..
- DOUBLE PRECISION ZERO
- PARAMETER ( ZERO=0.0D0 )
-C .. Scalar Arguments ..
- CHARACTER OVER
- INTEGER INFO, LDA, LDA1, LDA2, LDB, LDB1, LDB2, LDC,
- $ LDC1, LDC2, LDD, LDD1, LDD2, M, M1, M2, N, N1,
- $ N2, P, P1, P2
-C .. Array Arguments ..
- DOUBLE PRECISION A(LDA,*), A1(LDA1,*), A2(LDA2,*), B(LDB,*),
- $ B1(LDB1,*), B2(LDB2,*), C(LDC,*), C1(LDC1,*),
- $ C2(LDC2,*), D(LDD,*), D1(LDD1,*), D2(LDD2,*)
-C .. Local Scalars ..
- LOGICAL LOVER
- INTEGER I, J
-C .. External Functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-C .. External Subroutines ..
- EXTERNAL DLACPY, DLASET, XERBLA
-C .. Intrinsic Functions ..
- INTRINSIC MAX, MIN
-C .. Executable Statements ..
-C
- LOVER = LSAME( OVER, 'O' )
- N = N1 + N2
- M = M1 + M2
- P = P1 + P2
- INFO = 0
-C
-C Test the input scalar arguments.
-C
- IF( .NOT.LOVER .AND. .NOT.LSAME( OVER, 'N' ) ) THEN
- INFO = -1
- ELSE IF( N1.LT.0 ) THEN
- INFO = -2
- ELSE IF( M1.LT.0 ) THEN
- INFO = -3
- ELSE IF( P1.LT.0 ) THEN
- INFO = -4
- ELSE IF( N2.LT.0 ) THEN
- INFO = -5
- ELSE IF( M2.LT.0 ) THEN
- INFO = -6
- ELSE IF( P2.LT.0 ) THEN
- INFO = -7
- ELSE IF( LDA1.LT.MAX( 1, N1 ) ) THEN
- INFO = -9
- ELSE IF( LDB1.LT.MAX( 1, N1 ) ) THEN
- INFO = -11
- ELSE IF( ( N1.GT.0 .AND. LDC1.LT.MAX( 1, P1 ) ) .OR.
- $ ( N1.EQ.0 .AND. LDC1.LT.1 ) ) THEN
- INFO = -13
- ELSE IF( LDD1.LT.MAX( 1, P1 ) ) THEN
- INFO = -15
- ELSE IF( LDA2.LT.MAX( 1, N2 ) ) THEN
- INFO = -17
- ELSE IF( LDB2.LT.MAX( 1, N2 ) ) THEN
- INFO = -19
- ELSE IF( ( N2.GT.0 .AND. LDC2.LT.MAX( 1, P2 ) ) .OR.
- $ ( N2.EQ.0 .AND. LDC2.LT.1 ) ) THEN
- INFO = -21
- ELSE IF( LDD2.LT.MAX( 1, P2 ) ) THEN
- INFO = -23
- ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
- INFO = -28
- ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
- INFO = -30
- ELSE IF( ( N.GT.0 .AND. LDC.LT.MAX( 1, P ) ) .OR.
- $ ( N.EQ.0 .AND. LDC.LT.1 ) ) THEN
- INFO = -32
- ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
- INFO = -34
- END IF
-C
- IF ( INFO.NE.0 ) THEN
-C
-C Error return.
-C
- CALL XERBLA( 'AB05QD', -INFO )
- RETURN
- END IF
-C
-C Quick return if possible.
-C
- IF ( MAX( N, MIN( M, P ) ).EQ.0 )
- $ RETURN
-C ( A1 0 )
-C Construct A = ( ) .
-C ( 0 A2 )
-C
- IF ( LOVER .AND. LDA1.LE.LDA ) THEN
- IF ( LDA1.LT.LDA ) THEN
-C
- DO 20 J = N1, 1, -1
- DO 10 I = N1, 1, -1
- A(I,J) = A1(I,J)
- 10 CONTINUE
- 20 CONTINUE
-C
- END IF
- ELSE
- CALL DLACPY( 'F', N1, N1, A1, LDA1, A, LDA )
- END IF
-C
- IF ( N2.GT.0 ) THEN
- CALL DLASET( 'F', N1, N2, ZERO, ZERO, A(1,N1+1), LDA )
- CALL DLASET( 'F', N2, N1, ZERO, ZERO, A(N1+1,1), LDA )
- CALL DLACPY( 'F', N2, N2, A2, LDA2, A(N1+1,N1+1), LDA )
- END IF
-C
-C ( B1 0 )
-C Construct B = ( ) .
-C ( 0 B2 )
-C
- IF ( LOVER .AND. LDB1.LE.LDB ) THEN
- IF ( LDB1.LT.LDB ) THEN
-C
- DO 40 J = M1, 1, -1
- DO 30 I = N1, 1, -1
- B(I,J) = B1(I,J)
- 30 CONTINUE
- 40 CONTINUE
-C
- END IF
- ELSE
- CALL DLACPY( 'F', N1, M1, B1, LDB1, B, LDB )
- END IF
-C
- IF ( M2.GT.0 )
- $ CALL DLASET( 'F', N1, M2, ZERO, ZERO, B(1,M1+1), LDB )
- IF ( N2.GT.0 ) THEN
- CALL DLASET( 'F', N2, M1, ZERO, ZERO, B(N1+1,1), LDB )
- IF ( M2.GT.0 )
- $ CALL DLACPY( 'F', N2, M2, B2, LDB2, B(N1+1,M1+1), LDB )
- END IF
-C
-C ( C1 0 )
-C Construct C = ( ) .
-C ( 0 C2 )
-C
- IF ( LOVER .AND. LDC1.LE.LDC ) THEN
- IF ( LDC1.LT.LDC ) THEN
-C
- DO 60 J = N1, 1, -1
- DO 50 I = P1, 1, -1
- C(I,J) = C1(I,J)
- 50 CONTINUE
- 60 CONTINUE
-C
- END IF
- ELSE
- CALL DLACPY( 'F', P1, N1, C1, LDC1, C, LDC )
- END IF
-C
- IF ( N2.GT.0 )
- $ CALL DLASET( 'F', P1, N2, ZERO, ZERO, C(1,N1+1), LDC )
- IF ( P2.GT.0 ) THEN
- IF ( N1.GT.0 )
- $ CALL DLASET( 'F', P2, N1, ZERO, ZERO, C(P1+1,1), LDC )
- IF ( N2.GT.0 )
- $ CALL DLACPY( 'F', P2, N2, C2, LDC2, C(P1+1,N1+1), LDC )
- END IF
-C
-C ( D1 0 )
-C Construct D = ( ) .
-C ( 0 D2 )
-C
- IF ( LOVER .AND. LDD1.LE.LDD ) THEN
- IF ( LDD1.LT.LDD ) THEN
-C
- DO 80 J = M1, 1, -1
- DO 70 I = P1, 1, -1
- D(I,J) = D1(I,J)
- 70 CONTINUE
- 80 CONTINUE
-C
- END IF
- ELSE
- CALL DLACPY( 'F', P1, M1, D1, LDD1, D, LDD )
- END IF
-C
- IF ( M2.GT.0 )
- $ CALL DLASET( 'F', P1, M2, ZERO, ZERO, D(1,M1+1), LDD )
- IF ( P2.GT.0 ) THEN
- CALL DLASET( 'F', P2, M1, ZERO, ZERO, D(P1+1,1), LDD )
- IF ( M2.GT.0 )
- $ CALL DLACPY( 'F', P2, M2, D2, LDD2, D(P1+1,M1+1), LDD )
- END IF
-C
- RETURN
-C *** Last line of AB05QD ***
- END
Deleted: trunk/octave-forge/extra/control-devel/src/AB07MD.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB07MD.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB07MD.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,224 +0,0 @@
- SUBROUTINE AB07MD( JOBD, N, M, P, A, LDA, B, LDB, C, LDC, D, LDD,
- $ INFO )
-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 find the dual of a given state-space representation.
-C
-C ARGUMENTS
-C
-C Mode Parameters
-C
-C JOBD CHARACTER*1
-C Specifies whether or not a non-zero matrix D appears in
-C the given state space model:
-C = 'D': D is present;
-C = 'Z': D is assumed a zero matrix.
-C
-C Input/Output Parameters
-C
-C N (input) INTEGER
-C The order of the state-space representation. N >= 0.
-C
-C M (input) INTEGER
-C The number of system inputs. M >= 0.
-C
-C P (input) INTEGER
-C The number of system outputs. P >= 0.
-C
-C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C On entry, the leading N-by-N part of this array must
-C contain the original state dynamics matrix A.
-C On exit, the leading N-by-N part of this array contains
-C the dual state dynamics matrix A'.
-C
-C LDA INTEGER
-C The leading dimension of array A. LDA >= MAX(1,N).
-C
-C B (input/output) DOUBLE PRECISION array, dimension
-C (LDB,MAX(M,P))
-C On entry, the leading N-by-M part of this array must
-C contain the original input/state matrix B.
-C On exit, the leading N-by-P part of this array contains
-C the dual input/state matrix C'.
-C
-C LDB INTEGER
-C The leading dimension of array B. LDB >= MAX(1,N).
-C
-C C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C On entry, the leading P-by-N part of this array must
-C contain the original state/output matrix C.
-C On exit, the leading M-by-N part of this array contains
-C the dual state/output matrix B'.
-C
-C LDC INTEGER
-C The leading dimension of array C.
-C LDC >= MAX(1,M,P) if N > 0.
-C LDC >= 1 if N = 0.
-C
-C D (input/output) DOUBLE PRECISION array, dimension
-C (LDD,MAX(M,P))
-C On entry, if JOBD = 'D', the leading P-by-M part of this
-C array must contain the original direct transmission
-C matrix D.
-C On exit, if JOBD = 'D', the leading M-by-P part of this
-C array contains the dual direct transmission matrix D'.
-C The array D is not referenced if JOBD = 'Z'.
-C
-C LDD INTEGER
-C The leading dimension of array D.
-C LDD >= MAX(1,M,P) if JOBD = 'D'.
-C LDD >= 1 if JOBD = 'Z'.
-C
-C Error Indicator
-C
-C INFO INTEGER
-C = 0: successful exit;
-C < 0: if INFO = -i, the i-th argument had an illegal
-C value.
-C
-C METHOD
-C
-C If the given state-space representation is the M-input/P-output
-C (A,B,C,D), its dual is simply the P-input/M-output (A',C',B',D').
-C
-C REFERENCES
-C
-C None
-C
-C NUMERICAL ASPECTS
-C
-C None
-C
-C CONTRIBUTOR
-C
-C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Dec. 1996.
-C Supersedes Release 2.0 routine AB07AD by T.W.C.Williams, Kingston
-C Polytechnic, United Kingdom, March 1982.
-C
-C REVISIONS
-C
-C V. Sima, Research Institute for Informatics, Bucharest, Feb. 2004.
-C
-C KEYWORDS
-C
-C Dual system, state-space model, state-space representation.
-C
-C ******************************************************************
-C
-C .. Scalar Arguments ..
- CHARACTER JOBD
- INTEGER INFO, LDA, LDB, LDC, LDD, M, N, P
-C .. Array Arguments ..
- DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*)
-C .. Local Scalars ..
- LOGICAL LJOBD
- INTEGER J, MINMP, MPLIM
-C .. External functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-C .. External subroutines ..
- EXTERNAL DCOPY, DSWAP, XERBLA
-C .. Intrinsic Functions ..
- INTRINSIC MAX, MIN
-C .. Executable Statements ..
-C
- INFO = 0
- LJOBD = LSAME( JOBD, 'D' )
- MPLIM = MAX( M, P )
- MINMP = MIN( M, P )
-C
-C Test the input scalar arguments.
-C
- IF( .NOT.LJOBD .AND. .NOT.LSAME( JOBD, 'Z' ) ) THEN
- INFO = -1
- ELSE IF( N.LT.0 ) THEN
- INFO = -2
- ELSE IF( M.LT.0 ) THEN
- INFO = -3
- ELSE IF( P.LT.0 ) THEN
- INFO = -4
- ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
- INFO = -6
- ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
- INFO = -8
- ELSE IF( ( N.GT.0 .AND. LDC.LT.MAX( 1, MPLIM ) ) .OR.
- $ ( N.EQ.0 .AND. LDC.LT.1 ) ) THEN
- INFO = -10
- ELSE IF( ( LJOBD .AND. LDD.LT.MAX( 1, MPLIM ) ) .OR.
- $ ( .NOT.LJOBD .AND. LDD.LT.1 ) ) THEN
- INFO = -12
- END IF
-C
- IF ( INFO.NE.0 ) THEN
-C
-C Error return.
-C
- CALL XERBLA( 'AB07MD', -INFO )
- RETURN
- END IF
-C
-C Quick return if possible.
-C
- IF ( MAX( N, MINMP ).EQ.0 )
- $ RETURN
-C
- IF ( N.GT.0 ) THEN
-C
-C Transpose A, if non-scalar.
-C
- DO 10 J = 1, N - 1
- CALL DSWAP( N-J, A(J+1,J), 1, A(J,J+1), LDA )
- 10 CONTINUE
-C
-C Replace B by C' and C by B'.
-C
- DO 20 J = 1, MPLIM
- IF ( J.LE.MINMP ) THEN
- CALL DSWAP( N, B(1,J), 1, C(J,1), LDC )
- ELSE IF ( J.GT.P ) THEN
- CALL DCOPY( N, B(1,J), 1, C(J,1), LDC )
- ELSE
- CALL DCOPY( N, C(J,1), LDC, B(1,J), 1 )
- END IF
- 20 CONTINUE
-C
- END IF
-C
- IF ( LJOBD .AND. MINMP.GT.0 ) THEN
-C
-C Transpose D, if non-scalar.
-C
- DO 30 J = 1, MPLIM
- IF ( J.LT.MINMP ) THEN
- CALL DSWAP( MINMP-J, D(J+1,J), 1, D(J,J+1), LDD )
- ELSE IF ( J.GT.P ) THEN
- CALL DCOPY( P, D(1,J), 1, D(J,1), LDD )
- ELSE IF ( J.GT.M ) THEN
- CALL DCOPY( M, D(J,1), LDD, D(1,J), 1 )
- END IF
- 30 CONTINUE
-C
- END IF
-C
- RETURN
-C *** Last line of AB07MD ***
- END
Deleted: trunk/octave-forge/extra/control-devel/src/AB07ND.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB07ND.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB07ND.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,303 +0,0 @@
- SUBROUTINE AB07ND( N, M, A, LDA, B, LDB, C, LDC, D, LDD, RCOND,
- $ IWORK, DWORK, LDWORK, INFO )
-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 inverse (Ai,Bi,Ci,Di) of a given system (A,B,C,D).
-C
-C ARGUMENTS
-C
-C Input/Output Parameters
-C
-C N (input) INTEGER
-C The order of the state matrix A. N >= 0.
-C
-C M (input) INTEGER
-C The number of system inputs and outputs. M >= 0.
-C
-C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C On entry, the leading N-by-N part of this array must
-C contain the state matrix A of the original system.
-C On exit, the leading N-by-N part of this array contains
-C the state matrix Ai of the inverse system.
-C
-C LDA INTEGER
-C The leading dimension of the array A. LDA >= MAX(1,N).
-C
-C B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
-C On entry, the leading N-by-M part of this array must
-C contain the input matrix B of the original system.
-C On exit, the leading N-by-M part of this array contains
-C the input matrix Bi of the inverse system.
-C
-C LDB INTEGER
-C The leading dimension of the array B. LDB >= MAX(1,N).
-C
-C C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C On entry, the leading M-by-N part of this array must
-C contain the output matrix C of the original system.
-C On exit, the leading M-by-N part of this array contains
-C the output matrix Ci of the inverse system.
-C
-C LDC INTEGER
-C The leading dimension of the array C. LDC >= MAX(1,M).
-C
-C D (input/output) DOUBLE PRECISION array, dimension (LDD,M)
-C On entry, the leading M-by-M part of this array must
-C contain the feedthrough matrix D of the original system.
-C On exit, the leading M-by-M part of this array contains
-C the feedthrough matrix Di of the inverse system.
-C
-C LDD INTEGER
-C The leading dimension of the array D. LDD >= MAX(1,M).
-C
-C RCOND (output) DOUBLE PRECISION
-C The estimated reciprocal condition number of the
-C feedthrough matrix D of the original system.
-C
-C Workspace
-C
-C IWORK INTEGER array, dimension (2*M)
-C
-C DWORK DOUBLE PRECISION array, dimension (LDWORK)
-C On exit, if INFO = 0 or M+1, DWORK(1) returns the optimal
-C value of LDWORK.
-C
-C LDWORK INTEGER
-C The length of the array DWORK. LDWORK >= MAX(1,4*M).
-C For good performance, LDWORK should be larger.
-C
-C Error Indicator
-C
-C INFO INTEGER
-C = 0: successful exit;
-C < 0: if INFO = -i, the i-th argument had an illegal
-C value;
-C = i: the matrix D is exactly singular; the (i,i) diagonal
-C element is zero, i <= M; RCOND was set to zero;
-C = M+1: the matrix D is numerically singular, i.e., RCOND
-C is less than the relative machine precision, EPS
-C (see LAPACK Library routine DLAMCH). The
-C calculations have been completed, but the results
-C could be very inaccurate.
-C
-C METHOD
-C
-C The matrices of the inverse system are computed with the formulas:
-C -1 -1 -1 -1
-C Ai = A - B*D *C, Bi = -B*D , Ci = D *C, Di = D .
-C
-C NUMERICAL ASPECTS
-C
-C The accuracy depends mainly on the condition number of the matrix
-C D to be inverted. The estimated reciprocal condition number is
-C returned in RCOND.
-C
-C CONTRIBUTORS
-C
-C A. Varga, German Aerospace Center, Oberpfaffenhofen, March 2000.
-C D. Sima, University of Bucharest, April 2000.
-C V. Sima, Research Institute for Informatics, Bucharest, Apr. 2000.
-C Based on the routine SYSINV, A. Varga, 1992.
-C
-C REVISIONS
-C
-C A. Varga, German Aerospace Center, Oberpfaffenhofen, July 2000.
-C
-C KEYWORDS
-C
-C Inverse system, state-space model, state-space representation.
-C
-C ******************************************************************
-C
-C .. Parameters ..
- DOUBLE PRECISION ZERO, ONE
- PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
-C .. Scalar Arguments ..
- DOUBLE PRECISION RCOND
- INTEGER INFO, LDA, LDB, LDC, LDD, LDWORK, M, N
-C .. Array Arguments ..
- DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*),
- $ DWORK(*)
- INTEGER IWORK(*)
-C .. Local Scalars ..
- DOUBLE PRECISION DNORM
- INTEGER BL, CHUNK, I, IERR, J, MAXWRK
- LOGICAL BLAS3, BLOCK
-C .. External Functions ..
- DOUBLE PRECISION DLAMCH, DLANGE
- INTEGER ILAENV
- EXTERNAL DLAMCH, DLANGE, ILAENV
-C .. External Subroutines ..
- EXTERNAL DCOPY, DGECON, DGEMM, DGEMV, DGETRF, DGETRI,
- $ DLACPY, XERBLA
-C .. Intrinsic Functions ..
- INTRINSIC DBLE, MAX, MIN
-C .. Executable Statements ..
-C
- INFO = 0
-C
-C Test the input scalar arguments.
-C
- IF( N.LT.0 ) THEN
- INFO = -1
- ELSE IF( M.LT.0 ) THEN
- INFO = -2
- ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
- INFO = -4
- ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
- INFO = -6
- ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
- INFO = -8
- ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
- INFO = -10
- ELSE IF( LDWORK.LT.MAX( 1, 4*M ) ) THEN
- INFO = -14
- END IF
-C
- IF ( INFO.NE.0 ) THEN
-C
-C Error return.
-C
- CALL XERBLA( 'AB07ND', -INFO )
- RETURN
- END IF
-C
-C Quick return if possible.
-C
- IF ( M.EQ.0 ) THEN
- RCOND = ONE
- DWORK(1) = ONE
- RETURN
- END IF
-C
-C Factorize D.
-C
- CALL DGETRF( M, M, D, LDD, IWORK, INFO )
- IF ( INFO.NE.0 ) THEN
- RCOND = ZERO
- RETURN
- END IF
-C
-C Compute the reciprocal condition number of the matrix D.
-C Workspace: need 4*M.
-C (Note: Comments in the code beginning "Workspace:" describe the
-C minimal amount of workspace needed at that point in the code,
-C as well as the preferred amount for good performance.
-C NB refers to the optimal block size for the immediately
-C following subroutine, as returned by ILAENV.)
-C
- DNORM = DLANGE( '1-norm', M, M, D, LDD, DWORK )
- CALL DGECON( '1-norm', M, D, LDD, DNORM, RCOND, DWORK, IWORK(M+1),
- $ IERR )
- IF ( RCOND.LT.DLAMCH( 'Epsilon' ) )
- $ INFO = M + 1
-C -1
-C Compute Di = D .
-C Workspace: need M;
-C prefer M*NB.
-C
- MAXWRK = MAX( 4*M, M*ILAENV( 1, 'DGETRI', ' ', M, -1, -1, -1 ) )
- CALL DGETRI( M, D, LDD, IWORK, DWORK, LDWORK, IERR )
- IF ( N.GT.0 ) THEN
- CHUNK = LDWORK / M
- BLAS3 = CHUNK.GE.N .AND. M.GT.1
- BLOCK = MIN( CHUNK, M ).GT.1
-C -1
-C Compute Bi = -B*D .
-C
- IF ( BLAS3 ) THEN
-C
-C Enough workspace for a fast BLAS 3 algorithm.
-C
- CALL DLACPY( 'Full', N, M, B, LDB, DWORK, N )
- CALL DGEMM( 'NoTranspose', 'NoTranspose', N, M, M, -ONE,
- $ DWORK, N, D, LDD, ZERO, B, LDB )
-C
- ELSE IF( BLOCK ) THEN
-C
-C Use as many rows of B as possible.
-C
- DO 10 I = 1, N, CHUNK
- BL = MIN( N-I+1, CHUNK )
- CALL DLACPY( 'Full', BL, M, B(I,1), LDB, DWORK, BL )
- CALL DGEMM( 'NoTranspose', 'NoTranspose', BL, M, M, -ONE,
- $ DWORK, BL, D, LDD, ZERO, B(I,1), LDB )
- 10 CONTINUE
-C
- ELSE
-C
-C Use a BLAS 2 algorithm.
-C
- DO 20 I = 1, N
- CALL DCOPY( M, B(I,1), LDB, DWORK, 1 )
- CALL DGEMV( 'Transpose', M, M, -ONE, D, LDD, DWORK, 1,
- $ ZERO, B(I,1), LDB )
- 20 CONTINUE
-C
- END IF
-C
-C Compute Ai = A + Bi*C.
-C
- CALL DGEMM( 'NoTranspose', 'NoTranspose', N, N, M, ONE, B, LDB,
- $ C, LDC, ONE, A, LDA )
-C -1
-C Compute C <-- D *C.
-C
- IF ( BLAS3 ) THEN
-C
-C Enough workspace for a fast BLAS 3 algorithm.
-C
- CALL DLACPY( 'Full', M, N, C, LDC, DWORK, M )
- CALL DGEMM( 'NoTranspose', 'NoTranspose', M, N, M, ONE,
- $ D, LDD, DWORK, M, ZERO, C, LDC )
-C
- ELSE IF( BLOCK ) THEN
-C
-C Use as many columns of C as possible.
-C
- DO 30 J = 1, N, CHUNK
- BL = MIN( N-J+1, CHUNK )
- CALL DLACPY( 'Full', M, BL, C(1,J), LDC, DWORK, M )
- CALL DGEMM( 'NoTranspose', 'NoTranspose', M, BL, M, ONE,
- $ D, LDD, DWORK, M, ZERO, C(1,J), LDC )
- 30 CONTINUE
-C
- ELSE
-C
-C Use a BLAS 2 algorithm.
-C
- DO 40 J = 1, N
- CALL DCOPY( M, C(1,J), 1, DWORK, 1 )
- CALL DGEMV( 'NoTranspose', M, M, ONE, D, LDD, DWORK, 1,
- $ ZERO, C(1,J), 1 )
- 40 CONTINUE
-C
- END IF
- END IF
-C
-C Return optimal workspace in DWORK(1).
-C
- DWORK(1) = DBLE( MAX( MAXWRK, N*M ) )
- RETURN
-C
-C *** Last line of AB07ND ***
- END
Deleted: trunk/octave-forge/extra/control-devel/src/AB08MD.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB08MD.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB08MD.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,299 +0,0 @@
- SUBROUTINE AB08MD( EQUIL, N, M, P, A, LDA, B, LDB, C, LDC, D, LDD,
- $ RANK, TOL, IWORK, DWORK, LDWORK, INFO )
-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 normal rank of the transfer-function matrix of a
-C state-space model (A,B,C,D).
-C
-C ARGUMENTS
-C
-C Mode Parameters
-C
-C EQUIL CHARACTER*1
-C Specifies whether the user wishes to balance the compound
-C matrix (see METHOD) as follows:
-C = 'S': Perform balancing (scaling);
-C = 'N': Do not perform balancing.
-C
-C Input/Output Parameters
-C
-C N (input) INTEGER
-C The number of state variables, i.e., the order of the
-C matrix A. N >= 0.
-C
-C M (input) INTEGER
-C The number of system inputs. M >= 0.
-C
-C P (input) INTEGER
-C The number of system outputs. P >= 0.
-C
-C A (input) DOUBLE PRECISION array, dimension (LDA,N)
-C The leading N-by-N part of this array must contain the
-C state dynamics matrix A of the system.
-C
-C LDA INTEGER
-C The leading dimension of array A. LDA >= MAX(1,N).
-C
-C B (input) DOUBLE PRECISION array, dimension (LDB,M)
-C The leading N-by-M part of this array must contain the
-C input/state matrix B of the system.
-C
-C LDB INTEGER
-C The leading dimension of array B. LDB >= MAX(1,N).
-C
-C C (input) DOUBLE PRECISION array, dimension (LDC,N)
-C The leading P-by-N part of this array must contain the
-C state/output matrix C of the system.
-C
-C LDC INTEGER
-C The leading dimension of array C. LDC >= MAX(1,P).
-C
-C D (input) DOUBLE PRECISION array, dimension (LDD,M)
-C The leading P-by-M part of this array must contain the
-C direct transmission matrix D of the system.
-C
-C LDD INTEGER
-C The leading dimension of array D. LDD >= MAX(1,P).
-C
-C RANK (output) INTEGER
-C The normal rank of the transfer-function matrix.
-C
-C Tolerances
-C
-C TOL DOUBLE PRECISION
-C A tolerance used in rank decisions to determine the
-C effective rank, which is defined as the order of the
-C largest leading (or trailing) triangular submatrix in the
-C QR (or RQ) factorization with column (or row) pivoting
-C whose estimated condition number is less than 1/TOL.
-C If the user sets TOL to be less than SQRT((N+P)*(N+M))*EPS
-C then the tolerance is taken as SQRT((N+P)*(N+M))*EPS,
-C where EPS is the machine precision (see LAPACK Library
-C Routine DLAMCH).
-C
-C Workspace
-C
-C IWORK INTEGER array, dimension (2*N+MAX(M,P)+1)
-C
-C DWORK DOUBLE PRECISION array, dimension (LDWORK)
-C On exit, if INFO = 0, DWORK(1) returns the optimal value
-C of LDWORK.
-C
-C LDWORK INTEGER
-C The length of the array DWORK.
-C LDWORK >= (N+P)*(N+M) +
-C MAX( MIN(P,M) + MAX(3*M-1,N), 1,
-C MIN(P,N) + MAX(3*P-1,N+P,N+M) )
-C For optimum performance LDWORK should be larger.
-C
-C If LDWORK = -1, then a workspace query is assumed;
-C the routine only calculates the optimal size of the
-C DWORK array, returns this value as the first entry of
-C the DWORK array, and no error message related to LDWORK
-C is issued by XERBLA.
-C
-C Error Indicator
-C
-C INFO INTEGER
-C = 0: successful exit;
-C < 0: if INFO = -i, the i-th argument had an illegal
-C value.
-C
-C METHOD
-C
-C The routine reduces the (N+P)-by-(M+N) compound matrix (B A)
-C (D C)
-C
-C to one with the same invariant zeros and with D of full row rank.
-C The normal rank of the transfer-function matrix is the rank of D.
-C
-C REFERENCES
-C
-C [1] Svaricek, F.
-C Computation of the Structural Invariants of Linear
-C Multivariable Systems with an Extended Version of
-C the Program ZEROS.
-C System & Control Letters, 6, pp. 261-266, 1985.
-C
-C [2] Emami-Naeini, A. and Van Dooren, P.
-C Computation of Zeros of Linear Multivariable Systems.
-C Automatica, 18, pp. 415-430, 1982.
-C
-C NUMERICAL ASPECTS
-C
-C The algorithm is backward stable (see [2] and [1]).
-C
-C CONTRIBUTOR
-C
-C A. Varga, German Aerospace Center, Oberpfaffenhofen, May 2001.
-C
-C REVISIONS
-C
-C V. Sima, Research Institute for Informatics, Bucharest, June 2001,
-C Dec. 2003, Jan. 2009, Mar. 2009, Apr. 2009.
-C
-C KEYWORDS
-C
-C Multivariable system, orthogonal transformation,
-C structural invariant.
-C
-C ******************************************************************
-C
-C .. Parameters ..
- DOUBLE PRECISION ZERO, ONE
- PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
-C .. Scalar Arguments ..
- CHARACTER EQUIL
- INTEGER INFO, LDA, LDB, LDC, LDD, LDWORK, M, N, P, RANK
- DOUBLE PRECISION TOL
-C .. Array Arguments ..
- INTEGER IWORK(*)
- DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(LDD,*), DWORK(*)
-C .. Local Scalars ..
- LOGICAL LEQUIL, LQUERY
- INTEGER I, KW, MU, NB, NINFZ, NKROL, NM, NP, NU, RO,
- $ SIGMA, WRKOPT
- DOUBLE PRECISION MAXRED, SVLMAX, THRESH, TOLER
-C .. External Functions ..
- LOGICAL LSAME
- DOUBLE PRECISION DLAMCH, DLANGE
- EXTERNAL DLAMCH, DLANGE, LSAME
-C .. External Subroutines ..
- EXTERNAL AB08NX, DLACPY, TB01ID, XERBLA
-C .. Intrinsic Functions ..
- INTRINSIC DBLE, INT, MAX, MIN, SQRT
-C .. Executable Statements ..
-C
- NP = N + P
- NM = N + M
- INFO = 0
- LEQUIL = LSAME( EQUIL, 'S' )
- LQUERY = ( LDWORK.EQ.-1 )
- WRKOPT = NP*NM
-C
-C Test the input scalar arguments.
-C
- IF( .NOT.LEQUIL .AND. .NOT.LSAME( EQUIL, 'N' ) ) THEN
- INFO = -1
- ELSE IF( N.LT.0 ) THEN
- INFO = -2
- ELSE IF( M.LT.0 ) THEN
- INFO = -3
- ELSE IF( P.LT.0 ) THEN
- INFO = -4
- ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
- INFO = -6
- ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
- INFO = -8
- ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
- INFO = -10
- ELSE IF( LDD.LT.MAX( 1, P ) ) THEN
- INFO = -12
- ELSE
- KW = WRKOPT + MAX( MIN( P, M ) + MAX( 3*M-1, N ), 1,
- $ MIN( P, N ) + MAX( 3*P-1, NP, NM ) )
- IF( LQUERY ) THEN
- SVLMAX = ZERO
- NINFZ = 0
- CALL AB08NX( N, M, P, P, 0, SVLMAX, DWORK, MAX( 1, NP ),
- $ NINFZ, IWORK, IWORK, MU, NU, NKROL, TOL, IWORK,
- $ DWORK, -1, INFO )
- WRKOPT = MAX( KW, WRKOPT + INT( DWORK(1) ) )
- ELSE IF( LDWORK.LT.KW ) THEN
- INFO = -17
- END IF
- END IF
-C
- IF ( INFO.NE.0 ) THEN
-C
-C Error return.
-C
- CALL XERBLA( 'AB08MD', -INFO )
- RETURN
- ELSE IF( LQUERY ) THEN
- DWORK(1) = WRKOPT
- RETURN
- END IF
-C
-C Quick return if possible.
-C
- IF ( MIN( M, P ).EQ.0 ) THEN
- RANK = 0
- DWORK(1) = ONE
- RETURN
- END IF
-C
- DO 10 I = 1, 2*N+1
- IWORK(I) = 0
- 10 CONTINUE
-C
-C (Note: Comments in the code beginning "Workspace:" describe the
-C minimal amount of real workspace needed at that point in the
-C code, as well as the preferred amount for good performance.)
-C
-C Construct the compound matrix ( B A ), dimension (N+P)-by-(M+N).
-C ( D C )
-C Workspace: need (N+P)*(N+M).
-C
- CALL DLACPY( 'Full', N, M, B, LDB, DWORK, NP )
- CALL DLACPY( 'Full', P, M, D, LDD, DWORK(N+1), NP )
- CALL DLACPY( 'Full', N, N, A, LDA, DWORK(NP*M+1), NP )
- CALL DLACPY( 'Full', P, N, C, LDC, DWORK(NP*M+N+1), NP )
-C
-C If required, balance the compound matrix (default MAXRED).
-C Workspace: need N.
-C
- KW = WRKOPT + 1
- IF ( LEQUIL ) THEN
- MAXRED = ZERO
- CALL TB01ID( 'A', N, M, P, MAXRED, DWORK(NP*M+1), NP, DWORK,
- $ NP, DWORK(NP*M+N+1), NP, DWORK(KW), INFO )
- WRKOPT = WRKOPT + N
- END IF
-C
-C If required, set tolerance.
-C
- THRESH = SQRT( DBLE( NP*NM ) )*DLAMCH( 'Precision' )
- TOLER = TOL
- IF ( TOLER.LT.THRESH ) TOLER = THRESH
- SVLMAX = DLANGE( 'Frobenius', NP, NM, DWORK, NP, DWORK(KW) )
-C
-C Reduce this system to one with the same invariant zeros and with
-C D full row rank MU (the normal rank of the original system).
-C Real workspace: need (N+P)*(N+M) +
-C MAX( 1, MIN(P,M) + MAX(3*M-1,N),
-C MIN(P,N) + MAX(3*P-1,N+P,N+M) );
-C prefer larger.
-C Integer workspace: 2*N+MAX(M,P)+1.
-C
- RO = P
- SIGMA = 0
- NINFZ = 0
- CALL AB08NX( N, M, P, RO, SIGMA, SVLMAX, DWORK, NP, NINFZ, IWORK,
- $ IWORK(N+1), MU, NU, NKROL, TOLER, IWORK(2*N+2),
- $ DWORK(KW), LDWORK-KW+1, INFO )
- RANK = MU
-C
- DWORK(1) = MAX( WRKOPT, INT( DWORK(KW) ) + KW - 1 )
- RETURN
-C *** Last line of AB08MD ***
- END
Deleted: trunk/octave-forge/extra/control-devel/src/AB08NX.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB08NX.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB08NX.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,446 +0,0 @@
- SUBROUTINE AB08NX( N, M, P, RO, SIGMA, SVLMAX, ABCD, LDABCD,
- $ NINFZ, INFZ, KRONL, MU, NU, NKROL, TOL, IWORK,
- $ DWORK, LDWORK, INFO )
-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 extract from the (N+P)-by-(M+N) system
-C ( B A )
-C ( D C )
-C an (NU+MU)-by-(M+NU) "reduced" system
-C ( B' A')
-C ( D' C')
-C having the same transmission zeros but with D' of full row rank.
-C
-C ARGUMENTS
-C
-C Input/Output Parameters
-C
-C N (input) INTEGER
-C The number of state variables. N >= 0.
-C
-C M (input) INTEGER
-C The number of system inputs. M >= 0.
-C
-C P (input) INTEGER
-C The number of system outputs. P >= 0.
-C
-C RO (input/output) INTEGER
-C On entry,
-C = P for the original system;
-C = MAX(P-M, 0) for the pertransposed system.
-C On exit, RO contains the last computed rank.
-C
-C SIGMA (input/output) INTEGER
-C On entry,
-C = 0 for the original system;
-C = M for the pertransposed system.
-C On exit, SIGMA contains the last computed value sigma in
-C the algorithm.
-C
-C SVLMAX (input) DOUBLE PRECISION
-C During each reduction step, the rank-revealing QR
-C factorization of a matrix stops when the estimated minimum
-C singular value is smaller than TOL * MAX(SVLMAX,EMSV),
-C where EMSV is the estimated maximum singular value.
-C SVLMAX >= 0.
-C
-C ABCD (input/output) DOUBLE PRECISION array, dimension
-C (LDABCD,M+N)
-C On entry, the leading (N+P)-by-(M+N) part of this array
-C must contain the compound input matrix of the system.
-C On exit, the leading (NU+MU)-by-(M+NU) part of this array
-C contains the reduced compound input matrix of the system.
-C
-C LDABCD INTEGER
-C The leading dimension of array ABCD.
-C LDABCD >= MAX(1,N+P).
-C
-C NINFZ (input/output) INTEGER
-C On entry, the currently computed number of infinite zeros.
-C It should be initialized to zero on the first call.
-C NINFZ >= 0.
-C On exit, the number of infinite zeros.
-C
-C INFZ (input/output) INTEGER array, dimension (N)
-C On entry, INFZ(i) must contain the current number of
-C infinite zeros of degree i, where i = 1,2,...,N, found in
-C the previous call(s) of the routine. It should be
-C initialized to zero on the first call.
-C On exit, INFZ(i) contains the number of infinite zeros of
-C degree i, where i = 1,2,...,N.
-C
-C KRONL (input/output) INTEGER array, dimension (N+1)
-C On entry, this array must contain the currently computed
-C left Kronecker (row) indices found in the previous call(s)
-C of the routine. It should be initialized to zero on the
-C first call.
-C On exit, the leading NKROL elements of this array contain
-C the left Kronecker (row) indices.
-C
-C MU (output) INTEGER
-C The normal rank of the transfer function matrix of the
-C original system.
-C
-C NU (output) INTEGER
-C The dimension of the reduced system matrix and the number
-C of (finite) invariant zeros if D' is invertible.
-C
-C NKROL (output) INTEGER
-C The number of left Kronecker indices.
-C
-C Tolerances
-C
-C TOL DOUBLE PRECISION
-C A tolerance used in rank decisions to determine the
-C effective rank, which is defined as the order of the
-C largest leading (or trailing) triangular submatrix in the
-C QR (or RQ) factorization with column (or row) pivoting
-C whose estimated condition number is less than 1/TOL.
-C NOTE that when SVLMAX > 0, the estimated ranks could be
-C less than those defined above (see SVLMAX).
-C
-C Workspace
-C
-C IWORK INTEGER array, dimension (MAX(M,P))
-C
-C DWORK DOUBLE PRECISION array, dimension (LDWORK)
-C On exit, if INFO = 0, DWORK(1) returns the optimal value
-C of LDWORK.
-C
-C LDWORK INTEGER
-C The length of the array DWORK.
-C LDWORK >= MAX( 1, MIN(P,M) + MAX(3*M-1,N),
-C MIN(P,N) + MAX(3*P-1,N+P,N+M) ).
-C For optimum performance LDWORK should be larger.
-C
-C If LDWORK = -1, then a workspace query is assumed;
-C the routine only calculates the optimal size of the
-C DWORK array, returns this value as the first entry of
-C the DWORK array, and no error message related to LDWORK
-C is issued by XERBLA.
-C
-C Error Indicator
-C
-C INFO INTEGER
-C = 0: successful exit;
-C < 0: if INFO = -i, the i-th argument had an illegal
-C value.
-C
-C REFERENCES
-C
-C [1] Svaricek, F.
-C Computation of the Structural Invariants of Linear
-C Multivariable Systems with an Extended Version of
-C the Program ZEROS.
-C System & Control Letters, 6, pp. 261-266, 1985.
-C
-C [2] Emami-Naeini, A. and Van Dooren, P.
-C Computation of Zeros of Linear Multivariable Systems.
-C Automatica, 18, pp. 415-430, 1982.
-C
-C NUMERICAL ASPECTS
-C
-C The algorithm is backward stable.
-C
-C CONTRIBUTOR
-C
-C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996.
-C Supersedes Release 2.0 routine AB08BZ by F. Svaricek.
-C
-C REVISIONS
-C
-C V. Sima, Oct. 1997; Feb. 1998, Jan. 2009, Apr. 2009.
-C A. Varga, May 1999; May 2001.
-C
-C KEYWORDS
-C
-C Generalized eigenvalue problem, Kronecker indices, multivariable
-C system, orthogonal transformation, structural invariant.
-C
-C ******************************************************************
-C
-C .. Parameters ..
- DOUBLE PRECISION ZERO
- PARAMETER ( ZERO = 0.0D0 )
-C .. Scalar Arguments ..
- INTEGER INFO, LDABCD, LDWORK, M, MU, N, NINFZ, NKROL,
- $ NU, P, RO, SIGMA
- DOUBLE PRECISION SVLMAX, TOL
-C .. Array Arguments ..
- INTEGER INFZ(*), IWORK(*), KRONL(*)
- DOUBLE PRECISION ABCD(LDABCD,*), DWORK(*)
-C .. Local Scalars ..
- LOGICAL LQUERY
- INTEGER I1, IK, IROW, ITAU, IZ, JWORK, MM1, MNTAU, MNU,
- $ MPM, NB, NP, RANK, RO1, TAU, WRKOPT
- DOUBLE PRECISION T
-C .. Local Arrays ..
- DOUBLE PRECISION SVAL(3)
-C .. External Functions ..
- INTEGER ILAENV
- EXTERNAL ILAENV
-C .. External Subroutines ..
- EXTERNAL DLAPMT, DLARFG, DLASET, DLATZM, DORMQR, DORMRQ,
- $ MB03OY, MB03PY, XERBLA
-C .. Intrinsic Functions ..
- INTRINSIC INT, MAX, MIN
-C .. Executable Statements ..
-C
- NP = N + P
- MPM = MIN( P, M )
- INFO = 0
- LQUERY = ( LDWORK.EQ.-1 )
-C
-C Test the input scalar arguments.
-C
- IF( N.LT.0 ) THEN
- INFO = -1
- ELSE IF( M.LT.0 ) THEN
- INFO = -2
- ELSE IF( P.LT.0 ) THEN
- INFO = -3
- ELSE IF( RO.NE.P .AND. RO.NE.MAX( P-M, 0 ) ) THEN
- INFO = -4
- ELSE IF( SIGMA.NE.0 .AND. SIGMA.NE.M ) THEN
- INFO = -5
- ELSE IF( SVLMAX.LT.ZERO ) THEN
- INFO = -6
- ELSE IF( LDABCD.LT.MAX( 1, NP ) ) THEN
- INFO = -8
- ELSE IF( NINFZ.LT.0 ) THEN
- INFO = -9
- ELSE
- JWORK = MAX( 1, MPM + MAX( 3*M - 1, N ),
- $ MIN( P, N ) + MAX( 3*P - 1, NP, N+M ) )
- IF( LQUERY ) THEN
- IF( M.GT.0 ) THEN
- NB = MIN( 64, ILAENV( 1, 'DORMQR', 'LT', P, N, MPM,
- $ -1 ) )
- WRKOPT = MAX( JWORK, MPM + MAX( 1, N )*NB )
- ELSE
- WRKOPT = JWORK
- END IF
- NB = MIN( 64, ILAENV( 1, 'DORMRQ', 'RT', NP, N, MIN( P, N ),
- $ -1 ) )
- WRKOPT = MAX( WRKOPT, MIN( P, N ) + MAX( 1, NP )*NB )
- NB = MIN( 64, ILAENV( 1, 'DORMRQ', 'LN', N, M+N,
- $ MIN( P, N ), -1 ) )
- WRKOPT = MAX( WRKOPT, MIN( P, N ) + MAX( 1, M+N )*NB )
- ELSE IF( LDWORK.LT.JWORK ) THEN
- INFO = -18
- END IF
- END IF
-C
- IF ( INFO.NE.0 ) THEN
-C
-C Error return.
-C
- CALL XERBLA( 'AB08NX', -INFO )
- RETURN
- ELSE IF( LQUERY ) THEN
- DWORK(1) = WRKOPT
- RETURN
- END IF
-C
- MU = P
- NU = N
-C
- IZ = 0
- IK = 1
- MM1 = M + 1
- ITAU = 1
- NKROL = 0
- WRKOPT = 1
-C
-C Main reduction loop:
-C
-C M NU M NU
-C NU [ B A ] NU [ B A ]
-C MU [ D C ] --> SIGMA [ RD C1 ] (SIGMA = rank(D) =
-C TAU [ 0 C2 ] row size of RD)
-C
-C M NU-RO RO
-C NU-RO [ B1 A11 A12 ]
-C --> RO [ B2 A21 A22 ] (RO = rank(C2) =
-C SIGMA [ RD C11 C12 ] col size of LC)
-C TAU [ 0 0 LC ]
-C
-C M NU-RO
-C NU-RO [ B1 A11 ] NU := NU - RO
-C [----------] MU := RO + SIGMA
-C --> RO [ B2 A21 ] D := [B2;RD]
-C SIGMA [ RD C11 ] C := [A21;C11]
-C
- 20 IF ( MU.EQ.0 )
- $ GO TO 80
-C
-C (Note: Comments in the code beginning "Workspace:" describe the
-C minimal amount of real workspace needed at that point in the
-C code, as well as the preferred amount for good performance.)
-C
- RO1 = RO
- MNU = M + NU
- IF ( M.GT.0 ) THEN
- IF ( SIGMA.NE.0 ) THEN
- IROW = NU + 1
-C
-C Compress rows of D. First exploit triangular shape.
-C Workspace: need M+N-1.
-C
- DO 40 I1 = 1, SIGMA
- CALL DLARFG( RO+1, ABCD(IROW,I1), ABCD(IROW+1,I1), 1, T )
- CALL DLATZM( 'L', RO+1, MNU-I1, ABCD(IROW+1,I1), 1, T,
- $ ABCD(IROW,I1+1), ABCD(IROW+1,I1+1), LDABCD,
- $ DWORK )
- IROW = IROW + 1
- 40 CONTINUE
- CALL DLASET( 'Lower', RO+SIGMA-1, SIGMA, ZERO, ZERO,
- $ ABCD(NU+2,1), LDABCD )
- END IF
-C
-C Continue with Householder with column pivoting.
-C
-C The rank of D is the number of (estimated) singular values
-C that are greater than TOL * MAX(SVLMAX,EMSV). This number
-C includes the singular values of the first SIGMA columns.
-C Integer workspace: need M;
-C Workspace: need min(RO1,M) + 3*M - 1. RO1 <= P.
-C
- IF ( SIGMA.LT.M ) THEN
- JWORK = ITAU + MIN( RO1, M )
- I1 = SIGMA + 1
- IROW = NU + I1
- CALL MB03OY( RO1, M-SIGMA, ABCD(IROW,I1), LDABCD, TOL,
- $ SVLMAX, RANK, SVAL, IWORK, DWORK(ITAU),
- $ DWORK(JWORK), INFO )
- WRKOPT = MAX( WRKOPT, JWORK + 3*M - 2 )
-C
-C Apply the column permutations to matrices B and part of D.
-C
- CALL DLAPMT( .TRUE., NU+SIGMA, M-SIGMA, ABCD(1,I1), LDABCD,
- $ IWORK )
-C
- IF ( RANK.GT.0 ) THEN
-C
-C Apply the Householder transformations to the submatrix C.
-C Workspace: need min(RO1,M) + NU;
-C prefer min(RO1,M) + NU*NB.
-C
- CALL DORMQR( 'Left', 'Transpose', RO1, NU, RANK,
- $ ABCD(IROW,I1), LDABCD, DWORK(ITAU),
- $ ABCD(IROW,MM1), LDABCD, DWORK(JWORK),
- $ LDWORK-JWORK+1, INFO )
- WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 )
- IF ( RO1.GT.1 )
- $ CALL DLASET( 'Lower', RO1-1, MIN( RO1-1, RANK ), ZERO,
- $ ZERO, ABCD(IROW+1,I1), LDABCD )
- RO1 = RO1 - RANK
- END IF
- END IF
- END IF
-C
- TAU = RO1
- SIGMA = MU - TAU
-C
-C Determination of the orders of the infinite zeros.
-C
- IF ( IZ.GT.0 ) THEN
- INFZ(IZ) = INFZ(IZ) + RO - TAU
- NINFZ = NINFZ + IZ*( RO - TAU )
- END IF
- IF ( RO1.EQ.0 )
- $ GO TO 80
- IZ = IZ + 1
-C
- IF ( NU.LE.0 ) THEN
- MU = SIGMA
- NU = 0
- RO = 0
- ELSE
-C
-C Compress the columns of C2 using RQ factorization with row
-C pivoting, P * C2 = R * Q.
-C
- I1 = NU + SIGMA + 1
- MNTAU = MIN( TAU, NU )
- JWORK = ITAU + MNTAU
-C
-C The rank of C2 is the number of (estimated) singular values
-C greater than TOL * MAX(SVLMAX,EMSV).
-C Integer Workspace: need TAU;
-C Workspace: need min(TAU,NU) + 3*TAU - 1.
-C
- CALL MB03PY( TAU, NU, ABCD(I1,MM1), LDABCD, TOL, SVLMAX, RANK,
- $ SVAL, IWORK, DWORK(ITAU), DWORK(JWORK), INFO )
- WRKOPT = MAX( WRKOPT, JWORK + 3*TAU - 1 )
- IF ( RANK.GT.0 ) THEN
- IROW = I1 + TAU - RANK
-C
-C Apply Q' to the first NU columns of [A; C1] from the right.
-C Workspace: need min(TAU,NU) + NU + SIGMA; SIGMA <= P;
-C prefer min(TAU,NU) + (NU + SIGMA)*NB.
-C
- CALL DORMRQ( 'Right', 'Transpose', I1-1, NU, RANK,
- $ ABCD(IROW,MM1), LDABCD, DWORK(MNTAU-RANK+1),
- $ ABCD(1,MM1), LDABCD, DWORK(JWORK),
- $ LDWORK-JWORK+1, INFO )
- WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 )
-C
-C Apply Q to the first NU rows and M + NU columns of [ B A ]
-C from the left.
-C Workspace: need min(TAU,NU) + M + NU;
-C prefer min(TAU,NU) + (M + NU)*NB.
-C
- CALL DORMRQ( 'Left', 'NoTranspose', NU, MNU, RANK,
- $ ABCD(IROW,MM1), LDABCD, DWORK(MNTAU-RANK+1),
- $ ABCD, LDABCD, DWORK(JWORK), LDWORK-JWORK+1,
- $ INFO )
- WRKOPT = MAX( WRKOPT, INT( DWORK(JWORK) ) + JWORK - 1 )
-C
- CALL DLASET( 'Full', RANK, NU-RANK, ZERO, ZERO,
- $ ABCD(IROW,MM1), LDABCD )
- IF ( RANK.GT.1 )
- $ CALL DLASET( 'Lower', RANK-1, RANK-1, ZERO, ZERO,
- $ ABCD(IROW+1,MM1+NU-RANK), LDABCD )
- END IF
-C
- RO = RANK
- END IF
-C
-C Determine the left Kronecker indices (row indices).
-C
- KRONL(IK) = KRONL(IK) + TAU - RO
- NKROL = NKROL + KRONL(IK)
- IK = IK + 1
-C
-C C and D are updated to [A21 ; C11] and [B2 ; RD].
-C
- NU = NU - RO
- MU = SIGMA + RO
- IF ( RO.NE.0 )
- $ GO TO 20
-C
- 80 CONTINUE
- DWORK(1) = WRKOPT
- RETURN
-C *** Last line of AB08NX ***
- END
Deleted: trunk/octave-forge/extra/control-devel/src/AB09AD.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB09AD.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB09AD.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,363 +0,0 @@
- SUBROUTINE AB09AD( DICO, JOB, EQUIL, ORDSEL, N, M, P, NR, A, LDA,
- $ B, LDB, C, LDC, HSV, TOL, IWORK, DWORK, LDWORK,
- $ IWARN, INFO )
-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 a reduced order model (Ar,Br,Cr) for a stable original
-C state-space representation (A,B,C) by using either the square-root
-C or the balancing-free square-root Balance & Truncate (B & T)
-C model reduction method.
-C
-C ARGUMENTS
-C
-C Mode Parameters
-C
-C DICO CHARACTER*1
-C Specifies the type of the original system as follows:
-C = 'C': continuous-time system;
-C = 'D': discrete-time system.
-C
-C JOB CHARACTER*1
-C Specifies the model reduction approach to be used
-C as follows:
-C = 'B': use the square-root Balance & Truncate method;
-C = 'N': use the balancing-free square-root
-C Balance & Truncate method.
-C
-C EQUIL CHARACTER*1
-C Specifies whether the user wishes to preliminarily
-C equilibrate the triplet (A,B,C) as follows:
-C = 'S': perform equilibration (scaling);
-C = 'N': do not perform equilibration.
-C
-C ORDSEL CHARACTER*1
-C Specifies the order selection method as follows:
-C = 'F': the resulting order NR is fixed;
-C = 'A': the resulting order NR is automatically determined
-C on basis of the given tolerance TOL.
-C
-C Input/Output Parameters
-C
-C N (input) INTEGER
-C The order of the original state-space representation, i.e.
-C the order of the matrix A. N >= 0.
-C
-C M (input) INTEGER
-C The number of system inputs. M >= 0.
-C
-C P (input) INTEGER
-C The number of system outputs. P >= 0.
-C
-C NR (input/output) INTEGER
-C On entry with ORDSEL = 'F', NR is the desired order of the
-C resulting reduced order system. 0 <= NR <= N.
-C On exit, if INFO = 0, NR is the order of the resulting
-C reduced order model. NR is set as follows:
-C if ORDSEL = 'F', NR is equal to MIN(NR,NMIN), where NR
-C is the desired order on entry and NMIN is the order of a
-C minimal realization of the given system; NMIN is
-C determined as the number of Hankel singular values greater
-C than N*EPS*HNORM(A,B,C), where EPS is the machine
-C precision (see LAPACK Library Routine DLAMCH) and
-C HNORM(A,B,C) is the Hankel norm of the system (computed
-C in HSV(1));
-C if ORDSEL = 'A', NR is equal to the number of Hankel
-C singular values greater than MAX(TOL,N*EPS*HNORM(A,B,C)).
-C
-C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
-C On entry, the leading N-by-N part of this array must
-C contain the state dynamics matrix A.
-C On exit, if INFO = 0, the leading NR-by-NR part of this
-C array contains the state dynamics matrix Ar of the reduced
-C order system.
-C
-C LDA INTEGER
-C The leading dimension of array A. LDA >= MAX(1,N).
-C
-C B (input/output) DOUBLE PRECISION array, dimension (LDB,M)
-C On entry, the leading N-by-M part of this array must
-C contain the original input/state matrix B.
-C On exit, if INFO = 0, the leading NR-by-M part of this
-C array contains the input/state matrix Br of the reduced
-C order system.
-C
-C LDB INTEGER
-C The leading dimension of array B. LDB >= MAX(1,N).
-C
-C C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
-C On entry, the leading P-by-N part of this array must
-C contain the original state/output matrix C.
-C On exit, if INFO = 0, the leading P-by-NR part of this
-C array contains the state/output matrix Cr of the reduced
-C order system.
-C
-C LDC INTEGER
-C The leading dimension of array C. LDC >= MAX(1,P).
-C
-C HSV (output) DOUBLE PRECISION array, dimension (N)
-C If INFO = 0, it contains the Hankel singular values of
-C the original system ordered decreasingly. HSV(1) is the
-C Hankel norm of the system.
-C
-C Tolerances
-C
-C TOL DOUBLE PRECISION
-C If ORDSEL = 'A', TOL contains the tolerance for
-C determining the order of reduced system.
-C For model reduction, the recommended value is
-C TOL = c*HNORM(A,B,C), where c is a constant in the
-C interval [0.00001,0.001], and HNORM(A,B,C) is the
-C Hankel-norm of the given system (computed in HSV(1)).
-C For computing a minimal realization, the recommended
-C value is TOL = N*EPS*HNORM(A,B,C), where EPS is the
-C machine precision (see LAPACK Library Routine DLAMCH).
-C This value is used by default if TOL <= 0 on entry.
-C If ORDSEL = 'F', the value of TOL is ignored.
-C
-C Workspace
-C
-C IWORK INTEGER array, dimension (LIWORK)
-C LIWORK = 0, if JOB = 'B';
-C LIWORK = N, if JOB = 'N'.
-C
-C DWORK DOUBLE PRECISION array, dimension (LDWORK)
-C On exit, if INFO = 0, DWORK(1) returns the optimal value
-C of LDWORK.
-C
-C LDWORK INTEGER
-C The length of the array DWORK.
-C LDWORK >= MAX(1,N*(2*N+MAX(N,M,P)+5)+N*(N+1)/2).
-C For optimum performance LDWORK should be larger.
-C
-C Warning Indicator
-C
-C IWARN INTEGER
-C = 0: no warning;
-C = 1: with ORDSEL = 'F', the selected order NR is greater
-C than the order of a minimal realization of the
-C given system. In this case, the resulting NR is
-C set automatically to a value corresponding to the
-C order of a minimal realization of the system.
-C
-C Error Indicator
-C
-C INFO INTEGER
-C = 0: successful exit;
-C < 0: if INFO = -i, the i-th argument had an illegal
-C value;
-C = 1: the reduction of A to the real Schur form failed;
-C = 2: the state matrix A is not stable (if DICO = 'C')
-C or not convergent (if DICO = 'D');
-C = 3: the computation of Hankel singular values failed.
-C
-C METHOD
-C
-C Let be the stable linear system
-C
-C d[x(t)] = Ax(t) + Bu(t)
-C y(t) = Cx(t) (1)
-C
-C where d[x(t)] is dx(t)/dt for a continuous-time system and x(t+1)
-C for a discrete-time system. The subroutine AB09AD determines for
-C the given system (1), the matrices of a reduced order system
-C
-C d[z(t)] = Ar*z(t) + Br*u(t)
-C yr(t) = Cr*z(t) (2)
-C
-C such that
-C
-C HSV(NR) <= INFNORM(G-Gr) <= 2*[HSV(NR+1) + ... + HSV(N)],
-C
-C where G and Gr are transfer-function matrices of the systems
-C (A,B,C) and (Ar,Br,Cr), respectively, and INFNORM(G) is the
-C infinity-norm of G.
-C
-C If JOB = 'B', the square-root Balance & Truncate method of [1]
-C is used and, for DICO = 'C', the resulting model is balanced.
-C By setting TOL <= 0, the routine can be used to compute balanced
-C minimal state-space realizations of stable systems.
-C
-C If JOB = 'N', the balancing-free square-root version of the
-C Balance & Truncate method [2] is used.
-C By setting TOL <= 0, the routine can be used to compute minimal
-C state-space realizations of stable systems.
-C
-C REFERENCES
-C
-C [1] Tombs M.S. and Postlethwaite I.
-C Truncated balanced realization of stable, non-minimal
-C state-space systems.
-C Int. J. Control, Vol. 46, pp. 1319-1330, 1987.
-C
-C [2] Varga A.
-C Efficient minimal realization procedure based on balancing.
-C Proc. of IMACS/IFAC Symp. MCTS, Lille, France, May 1991,
-C A. El Moudui, P. Borne, S. G. Tzafestas (Eds.),
-C Vol. 2, pp. 42-46.
-C
-C NUMERICAL ASPECTS
-C
-C The implemented methods rely on accuracy enhancing square-root or
-C balancing-free square-root techniques.
-C 3
-C The algorithms require less than 30N floating point operations.
-C
-C CONTRIBUTOR
-C
-C C. Oara and A. Varga, German Aerospace Center,
-C DLR Oberpfaffenhofen, March 1998.
-C Based on the RASP routines SRBT and SRBFT.
-C
-C REVISIONS
-C
-C May 2, 1998.
-C November 11, 1998, V. Sima, Research Institute for Informatics,
-C Bucharest.
-C
-C KEYWORDS
-C
-C Balancing, minimal state-space representation, model reduction,
-C multivariable system, state-space model.
-C
-C ******************************************************************
-C
-C .. Parameters ..
- DOUBLE PRECISION ONE, C100
- PARAMETER ( ONE = 1.0D0, C100 = 100.0D0 )
-C .. Scalar Arguments ..
- CHARACTER DICO, EQUIL, JOB, ORDSEL
- INTEGER INFO, IWARN, LDA, LDB, LDC, LDWORK, M, N, NR, P
- DOUBLE PRECISION TOL
-C .. Array Arguments ..
- INTEGER IWORK(*)
- DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), DWORK(*), HSV(*)
-C .. Local Scalars ..
- LOGICAL FIXORD
- INTEGER IERR, KI, KR, KT, KTI, KW, NN
- DOUBLE PRECISION MAXRED, WRKOPT
-C .. External Functions ..
- LOGICAL LSAME
- EXTERNAL LSAME
-C .. External Subroutines ..
- EXTERNAL AB09AX, TB01ID, TB01WD, XERBLA
-C .. Intrinsic Functions ..
- INTRINSIC DBLE, MAX, MIN
-C .. Executable Statements ..
-C
- INFO = 0
- IWARN = 0
- FIXORD = LSAME( ORDSEL, 'F' )
-C
-C Test the input scalar arguments.
-C
- IF( .NOT. ( LSAME( DICO, 'C' ) .OR. LSAME( DICO, 'D' ) ) ) THEN
- INFO = -1
- ELSE IF( .NOT. ( LSAME( JOB, 'B' ) .OR. LSAME( JOB, 'N' ) ) ) THEN
- INFO = -2
- ELSE IF( .NOT. ( LSAME( EQUIL, 'S' ) .OR.
- $ LSAME( EQUIL, 'N' ) ) ) THEN
- INFO = -3
- ELSE IF( .NOT. ( FIXORD .OR. LSAME( ORDSEL, 'A' ) ) ) THEN
- INFO = -4
- ELSE IF( N.LT.0 ) THEN
- INFO = -5
- ELSE IF( M.LT.0 ) THEN
- INFO = -6
- ELSE IF( P.LT.0 ) THEN
- INFO = -7
- ELSE IF( FIXORD .AND. ( NR.LT.0 .OR. NR.GT.N ) ) THEN
- INFO = -8
- ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
- INFO = -10
- ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
- INFO = -12
- ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
- INFO = -14
- ELSE IF( LDWORK.LT.MAX( 1, N*( 2*N + MAX( N, M, P ) + 5 ) +
- $ ( N*( N + 1 ) )/2 ) ) THEN
- INFO = -19
- END IF
-C
- IF( INFO.NE.0 ) THEN
-C
-C Error return.
-C
- CALL XERBLA( 'AB09AD', -INFO )
- RETURN
- END IF
-C
-C Quick return if possible.
-C
- IF( MIN( N, M, P ).EQ.0 .OR. ( FIXORD .AND. NR.EQ.0 ) ) THEN
- NR = 0
- DWORK(1) = ONE
- RETURN
- END IF
-C
-C Allocate working storage.
-C
- NN = N*N
- KT = 1
- KR = KT + NN
- KI = KR + N
- KW = KI + N
-C
- IF( LSAME( EQUIL, 'S' ) ) THEN
-C
-C Scale simultaneously the matrices A, B and C:
-C A <- inv(D)*A*D, B <- inv(D)*B and C <- C*D, where D is a
-C diagonal matrix.
-C
- MAXRED = C100
- CALL TB01ID( 'A', N, M, P, MAXRED, A, LDA, B, LDB, C, LDC,
- $ DWORK, INFO )
- END IF
-C
-C Reduce A to the real Schur form using an orthogonal similarity
-C transformation A <- T'*A*T and apply the transformation to
-C B and C: B <- T'*B and C <- C*T.
-C
- CALL TB01WD( N, M, P, A, LDA, B, LDB, C, LDC, DWORK(KT), N,
- $ DWORK(KR), DWORK(KI), DWORK(KW), LDWORK-KW+1, IERR )
- IF( IERR.NE.0 ) THEN
- INFO = 1
- RETURN
- END IF
-C
- WRKOPT = DWORK(KW) + DBLE( KW-1 )
- KTI = KT + NN
- KW = KTI + NN
-C
- CALL AB09AX( DICO, JOB, ORDSEL, N, M, P, NR, A, LDA, B, LDB, C,
- $ LDC, HSV, DWORK(KT), N, DWORK(KTI), N, TOL, IWORK,
- $ DWORK(KW), LDWORK-KW+1, IWARN, IERR )
-C
- IF( IERR.NE.0 ) THEN
- INFO = IERR + 1
- RETURN
- END IF
-C
- DWORK(1) = MAX( WRKOPT, DWORK(KW) + DBLE( KW-1 ) )
-C
- RETURN
-C *** Last line of AB09AD ***
- END
Deleted: trunk/octave-forge/extra/control-devel/src/AB09AX.f
===================================================================
--- trunk/octave-forge/extra/control-devel/src/AB09AX.f 2012-02-22 15:13:45 UTC (rev 9648)
+++ trunk/octave-forge/extra/control-devel/src/AB09AX.f 2012-02-22 16:37:14 UTC (rev 9649)
@@ -1,564 +0,0 @@
- SUBROUTINE AB09AX( DICO, JOB, ORDSEL, N, M, P, NR, A, LDA, B, LDB,
- $ C, LDC, HSV, T, LDT, TI, LDTI, TOL, IWORK,
- $ DWORK, LDWORK, IWARN, INFO )
-C
-C SLICOT RELEASE 5.0.
-C
-C Copyright (c) 2002-2009 NICONET e.V.
-C
@@ Diff output truncated at 100000 characters. @@
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|