Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[7817] trunk/octave-forge/main/control/devel
|
From: <par...@us...> - 2010-10-08 19:53:23
|
Revision: 7817
http://octave.svn.sourceforge.net/octave/?rev=7817&view=rev
Author: paramaniac
Date: 2010-10-08 19:53:14 +0000 (Fri, 08 Oct 2010)
Log Message:
-----------
control: add draft code
Added Paths:
-----------
trunk/octave-forge/main/control/devel/ss2tf/
trunk/octave-forge/main/control/devel/ss2tf/MA02AD.f
trunk/octave-forge/main/control/devel/ss2tf/MB01PD.f
trunk/octave-forge/main/control/devel/ss2tf/MB01QD.f
trunk/octave-forge/main/control/devel/ss2tf/MB02RD.f
trunk/octave-forge/main/control/devel/ss2tf/MB02SD.f
trunk/octave-forge/main/control/devel/ss2tf/MC01PD.f
trunk/octave-forge/main/control/devel/ss2tf/MC01PY.f
trunk/octave-forge/main/control/devel/ss2tf/TB01ID.f
trunk/octave-forge/main/control/devel/ss2tf/TB01ZD.f
trunk/octave-forge/main/control/devel/ss2tf/TB04BD.dat
trunk/octave-forge/main/control/devel/ss2tf/TB04BD.f
trunk/octave-forge/main/control/devel/ss2tf/TB04BD.html
trunk/octave-forge/main/control/devel/ss2tf/TB04BD.res
trunk/octave-forge/main/control/devel/ss2tf/TB04BX.f
trunk/octave-forge/main/control/devel/ss2tf/common.cc
trunk/octave-forge/main/control/devel/ss2tf/makefile_ss2tf.m
trunk/octave-forge/main/control/devel/ss2tf/sltb04bd.cc
trunk/octave-forge/main/control/devel/ss2tf/test_ss2tf.m
Added: trunk/octave-forge/main/control/devel/ss2tf/MA02AD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/MA02AD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/MA02AD.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,108 @@
+ SUBROUTINE MA02AD( JOB, M, N, A, LDA, B, LDB )
+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 transpose all or part of a two-dimensional matrix A into
+C another matrix B.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C JOB CHARACTER*1
+C Specifies the part of the matrix A to be transposed into B
+C as follows:
+C = 'U': Upper triangular part;
+C = 'L': Lower triangular part;
+C Otherwise: All of the matrix A.
+C
+C Input/Output Parameters
+C
+C M (input) INTEGER
+C The number of rows of the matrix A. M >= 0.
+C
+C N (input) INTEGER
+C The number of columns of the matrix A. N >= 0.
+C
+C A (input) DOUBLE PRECISION array, dimension (LDA,N)
+C The m-by-n matrix A. If JOB = 'U', only the upper
+C triangle or trapezoid is accessed; if JOB = 'L', only the
+C lower triangle or trapezoid is accessed.
+C
+C LDA INTEGER
+C The leading dimension of the array A. LDA >= max(1,M).
+C
+C B (output) DOUBLE PRECISION array, dimension (LDB,M)
+C B = A' in the locations specified by JOB.
+C
+C LDB INTEGER
+C The leading dimension of the array B. LDB >= max(1,N).
+C
+C CONTRIBUTOR
+C
+C A. Varga, German Aerospace Center,
+C DLR Oberpfaffenhofen, March 1998.
+C Based on the RASP routine DMTRA.
+C
+C REVISIONS
+C
+C -
+C
+C ******************************************************************
+C
+C .. Scalar Arguments ..
+ CHARACTER JOB
+ INTEGER LDA, LDB, M, N
+C .. Array Arguments ..
+ DOUBLE PRECISION A(LDA,*), B(LDB,*)
+C .. Local Scalars ..
+ INTEGER I, J
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. Intrinsic Functions ..
+ INTRINSIC MIN
+C
+C .. Executable Statements ..
+C
+ IF( LSAME( JOB, 'U' ) ) THEN
+ DO 20 J = 1, N
+ DO 10 I = 1, MIN( J, M )
+ B(J,I) = A(I,J)
+ 10 CONTINUE
+ 20 CONTINUE
+ ELSE IF( LSAME( JOB, 'L' ) ) THEN
+ DO 40 J = 1, N
+ DO 30 I = J, M
+ B(J,I) = A(I,J)
+ 30 CONTINUE
+ 40 CONTINUE
+ ELSE
+ DO 60 J = 1, N
+ DO 50 I = 1, M
+ B(J,I) = A(I,J)
+ 50 CONTINUE
+ 60 CONTINUE
+ END IF
+C
+ RETURN
+C *** Last line of MA02AD ***
+ END
Added: trunk/octave-forge/main/control/devel/ss2tf/MB01PD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/MB01PD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/MB01PD.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,271 @@
+ SUBROUTINE MB01PD( SCUN, TYPE, M, N, KL, KU, ANRM, NBL, NROWS, A,
+ $ LDA, 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 scale a matrix or undo scaling. Scaling is performed, if
+C necessary, so that the matrix norm will be in a safe range of
+C representable numbers.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C SCUN CHARACTER*1
+C SCUN indicates the operation to be performed.
+C = 'S': scale the matrix.
+C = 'U': undo scaling of the matrix.
+C
+C TYPE CHARACTER*1
+C TYPE indicates the storage type of the input matrix.
+C = 'G': A is a full matrix.
+C = 'L': A is a (block) lower triangular matrix.
+C = 'U': A is an (block) upper triangular matrix.
+C = 'H': A is an (block) upper Hessenberg matrix.
+C = 'B': A is a symmetric band matrix with lower bandwidth
+C KL and upper bandwidth KU and with the only the
+C lower half stored.
+C = 'Q': A is a symmetric band matrix with lower bandwidth
+C KL and upper bandwidth KU and with the only the
+C upper half stored.
+C = 'Z': A is a band matrix with lower bandwidth KL and
+C upper bandwidth KU.
+C
+C Input/Output Parameters
+C
+C M (input) INTEGER
+C The number of rows of the matrix A. M >= 0.
+C
+C N (input) INTEGER
+C The number of columns of the matrix A. N >= 0.
+C
+C KL (input) INTEGER
+C The lower bandwidth of A. Referenced only if TYPE = 'B',
+C 'Q' or 'Z'.
+C
+C KU (input) INTEGER
+C The upper bandwidth of A. Referenced only if TYPE = 'B',
+C 'Q' or 'Z'.
+C
+C ANRM (input) DOUBLE PRECISION
+C The norm of the initial matrix A. ANRM >= 0.
+C When ANRM = 0 then an immediate return is effected.
+C ANRM should be preserved between the call of the routine
+C with SCUN = 'S' and the corresponding one with SCUN = 'U'.
+C
+C NBL (input) INTEGER
+C The number of diagonal blocks of the matrix A, if it has a
+C block structure. To specify that matrix A has no block
+C structure, set NBL = 0. NBL >= 0.
+C
+C NROWS (input) INTEGER array, dimension max(1,NBL)
+C NROWS(i) contains the number of rows and columns of the
+C i-th diagonal block of matrix A. The sum of the values
+C NROWS(i), for i = 1: NBL, should be equal to min(M,N).
+C The elements of the array NROWS are not referenced if
+C NBL = 0.
+C
+C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+C On entry, the leading M by N part of this array must
+C contain the matrix to be scaled/unscaled.
+C On exit, the leading M by N part of A will contain
+C the modified matrix.
+C The storage mode of A is specified by TYPE.
+C
+C LDA (input) INTEGER
+C The leading dimension of the array A. LDA >= max(1,M).
+C
+C Error Indicator
+C
+C INFO (output) INTEGER
+C = 0: successful exit
+C < 0: if INFO = -i, the i-th argument had an illegal
+C value.
+C
+C METHOD
+C
+C Denote by ANRM the norm of the matrix, and by SMLNUM and BIGNUM,
+C two positive numbers near the smallest and largest safely
+C representable numbers, respectively. The matrix is scaled, if
+C needed, such that the norm of the result is in the range
+C [SMLNUM, BIGNUM]. The scaling factor is represented as a ratio
+C of two numbers, one of them being ANRM, and the other one either
+C SMLNUM or BIGNUM, depending on ANRM being less than SMLNUM or
+C larger than BIGNUM, respectively. For undoing the scaling, the
+C norm is again compared with SMLNUM or BIGNUM, and the reciprocal
+C of the previous scaling factor is used.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996.
+C
+C REVISIONS
+C
+C Oct. 2001, V. Sima, Research Institute for Informatics, Bucharest.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C .. Scalar Arguments ..
+ CHARACTER SCUN, TYPE
+ INTEGER INFO, KL, KU, LDA, M, MN, N, NBL
+ DOUBLE PRECISION ANRM
+C .. Array Arguments ..
+ INTEGER NROWS ( * )
+ DOUBLE PRECISION A( LDA, * )
+C .. Local Scalars ..
+ LOGICAL FIRST, LSCALE
+ INTEGER I, ISUM, ITYPE
+ DOUBLE PRECISION BIGNUM, SMLNUM
+C .. External Functions ..
+ LOGICAL LSAME
+ DOUBLE PRECISION DLAMCH
+ EXTERNAL DLAMCH, LSAME
+C ..
+C .. External Subroutines ..
+ EXTERNAL DLABAD, MB01QD, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC MAX, MIN
+C .. Save statement ..
+ SAVE BIGNUM, FIRST, SMLNUM
+C .. Data statements ..
+ DATA FIRST/.TRUE./
+C ..
+C .. Executable Statements ..
+C
+C Test the input scalar arguments.
+C
+ INFO = 0
+ LSCALE = LSAME( SCUN, 'S' )
+ IF( LSAME( TYPE, 'G' ) ) THEN
+ ITYPE = 0
+ ELSE IF( LSAME( TYPE, 'L' ) ) THEN
+ ITYPE = 1
+ ELSE IF( LSAME( TYPE, 'U' ) ) THEN
+ ITYPE = 2
+ ELSE IF( LSAME( TYPE, 'H' ) ) THEN
+ ITYPE = 3
+ ELSE IF( LSAME( TYPE, 'B' ) ) THEN
+ ITYPE = 4
+ ELSE IF( LSAME( TYPE, 'Q' ) ) THEN
+ ITYPE = 5
+ ELSE IF( LSAME( TYPE, 'Z' ) ) THEN
+ ITYPE = 6
+ ELSE
+ ITYPE = -1
+ END IF
+C
+ MN = MIN( M, N )
+C
+ ISUM = 0
+ IF( NBL.GT.0 ) THEN
+ DO 10 I = 1, NBL
+ ISUM = ISUM + NROWS(I)
+ 10 CONTINUE
+ END IF
+C
+ IF( .NOT.LSCALE .AND. .NOT.LSAME( SCUN, 'U' ) ) THEN
+ INFO = -1
+ ELSE IF( ITYPE.EQ.-1 ) THEN
+ INFO = -2
+ ELSE IF( M.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 .OR.
+ $ ( ( ITYPE.EQ.4 .OR. ITYPE.EQ.5 ) .AND. N.NE.M ) ) THEN
+ INFO = -4
+ ELSE IF( ANRM.LT.ZERO ) THEN
+ INFO = -7
+ ELSE IF( NBL.LT.0 ) THEN
+ INFO = -8
+ ELSE IF( NBL.GT.0 .AND. ISUM.NE.MN ) THEN
+ INFO = -9
+ ELSE IF( ITYPE.LE.3 .AND. LDA.LT.MAX( 1, M ) ) THEN
+ INFO = -11
+ ELSE IF( ITYPE.GE.4 ) THEN
+ IF( KL.LT.0 .OR. KL.GT.MAX( M-1, 0 ) ) THEN
+ INFO = -5
+ ELSE IF( KU.LT.0 .OR. KU.GT.MAX( N-1, 0 ) .OR.
+ $ ( ( ITYPE.EQ.4 .OR. ITYPE.EQ.5 ) .AND. KL.NE.KU ) )
+ $ THEN
+ INFO = -6
+ ELSE IF( ( ITYPE.EQ.4 .AND. LDA.LT.KL+1 ) .OR.
+ $ ( ITYPE.EQ.5 .AND. LDA.LT.KU+1 ) .OR.
+ $ ( ITYPE.EQ.6 .AND. LDA.LT.2*KL+KU+1 ) ) THEN
+ INFO = -11
+ END IF
+ END IF
+C
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'MB01PD', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF( MN.EQ.0 .OR. ANRM.EQ.ZERO )
+ $ RETURN
+C
+ IF ( FIRST ) THEN
+C
+C Get machine parameters.
+C
+ SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'P' )
+ BIGNUM = ONE / SMLNUM
+ CALL DLABAD( SMLNUM, BIGNUM )
+ FIRST = .FALSE.
+ END IF
+C
+ IF ( LSCALE ) THEN
+C
+C Scale A, if its norm is outside range [SMLNUM,BIGNUM].
+C
+ IF( ANRM.LT.SMLNUM ) THEN
+C
+C Scale matrix norm up to SMLNUM.
+C
+ CALL MB01QD( TYPE, M, N, KL, KU, ANRM, SMLNUM, NBL, NROWS,
+ $ A, LDA, INFO )
+ ELSE IF( ANRM.GT.BIGNUM ) THEN
+C
+C Scale matrix norm down to BIGNUM.
+C
+ CALL MB01QD( TYPE, M, N, KL, KU, ANRM, BIGNUM, NBL, NROWS,
+ $ A, LDA, INFO )
+ END IF
+C
+ ELSE
+C
+C Undo scaling.
+C
+ IF( ANRM.LT.SMLNUM ) THEN
+ CALL MB01QD( TYPE, M, N, KL, KU, SMLNUM, ANRM, NBL, NROWS,
+ $ A, LDA, INFO )
+ ELSE IF( ANRM.GT.BIGNUM ) THEN
+ CALL MB01QD( TYPE, M, N, KL, KU, BIGNUM, ANRM, NBL, NROWS,
+ $ A, LDA, INFO )
+ END IF
+ END IF
+C
+ RETURN
+C *** Last line of MB01PD ***
+ END
Added: trunk/octave-forge/main/control/devel/ss2tf/MB01QD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/MB01QD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/MB01QD.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,334 @@
+ SUBROUTINE MB01QD( TYPE, M, N, KL, KU, CFROM, CTO, NBL, NROWS, A,
+ $ LDA, 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 multiply the M by N real matrix A by the real scalar CTO/CFROM.
+C This is done without over/underflow as long as the final result
+C CTO*A(I,J)/CFROM does not over/underflow. TYPE specifies that
+C A may be full, (block) upper triangular, (block) lower triangular,
+C (block) upper Hessenberg, or banded.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C TYPE CHARACTER*1
+C TYPE indices the storage type of the input matrix.
+C = 'G': A is a full matrix.
+C = 'L': A is a (block) lower triangular matrix.
+C = 'U': A is a (block) upper triangular matrix.
+C = 'H': A is a (block) upper Hessenberg matrix.
+C = 'B': A is a symmetric band matrix with lower bandwidth
+C KL and upper bandwidth KU and with the only the
+C lower half stored.
+C = 'Q': A is a symmetric band matrix with lower bandwidth
+C KL and upper bandwidth KU and with the only the
+C upper half stored.
+C = 'Z': A is a band matrix with lower bandwidth KL and
+C upper bandwidth KU.
+C
+C Input/Output Parameters
+C
+C M (input) INTEGER
+C The number of rows of the matrix A. M >= 0.
+C
+C N (input) INTEGER
+C The number of columns of the matrix A. N >= 0.
+C
+C KL (input) INTEGER
+C The lower bandwidth of A. Referenced only if TYPE = 'B',
+C 'Q' or 'Z'.
+C
+C KU (input) INTEGER
+C The upper bandwidth of A. Referenced only if TYPE = 'B',
+C 'Q' or 'Z'.
+C
+C CFROM (input) DOUBLE PRECISION
+C CTO (input) DOUBLE PRECISION
+C The matrix A is multiplied by CTO/CFROM. A(I,J) is
+C computed without over/underflow if the final result
+C CTO*A(I,J)/CFROM can be represented without over/
+C underflow. CFROM must be nonzero.
+C
+C NBL (input) INTEGER
+C The number of diagonal blocks of the matrix A, if it has a
+C block structure. To specify that matrix A has no block
+C structure, set NBL = 0. NBL >= 0.
+C
+C NROWS (input) INTEGER array, dimension max(1,NBL)
+C NROWS(i) contains the number of rows and columns of the
+C i-th diagonal block of matrix A. The sum of the values
+C NROWS(i), for i = 1: NBL, should be equal to min(M,N).
+C The array NROWS is not referenced if NBL = 0.
+C
+C A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+C The matrix to be multiplied by CTO/CFROM. See TYPE for
+C the storage type.
+C
+C LDA (input) INTEGER
+C The leading dimension of the array A. LDA >= max(1,M).
+C
+C Error Indicator
+C
+C INFO INTEGER
+C Not used in this implementation.
+C
+C METHOD
+C
+C Matrix A is multiplied by the real scalar CTO/CFROM, taking into
+C account the specified storage mode of the matrix.
+C MB01QD is a version of the LAPACK routine DLASCL, modified for
+C dealing with block triangular, or block Hessenberg matrices.
+C For efficiency, no tests of the input scalar parameters are
+C performed.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, Nov. 1996.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C ..
+C .. Scalar Arguments ..
+ CHARACTER TYPE
+ INTEGER INFO, KL, KU, LDA, M, N, NBL
+ DOUBLE PRECISION CFROM, CTO
+C ..
+C .. Array Arguments ..
+ INTEGER NROWS ( * )
+ DOUBLE PRECISION A( LDA, * )
+C ..
+C .. Local Scalars ..
+ LOGICAL DONE, NOBLC
+ INTEGER I, IFIN, ITYPE, J, JFIN, JINI, K, K1, K2, K3,
+ $ K4
+ DOUBLE PRECISION BIGNUM, CFROM1, CFROMC, CTO1, CTOC, MUL, SMLNUM
+C ..
+C .. External Functions ..
+ LOGICAL LSAME
+ DOUBLE PRECISION DLAMCH
+ EXTERNAL LSAME, DLAMCH
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX, MIN
+C ..
+C .. Executable Statements ..
+C
+ IF( LSAME( TYPE, 'G' ) ) THEN
+ ITYPE = 0
+ ELSE IF( LSAME( TYPE, 'L' ) ) THEN
+ ITYPE = 1
+ ELSE IF( LSAME( TYPE, 'U' ) ) THEN
+ ITYPE = 2
+ ELSE IF( LSAME( TYPE, 'H' ) ) THEN
+ ITYPE = 3
+ ELSE IF( LSAME( TYPE, 'B' ) ) THEN
+ ITYPE = 4
+ ELSE IF( LSAME( TYPE, 'Q' ) ) THEN
+ ITYPE = 5
+ ELSE
+ ITYPE = 6
+ END IF
+C
+C Quick return if possible.
+C
+ IF( MIN( M, N ).EQ.0 )
+ $ RETURN
+C
+C Get machine parameters.
+C
+ SMLNUM = DLAMCH( 'S' )
+ BIGNUM = ONE / SMLNUM
+C
+ CFROMC = CFROM
+ CTOC = CTO
+C
+ 10 CONTINUE
+ CFROM1 = CFROMC*SMLNUM
+ CTO1 = CTOC / BIGNUM
+ IF( ABS( CFROM1 ).GT.ABS( CTOC ) .AND. CTOC.NE.ZERO ) THEN
+ MUL = SMLNUM
+ DONE = .FALSE.
+ CFROMC = CFROM1
+ ELSE IF( ABS( CTO1 ).GT.ABS( CFROMC ) ) THEN
+ MUL = BIGNUM
+ DONE = .FALSE.
+ CTOC = CTO1
+ ELSE
+ MUL = CTOC / CFROMC
+ DONE = .TRUE.
+ END IF
+C
+ NOBLC = NBL.EQ.0
+C
+ IF( ITYPE.EQ.0 ) THEN
+C
+C Full matrix
+C
+ DO 30 J = 1, N
+ DO 20 I = 1, M
+ A( I, J ) = A( I, J )*MUL
+ 20 CONTINUE
+ 30 CONTINUE
+C
+ ELSE IF( ITYPE.EQ.1 ) THEN
+C
+ IF ( NOBLC ) THEN
+C
+C Lower triangular matrix
+C
+ DO 50 J = 1, N
+ DO 40 I = J, M
+ A( I, J ) = A( I, J )*MUL
+ 40 CONTINUE
+ 50 CONTINUE
+C
+ ELSE
+C
+C Block lower triangular matrix
+C
+ JFIN = 0
+ DO 80 K = 1, NBL
+ JINI = JFIN + 1
+ JFIN = JFIN + NROWS( K )
+ DO 70 J = JINI, JFIN
+ DO 60 I = JINI, M
+ A( I, J ) = A( I, J )*MUL
+ 60 CONTINUE
+ 70 CONTINUE
+ 80 CONTINUE
+ END IF
+C
+ ELSE IF( ITYPE.EQ.2 ) THEN
+C
+ IF ( NOBLC ) THEN
+C
+C Upper triangular matrix
+C
+ DO 100 J = 1, N
+ DO 90 I = 1, MIN( J, M )
+ A( I, J ) = A( I, J )*MUL
+ 90 CONTINUE
+ 100 CONTINUE
+C
+ ELSE
+C
+C Block upper triangular matrix
+C
+ JFIN = 0
+ DO 130 K = 1, NBL
+ JINI = JFIN + 1
+ JFIN = JFIN + NROWS( K )
+ IF ( K.EQ.NBL ) JFIN = N
+ DO 120 J = JINI, JFIN
+ DO 110 I = 1, MIN( JFIN, M )
+ A( I, J ) = A( I, J )*MUL
+ 110 CONTINUE
+ 120 CONTINUE
+ 130 CONTINUE
+ END IF
+C
+ ELSE IF( ITYPE.EQ.3 ) THEN
+C
+ IF ( NOBLC ) THEN
+C
+C Upper Hessenberg matrix
+C
+ DO 150 J = 1, N
+ DO 140 I = 1, MIN( J+1, M )
+ A( I, J ) = A( I, J )*MUL
+ 140 CONTINUE
+ 150 CONTINUE
+C
+ ELSE
+C
+C Block upper Hessenberg matrix
+C
+ JFIN = 0
+ DO 180 K = 1, NBL
+ JINI = JFIN + 1
+ JFIN = JFIN + NROWS( K )
+C
+ IF ( K.EQ.NBL ) THEN
+ JFIN = N
+ IFIN = N
+ ELSE
+ IFIN = JFIN + NROWS( K+1 )
+ END IF
+C
+ DO 170 J = JINI, JFIN
+ DO 160 I = 1, MIN( IFIN, M )
+ A( I, J ) = A( I, J )*MUL
+ 160 CONTINUE
+ 170 CONTINUE
+ 180 CONTINUE
+ END IF
+C
+ ELSE IF( ITYPE.EQ.4 ) THEN
+C
+C Lower half of a symmetric band matrix
+C
+ K3 = KL + 1
+ K4 = N + 1
+ DO 200 J = 1, N
+ DO 190 I = 1, MIN( K3, K4-J )
+ A( I, J ) = A( I, J )*MUL
+ 190 CONTINUE
+ 200 CONTINUE
+C
+ ELSE IF( ITYPE.EQ.5 ) THEN
+C
+C Upper half of a symmetric band matrix
+C
+ K1 = KU + 2
+ K3 = KU + 1
+ DO 220 J = 1, N
+ DO 210 I = MAX( K1-J, 1 ), K3
+ A( I, J ) = A( I, J )*MUL
+ 210 CONTINUE
+ 220 CONTINUE
+C
+ ELSE IF( ITYPE.EQ.6 ) THEN
+C
+C Band matrix
+C
+ K1 = KL + KU + 2
+ K2 = KL + 1
+ K3 = 2*KL + KU + 1
+ K4 = KL + KU + 1 + M
+ DO 240 J = 1, N
+ DO 230 I = MAX( K1-J, K2 ), MIN( K3, K4-J )
+ A( I, J ) = A( I, J )*MUL
+ 230 CONTINUE
+ 240 CONTINUE
+C
+ END IF
+C
+ IF( .NOT.DONE )
+ $ GO TO 10
+C
+ RETURN
+C *** Last line of MB01QD ***
+ END
Added: trunk/octave-forge/main/control/devel/ss2tf/MB02RD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/MB02RD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/MB02RD.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,197 @@
+ SUBROUTINE MB02RD( TRANS, N, NRHS, H, LDH, IPIV, B, LDB, 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 solve a system of linear equations
+C H * X = B or H' * X = B
+C with an upper Hessenberg N-by-N matrix H using the LU
+C factorization computed by MB02SD.
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C TRANS CHARACTER*1
+C Specifies the form of the system of equations:
+C = 'N': H * X = B (No transpose)
+C = 'T': H'* X = B (Transpose)
+C = 'C': H'* X = B (Conjugate transpose = Transpose)
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the matrix H. N >= 0.
+C
+C NRHS (input) INTEGER
+C The number of right hand sides, i.e., the number of
+C columns of the matrix B. NRHS >= 0.
+C
+C H (input) DOUBLE PRECISION array, dimension (LDH,N)
+C The factors L and U from the factorization H = P*L*U
+C as computed by MB02SD.
+C
+C LDH INTEGER
+C The leading dimension of the array H. LDH >= max(1,N).
+C
+C IPIV (input) INTEGER array, dimension (N)
+C The pivot indices from MB02SD; for 1<=i<=N, row i of the
+C matrix was interchanged with row IPIV(i).
+C
+C B (input/output) DOUBLE PRECISION array, dimension
+C (LDB,NRHS)
+C On entry, the right hand side matrix B.
+C On exit, the solution matrix X.
+C
+C LDB INTEGER
+C The leading dimension of the array B. LDB >= max(1,N).
+C
+C Error Indicator
+C
+C INFO (output) 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 uses the factorization
+C H = P * L * U
+C where P is a permutation matrix, L is lower triangular with unit
+C diagonal elements (and one nonzero subdiagonal), and U is upper
+C triangular.
+C
+C REFERENCES
+C
+C -
+C
+C NUMERICAL ASPECTS
+C 2
+C The algorithm requires 0( N x NRHS ) operations.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, June 1998.
+C
+C REVISIONS
+C
+C -
+C
+C KEYWORDS
+C
+C Hessenberg form, matrix algebra.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ONE
+ PARAMETER ( ONE = 1.0D+0 )
+C .. Scalar Arguments ..
+ CHARACTER TRANS
+ INTEGER INFO, LDB, LDH, N, NRHS
+C ..
+C .. Array Arguments ..
+ INTEGER IPIV( * )
+ DOUBLE PRECISION B( LDB, * ), H( LDH, * )
+C .. Local Scalars ..
+ LOGICAL NOTRAN
+ INTEGER J, JP
+C .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+C .. External Subroutines ..
+ EXTERNAL DAXPY, DSWAP, DTRSM, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC MAX
+C .. Executable Statements ..
+C
+C Test the input parameters.
+C
+ INFO = 0
+ NOTRAN = LSAME( TRANS, 'N' )
+ IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) .AND. .NOT.
+ $ LSAME( TRANS, 'C' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( NRHS.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
+ INFO = -5
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -8
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'MB02RD', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF( N.EQ.0 .OR. NRHS.EQ.0 )
+ $ RETURN
+C
+ IF( NOTRAN ) THEN
+C
+C Solve H * X = B.
+C
+C Solve L * X = B, overwriting B with X.
+C
+C L is represented as a product of permutations and unit lower
+C triangular matrices L = P(1) * L(1) * ... * P(n-1) * L(n-1),
+C where each transformation L(i) is a rank-one modification of
+C the identity matrix.
+C
+ DO 10 J = 1, N - 1
+ JP = IPIV( J )
+ IF( JP.NE.J )
+ $ CALL DSWAP( NRHS, B( JP, 1 ), LDB, B( J, 1 ), LDB )
+ CALL DAXPY( NRHS, -H( J+1, J ), B( J, 1 ), LDB, B( J+1, 1 ),
+ $ LDB )
+ 10 CONTINUE
+C
+C Solve U * X = B, overwriting B with X.
+C
+ CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
+ $ NRHS, ONE, H, LDH, B, LDB )
+C
+ ELSE
+C
+C Solve H' * X = B.
+C
+C Solve U' * X = B, overwriting B with X.
+C
+ CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS,
+ $ ONE, H, LDH, B, LDB )
+C
+C Solve L' * X = B, overwriting B with X.
+C
+ DO 20 J = N - 1, 1, -1
+ CALL DAXPY( NRHS, -H( J+1, J ), B( J+1, 1 ), LDB, B( J, 1 ),
+ $ LDB )
+ JP = IPIV( J )
+ IF( JP.NE.J )
+ $ CALL DSWAP( NRHS, B( JP, 1 ), LDB, B( J, 1 ), LDB )
+ 20 CONTINUE
+ END IF
+C
+ RETURN
+C *** Last line of MB02RD ***
+ END
Added: trunk/octave-forge/main/control/devel/ss2tf/MB02SD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/MB02SD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/MB02SD.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,164 @@
+ SUBROUTINE MB02SD( N, H, LDH, IPIV, 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 an LU factorization of an n-by-n upper Hessenberg
+C matrix H using partial pivoting with row interchanges.
+C
+C ARGUMENTS
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the matrix H. N >= 0.
+C
+C H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
+C On entry, the n-by-n upper Hessenberg matrix to be
+C factored.
+C On exit, the factors L and U from the factorization
+C H = P*L*U; the unit diagonal elements of L are not stored,
+C and L is lower bidiagonal.
+C
+C LDH INTEGER
+C The leading dimension of the array H. LDH >= max(1,N).
+C
+C IPIV (output) INTEGER array, dimension (N)
+C The pivot indices; for 1 <= i <= N, row i of the matrix
+C was interchanged with row IPIV(i).
+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 > 0: if INFO = i, U(i,i) is exactly zero. The
+C factorization has been completed, but the factor U
+C is exactly singular, and division by zero will occur
+C if it is used to solve a system of equations.
+C
+C METHOD
+C
+C The factorization has the form
+C H = P * L * U
+C where P is a permutation matrix, L is lower triangular with unit
+C diagonal elements (and one nonzero subdiagonal), and U is upper
+C triangular.
+C
+C This is the right-looking Level 1 BLAS version of the algorithm
+C (adapted after DGETF2).
+C
+C REFERENCES
+C
+C -
+C
+C NUMERICAL ASPECTS
+C 2
+C The algorithm requires 0( N ) operations.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Katholieke Univ. Leuven, Belgium, June 1998.
+C
+C REVISIONS
+C
+C V. Sima, Research Institute for Informatics, Bucharest, Oct. 2000,
+C Jan. 2005.
+C
+C KEYWORDS
+C
+C Hessenberg form, matrix algebra.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO
+ PARAMETER ( ZERO = 0.0D+0 )
+C .. Scalar Arguments ..
+ INTEGER INFO, LDH, N
+C .. Array Arguments ..
+ INTEGER IPIV(*)
+ DOUBLE PRECISION H(LDH,*)
+C .. Local Scalars ..
+ INTEGER J, JP
+C .. External Subroutines ..
+ EXTERNAL DAXPY, DSWAP, XERBLA
+C .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX
+C ..
+C .. Executable Statements ..
+C
+C Check the scalar input parameters.
+C
+ INFO = 0
+ IF( N.LT.0 ) THEN
+ INFO = -1
+ ELSE IF( LDH.LT.MAX( 1, N ) ) THEN
+ INFO = -3
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'MB02SD', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ IF( N.EQ.0 )
+ $ RETURN
+C
+ DO 10 J = 1, N
+C
+C Find pivot and test for singularity.
+C
+ JP = J
+ IF ( J.LT.N ) THEN
+ IF ( ABS( H( J+1, J ) ).GT.ABS( H( J, J ) ) )
+ $ JP = J + 1
+ END IF
+ IPIV( J ) = JP
+ IF( H( JP, J ).NE.ZERO ) THEN
+C
+C Apply the interchange to columns J:N.
+C
+ IF( JP.NE.J )
+ $ CALL DSWAP( N-J+1, H( J, J ), LDH, H( JP, J ), LDH )
+C
+C Compute element J+1 of J-th column.
+C
+ IF( J.LT.N )
+ $ H( J+1, J ) = H( J+1, J )/H( J, J )
+C
+ ELSE IF( INFO.EQ.0 ) THEN
+C
+ INFO = J
+ END IF
+C
+ IF( J.LT.N ) THEN
+C
+C Update trailing submatrix.
+C
+ CALL DAXPY( N-J, -H( J+1, J ), H( J, J+1 ), LDH,
+ $ H( J+1, J+1 ), LDH )
+ END IF
+ 10 CONTINUE
+ RETURN
+C *** Last line of MB02SD ***
+ END
Added: trunk/octave-forge/main/control/devel/ss2tf/MC01PD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/MC01PD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/MC01PD.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,159 @@
+ SUBROUTINE MC01PD( K, REZ, IMZ, P, DWORK, 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 coefficients of a real polynomial P(x) from its
+C zeros.
+C
+C ARGUMENTS
+C
+C Input/Output Parameters
+C
+C K (input) INTEGER
+C The number of zeros (and hence the degree) of P(x).
+C K >= 0.
+C
+C REZ (input) DOUBLE PRECISION array, dimension (K)
+C IMZ (input) DOUBLE PRECISION array, dimension (K)
+C The real and imaginary parts of the i-th zero of P(x)
+C must be stored in REZ(i) and IMZ(i), respectively, where
+C i = 1, 2, ..., K. The zeros may be supplied in any order,
+C except that complex conjugate zeros must appear
+C consecutively.
+C
+C P (output) DOUBLE PRECISION array, dimension (K+1)
+C This array contains the coefficients of P(x) in increasing
+C powers of x. If K = 0, then P(1) is set to one.
+C
+C Workspace
+C
+C DWORK DOUBLE PRECISION array, dimension (K+1)
+C If K = 0, this array is not referenced.
+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 > 0: if INFO = i, (REZ(i),IMZ(i)) is a complex zero but
+C (REZ(i-1),IMZ(i-1)) is not its conjugate.
+C
+C METHOD
+C
+C The routine computes the coefficients of the real K-th degree
+C polynomial P(x) as
+C
+C P(x) = (x - r(1)) * (x - r(2)) * ... * (x - r(K))
+C
+C where r(i) = (REZ(i),IMZ(i)).
+C
+C Note that REZ(i) = REZ(j) and IMZ(i) = -IMZ(j) if r(i) and r(j)
+C form a complex conjugate pair (where i <> j), and that IMZ(i) = 0
+C if r(i) is real.
+C
+C NUMERICAL ASPECTS
+C
+C None.
+C
+C CONTRIBUTOR
+C
+C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Mar. 1997.
+C Supersedes Release 2.0 routine MC01DD by A.J. Geurts.
+C
+C REVISIONS
+C
+C V. Sima, May 2002.
+C
+C KEYWORDS
+C
+C Elementary polynomial operations, polynomial operations.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C .. Scalar Arguments ..
+ INTEGER INFO, K
+C .. Array Arguments ..
+ DOUBLE PRECISION DWORK(*), IMZ(*), P(*), REZ(*)
+C .. Local Scalars ..
+ INTEGER I
+ DOUBLE PRECISION U, V
+C .. External Subroutines ..
+ EXTERNAL DAXPY, DCOPY, XERBLA
+C .. Executable Statements ..
+C
+C Test the input scalar arguments.
+C
+ IF( K.LT.0 ) THEN
+ INFO = -1
+C
+C Error return.
+C
+ CALL XERBLA( 'MC01PD', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ INFO = 0
+ P(1) = ONE
+ IF ( K.EQ.0 )
+ $ RETURN
+C
+ I = 1
+C WHILE ( I <= K ) DO
+ 20 IF ( I.LE.K ) THEN
+ U = REZ(I)
+ V = IMZ(I)
+ DWORK(1) = ZERO
+C
+ IF ( V.EQ.ZERO ) THEN
+ CALL DCOPY( I, P, 1, DWORK(2), 1 )
+ CALL DAXPY( I, -U, P, 1, DWORK, 1 )
+ I = I + 1
+C
+ ELSE
+ IF ( I.EQ.K ) THEN
+ INFO = K
+ RETURN
+ ELSE IF ( ( U.NE.REZ(I+1) ) .OR. ( V.NE.-IMZ(I+1) ) ) THEN
+ INFO = I + 1
+ RETURN
+ END IF
+C
+ DWORK(2) = ZERO
+ CALL DCOPY( I, P, 1, DWORK(3), 1 )
+ CALL DAXPY( I, -(U + U), P, 1, DWORK(2), 1 )
+ CALL DAXPY( I, U**2+V**2, P, 1, DWORK, 1 )
+ I = I + 2
+ END IF
+C
+ CALL DCOPY( I, DWORK, 1, P, 1 )
+ GO TO 20
+ END IF
+C END WHILE 20
+C
+ RETURN
+C *** Last line of MC01PD ***
+ END
Added: trunk/octave-forge/main/control/devel/ss2tf/MC01PY.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/MC01PY.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/MC01PY.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,157 @@
+ SUBROUTINE MC01PY( K, REZ, IMZ, P, DWORK, 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 coefficients of a real polynomial P(x) from its
+C zeros. The coefficients are stored in decreasing order of the
+C powers of x.
+C
+C ARGUMENTS
+C
+C Input/Output Parameters
+C
+C K (input) INTEGER
+C The number of zeros (and hence the degree) of P(x).
+C K >= 0.
+C
+C REZ (input) DOUBLE PRECISION array, dimension (K)
+C IMZ (input) DOUBLE PRECISION array, dimension (K)
+C The real and imaginary parts of the i-th zero of P(x)
+C must be stored in REZ(i) and IMZ(i), respectively, where
+C i = 1, 2, ..., K. The zeros may be supplied in any order,
+C except that complex conjugate zeros must appear
+C consecutively.
+C
+C P (output) DOUBLE PRECISION array, dimension (K+1)
+C This array contains the coefficients of P(x) in decreasing
+C powers of x.
+C
+C Workspace
+C
+C DWORK DOUBLE PRECISION array, dimension (K)
+C If K = 0, this array is not referenced.
+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 > 0: if INFO = i, (REZ(i),IMZ(i)) is a complex zero but
+C (REZ(i-1),IMZ(i-1)) is not its conjugate.
+C
+C METHOD
+C
+C The routine computes the coefficients of the real K-th degree
+C polynomial P(x) as
+C
+C P(x) = (x - r(1)) * (x - r(2)) * ... * (x - r(K))
+C
+C where r(i) = (REZ(i),IMZ(i)).
+C
+C Note that REZ(i) = REZ(j) and IMZ(i) = -IMZ(j) if r(i) and r(j)
+C form a complex conjugate pair (where i <> j), and that IMZ(i) = 0
+C if r(i) is real.
+C
+C NUMERICAL ASPECTS
+C
+C None.
+C
+C CONTRIBUTOR
+C
+C V. Sima, Research Institute for Informatics, Bucharest, May 2002.
+C
+C REVISIONS
+C
+C -
+C
+C KEYWORDS
+C
+C Elementary polynomial operations, polynomial operations.
+C
+C ******************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
+C .. Scalar Arguments ..
+ INTEGER INFO, K
+C .. Array Arguments ..
+ DOUBLE PRECISION DWORK(*), IMZ(*), P(*), REZ(*)
+C .. Local Scalars ..
+ INTEGER I
+ DOUBLE PRECISION U, V
+C .. External Subroutines ..
+ EXTERNAL DAXPY, DCOPY, XERBLA
+C .. Executable Statements ..
+C
+C Test the input scalar arguments.
+C
+ IF( K.LT.0 ) THEN
+ INFO = -1
+C
+C Error return.
+C
+ CALL XERBLA( 'MC01PY', -INFO )
+ RETURN
+ END IF
+C
+C Quick return if possible.
+C
+ INFO = 0
+ P(1) = ONE
+ IF ( K.EQ.0 )
+ $ RETURN
+C
+ I = 1
+C WHILE ( I <= K ) DO
+ 20 IF ( I.LE.K ) THEN
+ U = REZ(I)
+ V = IMZ(I)
+ DWORK(I) = ZERO
+C
+ IF ( V.EQ.ZERO ) THEN
+ CALL DAXPY( I, -U, P, 1, DWORK, 1 )
+C
+ ELSE
+ IF ( I.EQ.K ) THEN
+ INFO = K
+ RETURN
+ ELSE IF ( ( U.NE.REZ(I+1) ) .OR. ( V.NE.-IMZ(I+1) ) ) THEN
+ INFO = I + 1
+ RETURN
+ END IF
+C
+ DWORK(I+1) = ZERO
+ CALL DAXPY( I, -(U + U), P, 1, DWORK, 1 )
+ CALL DAXPY( I, U**2+V**2, P, 1, DWORK(2), 1 )
+ I = I + 1
+ END IF
+C
+ CALL DCOPY( I, DWORK, 1, P(2), 1 )
+ I = I + 1
+ GO TO 20
+ END IF
+C END WHILE 20
+C
+ RETURN
+C *** Last line of MC01PY ***
+ END
Added: trunk/octave-forge/main/control/devel/ss2tf/TB01ID.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/TB01ID.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/TB01ID.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,402 @@
+ SUBROUTINE TB01ID( JOB, N, M, P, MAXRED, A, LDA, B, LDB, C, LDC,
+ $ SCALE, 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 reduce the 1-norm of a system matrix
+C
+C S = ( A B )
+C ( C 0 )
+C
+C corresponding to the triple (A,B,C), by balancing. This involves
+C a diagonal similarity transformation inv(D)*A*D applied
+C iteratively to A to make the rows and columns of
+C -1
+C diag(D,I) * S * diag(D,I)
+C
+C as close in norm as possible.
+C
+C The balancing can be performed optionally on the following
+C particular system matrices
+C
+C S = A, S = ( A B ) or S = ( A )
+C ( C )
+C
+C ARGUMENTS
+C
+C Mode Parameters
+C
+C JOB CHARACTER*1
+C Indicates which matrices are involved in balancing, as
+C follows:
+C = 'A': All matrices are involved in balancing;
+C = 'B': B and A matrices are involved in balancing;
+C = 'C': C and A matrices are involved in balancing;
+C = 'N': B and C matrices are not involved in balancing.
+C
+C Input/Output Parameters
+C
+C N (input) INTEGER
+C The order of the matrix A, the number of rows of matrix B
+C and the number of columns of matrix C.
+C N represents the dimension of the state vector. N >= 0.
+C
+C M (input) INTEGER.
+C The number of columns of matrix B.
+C M represents the dimension of input vector. M >= 0.
+C
+C P (input) INTEGER.
+C The number of rows of matrix C.
+C P represents the dimension of output vector. P >= 0.
+C
+C MAXRED (input/output) DOUBLE PRECISION
+C On entry, the maximum allowed reduction in the 1-norm of
+C S (in an iteration) if zero rows or columns are
+C encountered.
+C If MAXRED > 0.0, MAXRED must be larger than one (to enable
+C the norm reduction).
+C If MAXRED <= 0.0, then the value 10.0 for MAXRED is
+C used.
+C On exit, if the 1-norm of the given matrix S is non-zero,
+C the ratio between the 1-norm of the given matrix and the
+C 1-norm of the balanced matrix.
+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 system state matrix A.
+C On exit, the leading N-by-N part of this array contains
+C the balanced matrix inv(D)*A*D.
+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, if M > 0, the leading N-by-M part of this array
+C must contain the system input matrix B.
+C On exit, if M > 0, the leading N-by-M part of this array
+C contains the balanced matrix inv(D)*B.
+C The array B is not referenced if M = 0.
+C
+C LDB INTEGER
+C The leading dimension of the array B.
+C LDB >= MAX(1,N) if M > 0.
+C LDB >= 1 if M = 0.
+C
+C C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
+C On entry, if P > 0, the leading P-by-N part of this array
+C must contain the system output matrix C.
+C On exit, if P > 0, the leading P-by-N part of this array
+C contains the balanced matrix C*D.
+C The array C is not referenced if P = 0.
+C
+C LDC INTEGER
+C The leading dimension of the array C. LDC >= MAX(1,P).
+C
+C SCALE (output) DOUBLE PRECISION array, dimension (N)
+C The scaling factors applied to S. If D(j) is the scaling
+C factor applied to row and column j, then SCALE(j) = D(j),
+C for j = 1,...,N.
+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 Balancing consists of applying a diagonal similarity
+C transformation
+C -1
+C diag(D,I) * S * diag(D,I)
+C
+C to make the 1-norms of each row of the first N rows of S and its
+C corresponding column nearly equal.
+C
+C Information about the diagonal matrix D is returned in the vector
+C SCALE.
+C
+C REFERENCES
+C
+C [1] Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J.,
+C Du Croz, J., Greenbaum, A., Hammarling, S., McKenney, A.,
+C Ostrouchov, S., and Sorensen, D.
+C LAPACK Users' Guide: Second Edition.
+C SIAM, Philadelphia, 1995.
+C
+C NUMERICAL ASPECTS
+C
+C None.
+C
+C CONTRIBUTOR
+C
+C Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Jan. 1998.
+C This subroutine is based on LAPACK routine DGEBAL, and routine
+C BALABC (A. Varga, German Aerospace Research Establishment, DLR).
+C
+C REVISIONS
+C
+C -
+C
+C KEYWORDS
+C
+C Balancing, eigenvalue, matrix algebra, matrix operations,
+C similarity transformation.
+C
+C *********************************************************************
+C
+C .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+ DOUBLE PRECISION SCLFAC
+ PARAMETER ( SCLFAC = 1.0D+1 )
+ DOUBLE PRECISION FACTOR, MAXR
+ PARAMETER ( FACTOR = 0.95D+0, MAXR = 10.0D+0 )
+C ..
+C .. Scalar Arguments ..
+ CHARACTER JOB
+ INTEGER INFO, LDA, LDB, LDC, M, N, P
+ DOUBLE PRECISION MAXRED
+C ..
+C .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
+ $ SCALE( * )
+C ..
+C .. Local Scalars ..
+ LOGICAL NOCONV, WITHB, WITHC
+ INTEGER I, ICA, IRA, J
+ DOUBLE PRECISION CA, CO, F, G, MAXNRM, RA, RO, S, SFMAX1,
+ $ SFMAX2, SFMIN1, SFMIN2, SNORM, SRED
+C ..
+C .. External Functions ..
+ LOGICAL LSAME
+ INTEGER IDAMAX
+ DOUBLE PRECISION DASUM, DLAMCH
+ EXTERNAL DASUM, DLAMCH, IDAMAX, LSAME
+C ..
+C .. External Subroutines ..
+ EXTERNAL DSCAL, XERBLA
+C ..
+C .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX, MIN
+C ..
+C .. Executable Statements ..
+C
+C Test the scalar input arguments.
+C
+ INFO = 0
+ WITHB = LSAME( JOB, 'A' ) .OR. LSAME( JOB, 'B' )
+ WITHC = LSAME( JOB, 'A' ) .OR. LSAME( JOB, 'C' )
+C
+ IF( .NOT.WITHB .AND. .NOT.WITHC .AND. .NOT.LSAME( JOB, '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( MAXRED.GT.ZERO .AND. MAXRED.LT.ONE ) THEN
+ INFO = -5
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ ELSE IF( ( M.GT.0 .AND. LDB.LT.MAX( 1, N ) ) .OR.
+ $ ( M.EQ.0 .AND. LDB.LT.1 ) ) THEN
+ INFO = -9
+ ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
+ INFO = -11
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'TB01ID', -INFO )
+ RETURN
+ END IF
+C
+ IF( N.EQ.0 )
+ $ RETURN
+C
+C Compute the 1-norm of the required part of matrix S and exit if
+C it is zero.
+C
+ SNORM = ZERO
+C
+ DO 10 J = 1, N
+ SCALE( J ) = ONE
+ CO = DASUM( N, A( 1, J ), 1 )
+ IF( WITHC .AND. P.GT.0 )
+ $ CO = CO + DASUM( P, C( 1, J ), 1 )
+ SNORM = MAX( SNORM, CO )
+ 10 CONTINUE
+C
+ IF( WITHB ) THEN
+C
+ DO 20 J = 1, M
+ SNORM = MAX( SNORM, DASUM( N, B( 1, J ), 1 ) )
+ 20 CONTINUE
+C
+ END IF
+C
+ IF( SNORM.EQ.ZERO )
+ $ RETURN
+C
+C Set some machine parameters and the maximum reduction in the
+C 1-norm of S if zero rows or columns are encountered.
+C
+ SFMIN1 = DLAMCH( 'S' ) / DLAMCH( 'P' )
+ SFMAX1 = ONE / SFMIN1
+ SFMIN2 = SFMIN1*SCLFAC
+ SFMAX2 = ONE / SFMIN2
+C
+ SRED = MAXRED
+ IF( SRED.LE.ZERO ) SRED = MAXR
+C
+ MAXNRM = MAX( SNORM/SRED, SFMIN1 )
+C
+C Balance the matrix.
+C
+C Iterative loop for norm reduction.
+C
+ 30 CONTINUE
+ NOCONV = .FALSE.
+C
+ DO 90 I = 1, N
+ CO = ZERO
+ RO = ZERO
+C
+ DO 40 J = 1, N
+ IF( J.EQ.I )
+ $ GO TO 40
+ CO = CO + ABS( A( J, I ) )
+ RO = RO + ABS( A( I, J ) )
+ 40 CONTINUE
+C
+ ICA = IDAMAX( N, A( 1, I ), 1 )
+ CA = ABS( A( ICA, I ) )
+ IRA = IDAMAX( N, A( I, 1 ), LDA )
+ RA = ABS( A( I, IRA ) )
+C
+ IF( WITHC .AND. P.GT.0 ) THEN
+ CO = CO + DASUM( P, C( 1, I ), 1 )
+ ICA = IDAMAX( P, C( 1, I ), 1 )
+ CA = MAX( CA, ABS( C( ICA, I ) ) )
+ END IF
+C
+ IF( WITHB .AND. M.GT.0 ) THEN
+ RO = RO + DASUM( M, B( I, 1 ), LDB )
+ IRA = IDAMAX( M, B( I, 1 ), LDB )
+ RA = MAX( RA, ABS( B( I, IRA ) ) )
+ END IF
+C
+C Special case of zero CO and/or RO.
+C
+ IF( CO.EQ.ZERO .AND. RO.EQ.ZERO )
+ $ GO TO 90
+ IF( CO.EQ.ZERO ) THEN
+ IF( RO.LE.MAXNRM )
+ $ GO TO 90
+ CO = MAXNRM
+ END IF
+ IF( RO.EQ.ZERO ) THEN
+ IF( CO.LE.MAXNRM )
+ $ GO TO 90
+ RO = MAXNRM
+ END IF
+C
+C Guard against zero CO or RO due to underflow.
+C
+ G = RO / SCLFAC
+ F = ONE
+ S = CO + RO
+ 50 CONTINUE
+ IF( CO.GE.G .OR. MAX( F, CO, CA ).GE.SFMAX2 .OR.
+ $ MIN( RO, G, RA ).LE.SFMIN2 )GO TO 60
+ F = F*SCLFAC
+ CO = CO*SCLFAC
+ CA = CA*SCLFAC
+ G = G / SCLFAC
+ RO = RO / SCLFAC
+ RA = RA / SCLFAC
+ GO TO 50
+C
+ 60 CONTINUE
+ G = CO / SCLFAC
+ 70 CONTINUE
+ IF( G.LT.RO .OR. MAX( RO, RA ).GE.SFMAX2 .OR.
+ $ MIN( F, CO, G, CA ).LE.SFMIN2 )GO TO 80
+ F = F / SCLFAC
+ CO = CO / SCLFAC
+ CA = CA / SCLFAC
+ G = G / SCLFAC
+ RO = RO*SCLFAC
+ RA = RA*SCLFAC
+ GO TO 70
+C
+C Now balance.
+C
+ 80 CONTINUE
+ IF( ( CO+RO ).GE.FACTOR*S )
+ $ GO TO 90
+ IF( F.LT.ONE .AND. SCALE( I ).LT.ONE ) THEN
+ IF( F*SCALE( I ).LE.SFMIN1 )
+ $ GO TO 90
+ END IF
+ IF( F.GT.ONE .AND. SCALE( I ).GT.ONE ) THEN
+ IF( SCALE( I ).GE.SFMAX1 / F )
+ $ GO TO 90
+ END IF
+ G = ONE / F
+ SCALE( I ) = SCALE( I )*F
+ NOCONV = .TRUE.
+C
+ CALL DSCAL( N, G, A( I, 1 ), LDA )
+ CALL DSCAL( N, F, A( 1, I ), 1 )
+ IF( M.GT.0 ) CALL DSCAL( M, G, B( I, 1 ), LDB )
+ IF( P.GT.0 ) CALL DSCAL( P, F, C( 1, I ), 1 )
+C
+ 90 CONTINUE
+C
+ IF( NOCONV )
+ $ GO TO 30
+C
+C Set the norm reduction parameter.
+C
+ MAXRED = SNORM
+ SNORM = ZERO
+C
+ DO 100 J = 1, N
+ CO = DASUM( N, A( 1, J ), 1 )
+ IF( WITHC .AND. P.GT.0 )
+ $ CO = CO + DASUM( P, C( 1, J ), 1 )
+ SNORM = MAX( SNORM, CO )
+ 100 CONTINUE
+C
+ IF( WITHB ) THEN
+C
+ DO 110 J = 1, M
+ SNORM = MAX( SNORM, DASUM( N, B( 1, J ), 1 ) )
+ 110 CONTINUE
+C
+ END IF
+ MAXRED = MAXRED/SNORM
+ RETURN
+C *** Last line of TB01ID ***
+ END
Added: trunk/octave-forge/main/control/devel/ss2tf/TB01ZD.f
===================================================================
--- trunk/octave-forge/main/control/devel/ss2tf/TB01ZD.f (rev 0)
+++ trunk/octave-forge/main/control/devel/ss2tf/TB01ZD.f 2010-10-08 19:53:14 UTC (rev 7817)
@@ -0,0 +1,440 @@
+ SUBROUTINE TB01ZD( JOBZ, N, P, A, LDA, B, C, LDC, NCONT, Z, LDZ,
+ $ TAU, TOL, 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 find a controllable realization for the linear time-invariant
+C single-input system
+C
+C ...
[truncated message content] |
