Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[8381] trunk/octave-forge/main/control/devel
|
From: <par...@us...> - 2011-07-07 21:39:51
|
Revision: 8381
http://octave.svn.sourceforge.net/octave/?rev=8381&view=rev
Author: paramaniac
Date: 2011-07-07 21:39:42 +0000 (Thu, 07 Jul 2011)
Log Message:
-----------
control: add fortran files for model and controller reduction
Added Paths:
-----------
trunk/octave-forge/main/control/devel/slred/
trunk/octave-forge/main/control/devel/slred/reduction_methods.txt
trunk/octave-forge/main/control/devel/slred/slconred/
trunk/octave-forge/main/control/devel/slred/slconred/AB05PD.f
trunk/octave-forge/main/control/devel/slred/slconred/AB05QD.f
trunk/octave-forge/main/control/devel/slred/slconred/AB07ND.f
trunk/octave-forge/main/control/devel/slred/slconred/AB09AD.f
trunk/octave-forge/main/control/devel/slred/slconred/AB09AX.f
trunk/octave-forge/main/control/devel/slred/slconred/AB09BD.f
trunk/octave-forge/main/control/devel/slred/slconred/AB09BX.f
trunk/octave-forge/main/control/devel/slred/slconred/AB09DD.f
trunk/octave-forge/main/control/devel/slred/slconred/AB09IX.f
trunk/octave-forge/main/control/devel/slred/slconred/MA02AD.f
trunk/octave-forge/main/control/devel/slred/slconred/MA02DD.f
trunk/octave-forge/main/control/devel/slred/slconred/MA02GD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB01SD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB01WD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB01YD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB01ZD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB02UD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB03QD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB03QX.f
trunk/octave-forge/main/control/devel/slred/slconred/MB03QY.f
trunk/octave-forge/main/control/devel/slred/slconred/MB03UD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB04ND.f
trunk/octave-forge/main/control/devel/slred/slconred/MB04NY.f
trunk/octave-forge/main/control/devel/slred/slconred/MB04OD.f
trunk/octave-forge/main/control/devel/slred/slconred/MB04OY.f
trunk/octave-forge/main/control/devel/slred/slconred/SB03OD.f
trunk/octave-forge/main/control/devel/slred/slconred/SB03OR.f
trunk/octave-forge/main/control/devel/slred/slconred/SB03OT.f
trunk/octave-forge/main/control/devel/slred/slconred/SB03OU.f
trunk/octave-forge/main/control/devel/slred/slconred/SB03OV.f
trunk/octave-forge/main/control/devel/slred/slconred/SB03OY.f
trunk/octave-forge/main/control/devel/slred/slconred/SB04PX.f
trunk/octave-forge/main/control/devel/slred/slconred/SB08GD.f
trunk/octave-forge/main/control/devel/slred/slconred/SB08HD.f
trunk/octave-forge/main/control/devel/slred/slconred/SB16AD.f
trunk/octave-forge/main/control/devel/slred/slconred/SB16AY.f
trunk/octave-forge/main/control/devel/slred/slconred/SB16BD.f
trunk/octave-forge/main/control/devel/slred/slconred/SB16CD.f
trunk/octave-forge/main/control/devel/slred/slconred/SB16CY.f
trunk/octave-forge/main/control/devel/slred/slconred/TB01ID.f
trunk/octave-forge/main/control/devel/slred/slconred/TB01KD.f
trunk/octave-forge/main/control/devel/slred/slconred/TB01LD.f
trunk/octave-forge/main/control/devel/slred/slconred/TB01WD.f
trunk/octave-forge/main/control/devel/slred/slconred/makefile_slconred.m
trunk/octave-forge/main/control/devel/slred/slconred/select.f
trunk/octave-forge/main/control/devel/slred/slmodred/
trunk/octave-forge/main/control/devel/slred/slmodred/AB04MD.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB07MD.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB07ND.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB08MD.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB08NX.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09AX.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09CX.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09DD.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09HD.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09HY.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09ID.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09IX.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09IY.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09JD.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09JV.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09JW.f
trunk/octave-forge/main/control/devel/slred/slmodred/AB09JX.f
trunk/octave-forge/main/control/devel/slred/slmodred/AG07BD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MA02AD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MA02BD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MA02DD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB01PD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB01QD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB01SD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB01WD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB01YD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB01ZD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB03OY.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB03PY.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB03QD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB03QX.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB03QY.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB03UD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB04ND.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB04NY.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB04OD.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB04OX.f
trunk/octave-forge/main/control/devel/slred/slmodred/MB04OY.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB01FY.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB02MD.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB02MR.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB02MS.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB02MU.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB02MV.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB02MW.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB03OR.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB03OT.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB03OU.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB03OV.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB03OY.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB04PX.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB04PY.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB08CD.f
trunk/octave-forge/main/control/devel/slred/slmodred/SB08DD.f
trunk/octave-forge/main/control/devel/slred/slmodred/TB01ID.f
trunk/octave-forge/main/control/devel/slred/slmodred/TB01KD.f
trunk/octave-forge/main/control/devel/slred/slmodred/TB01LD.f
trunk/octave-forge/main/control/devel/slred/slmodred/TB01PD.f
trunk/octave-forge/main/control/devel/slred/slmodred/TB01UD.f
trunk/octave-forge/main/control/devel/slred/slmodred/TB01WD.f
trunk/octave-forge/main/control/devel/slred/slmodred/TB01XD.f
trunk/octave-forge/main/control/devel/slred/slmodred/delctg.f
trunk/octave-forge/main/control/devel/slred/slmodred/makefile_slmodred.m
trunk/octave-forge/main/control/devel/slred/slmodred/select.f
Added: trunk/octave-forge/main/control/devel/slred/reduction_methods.txt
===================================================================
--- trunk/octave-forge/main/control/devel/slred/reduction_methods.txt (rev 0)
+++ trunk/octave-forge/main/control/devel/slred/reduction_methods.txt 2011-07-07 21:39:42 UTC (rev 8381)
@@ -0,0 +1,57 @@
+Reduction Methods
+
+absolute/additive error
+ BTA balanced truncation
+ SPA singular perturbation approximation
+ HNA Hankel-norm approximation
+
+relative/multiplicative error
+ BST balanced stochastic truncation
+
+
+modred Model Reduction
+
+AB 09 ID - BTA SR
+ - BTA BFSR
+ - SPA SR
+ - SPA BFSR
+ FW BTA SR
+ FW BTA BFSR
+ FW SPA SR
+ FW SPA BFSR
+
+AB 09 JD - HNA
+ FW HNA
+
+AB 09 HD BST BTA SR
+ BST BTA BFSR
+ BST SPA SR
+ BST SPA BFSR
+
+
+conred Controller Reduction
+
+SB 16 AD - BTA SR
+ - BTA BFSR
+ - SPA SR
+ - SPA BFSR
+ FW BTA SR
+ FW BTA BFSR
+ FW SPA SR
+ FW SPA BFSR
+
+SB 16 BD CF BTA SR
+ CF BTA BFSR
+ CF SPA SR
+ CF SPA BFSR
+
+SB 16 CD FWCF BTA SR
+ FWCF BTA BFSR
+
+
+
+FW frequency-weighted
+CF coprime factorization
+
+SR square-root
+BFSR balancing-free square-root
\ No newline at end of file
Property changes on: trunk/octave-forge/main/control/devel/slred/slconred
___________________________________________________________________
Added: svn:ignore
+ *.oct
Added: trunk/octave-forge/main/control/devel/slred/slconred/AB05PD.f
===================================================================
--- trunk/octave-forge/main/control/devel/slred/slconred/AB05PD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/slred/slconred/AB05PD.f 2011-07-07 21:39:42 UTC (rev 8381)
@@ -0,0 +1,385 @@
+ 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
Added: trunk/octave-forge/main/control/devel/slred/slconred/AB05QD.f
===================================================================
--- trunk/octave-forge/main/control/devel/slred/slconred/AB05QD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/slred/slconred/AB05QD.f 2011-07-07 21:39:42 UTC (rev 8381)
@@ -0,0 +1,419 @@
+ 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
Added: trunk/octave-forge/main/control/devel/slred/slconred/AB07ND.f
===================================================================
--- trunk/octave-forge/main/control/devel/slred/slconred/AB07ND.f (rev 0)
+++ trunk/octave-forge/main/control/devel/slred/slconred/AB07ND.f 2011-07-07 21:39:42 UTC (rev 8381)
@@ -0,0 +1,303 @@
+ 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
Added: trunk/octave-forge/main/control/devel/slred/slconred/AB09AD.f
===================================================================
--- trunk/octave-forge/main/control/devel/slred/slconred/AB09AD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/slred/slconred/AB09AD.f 2011-07-07 21:39:42 UTC (rev 8381)
@@ -0,0 +1,363 @@
+ 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...
[truncated message content] |
