SUBROUTINE SB01DD( N, M, INDCON, A, LDA, B, LDB, NBLK, WR, WI,
$ Z, LDZ, Y, COUNT, G, LDG, TOL, IWORK, DWORK,
$ LDWORK, INFO )
C
C PURPOSE
C
C To compute for a controllable matrix pair ( A, B ) a matrix G
C such that the matrix A - B*G has the desired eigenstructure,
C specified by desired eigenvalues and free eigenvector elements.
C
C The pair ( A, B ) should be given in orthogonal canonical form
C as returned by the SLICOT Library routine AB01ND.
C
C ARGUMENTS
C
C Input/Output Parameters
C
C N (input) INTEGER
C The order of the matrix A and the number of rows of the
C matrix B. N >= 0.
C
C M (input) INTEGER
C The number of columns of the matrix B. M >= 0.
C
C INDCON (input) INTEGER
C The controllability index of the pair ( A, B ).
C 0 <= INDCON <= N.
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 N-by-N matrix A in orthogonal canonical form,
C as returned by SLICOT Library routine AB01ND.
C On exit, the leading N-by-N part of this array contains
C the real Schur form of the matrix A - B*G.
C The elements below the real Schur form of A are set to
C zero.
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 N-by-M matrix B in orthogonal canonical form,
C as returned by SLICOT Library routine AB01ND.
C On exit, the leading N-by-M part of this array contains
C the transformed matrix B.
C
C LDB INTEGER
C The leading dimension of the array B. LDB >= max(1,N).
C
C NBLK (input) INTEGER array, dimension (N)
C The leading INDCON elements of this array must contain the
C orders of the diagonal blocks in the orthogonal canonical
C form of A, as returned by SLICOT Library routine AB01ND.
C The values of these elements must satisfy the following
C conditions:
C NBLK(1) >= NBLK(2) >= ... >= NBLK(INDCON),
C NBLK(1) + NBLK(2) + ... + NBLK(INDCON) = N.
C
C WR (input) DOUBLE PRECISION array, dimension (N)
C WI (input) DOUBLE PRECISION array, dimension (N)
C These arrays must contain the real and imaginary parts,
C respectively, of the desired poles of the closed-loop
C system, i.e., the eigenvalues of A - B*G. The poles can be
C unordered, except that complex conjugate pairs of poles
C must appear consecutively.
C The elements of WI for complex eigenvalues are modified
C internally, but restored on exit.
C
C Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
C On entry, the leading N-by-N part of this array must
C contain the orthogonal matrix Z generated by SLICOT
C Library routine AB01ND in the reduction of ( A, B ) to
C orthogonal canonical form.
C On exit, the leading N-by-N part of this array contains
C the orthogonal transformation matrix which reduces A - B*G
C to real Schur form.
C
C LDZ INTEGER
C The leading dimension of the array Z. LDZ >= max(1,N).
C
C Y (input) DOUBLE PRECISION array, dimension (M*N)
C Y contains elements which are used as free parameters
C in the eigenstructure design. The values of these
C parameters are often set by an external optimization
C procedure.
C
C COUNT (output) INTEGER
C The actual number of elements in Y used as free
C eigenvector and feedback matrix elements in the
C eigenstructure design.
C
C G (output) DOUBLE PRECISION array, dimension (LDG,N)
C The leading M-by-N part of this array contains the
C feedback matrix which assigns the desired eigenstructure
C of A - B*G.
C
C LDG INTEGER
C The leading dimension of the array G. LDG >= max(1,M).
C
C Tolerances
C
C TOL DOUBLE PRECISION
C The tolerance to be used in rank determination when
C transforming (A, B). If the user sets TOL > 0, then
C the given value of TOL is used as a lower bound for the
C reciprocal condition number (see the description of the
C argument RCOND in the SLICOT routine MB03OD); a
C (sub)matrix whose estimated condition number is less than
C 1/TOL is considered to be of full rank. If the user sets
C TOL <= 0, then an implicitly computed, default tolerance,
C defined by TOLDEF = N*N*EPS, is used instead, where
C EPS is the machine precision (see LAPACK Library routine
C DLAMCH).
C
C Workspace
C
C IWORK INTEGER array, dimension (M)
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(M*N,M*M+2*N+4*M+1).
C For optimum 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 = 1: if the pair ( A, B ) is not controllable or the free
C parameters are not set appropriately.
C
C METHOD
C
C The routine implements the method proposed in [1], [2].
C
C REFERENCES
C
C [1] Petkov, P.Hr., Konstantinov, M.M., Gu, D.W. and
C Postlethwaite, I.
C Optimal pole assignment design of linear multi-input systems.
C Report 96-11, Department of Engineering, Leicester University,
C 1996.
C
C [2] Petkov, P.Hr., Christov, N.D. and Konstantinov, M.M.
C A computational algorithm for pole assignment of linear multi
C input systems. IEEE Trans. Automatic Control, vol. AC-31,
C pp. 1044-1047, 1986.
C
C NUMERICAL ASPECTS
C
C The method implemented is backward stable.
C
C FURTHER COMMENTS
C
C The eigenvalues of the real Schur form matrix As, returned in the
C array A, are very close to the desired eigenvalues WR+WI*i.
C However, the eigenvalues of the closed-loop matrix A - B*G,
C computed by the QR algorithm using the matrices A and B, given on
C entry, may be far from WR+WI*i, although the relative error
C norm( Z'*(A - B*G)*Z - As )/norm( As )
C is close to machine accuracy. This may happen when the eigenvalue
C problem for the matrix A - B*G is ill-conditioned.
C
C CONTRIBUTORS
C
C P.Hr. Petkov, Technical University of Sofia, Oct. 1998.
C V. Sima, Katholieke Universiteit Leuven, Jan. 1999, SLICOT Library
C version.
C
C REVISIONS
C
C V. Sima, Research Institute for Informatics, Bucharest, Mar. 2005,
C Apr. 2017.
C
C KEYWORDS
C
C Closed loop spectrum, closed loop systems, eigenvalue assignment,
C orthogonal canonical form, orthogonal transformation, pole
C placement, Schur form.
C
C ******************************************************************
C
C .. Parameters ..
DOUBLE PRECISION ZERO, ONE
PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
C
C .. Scalar Arguments ..
INTEGER COUNT, INDCON, INFO, LDA, LDB, LDG, LDWORK,
$ LDZ, M, N
DOUBLE PRECISION TOL
C ..
C .. Array Arguments ..
INTEGER IWORK( * ), NBLK( * )
DOUBLE PRECISION A( LDA, * ), B( LDB, * ), DWORK( * ),
$ G( LDG, * ), WI( * ), WR( * ), Y( * ),
$ Z( LDZ, * )
C ..
C .. Local Scalars ..
LOGICAL COMPLX
INTEGER I, IA, INDCN1, INDCN2, INDCRT, IP, IRMX, IWRK,
$ K, KK, KMR, L, LP1, M1, MAXWRK, MI, MP1, MR,
$ MR1, NBLKCR, NC, NI, NJ, NP1, NR, NR1, RANK
DOUBLE PRECISION P, Q, R, S, SVLMXA, SVLMXB, TOLDEF
C ..
C .. Local Arrays ..
DOUBLE PRECISION SVAL( 3 )
C ..
C .. External Functions ..
DOUBLE PRECISION DASUM, DLAMCH, DLANGE, DLAPY2
EXTERNAL DASUM, DLAMCH, DLANGE, DLAPY2
C ..
C .. External Subroutines ..
EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DLACPY, DLARF,
$ DLARFG, DLARTG, DLASET, DROT, DSCAL, MB02QD,
$ XERBLA
C ..
C .. Intrinsic Functions ..
INTRINSIC DBLE, INT, MAX, MIN
C ..
C .. Executable Statements ..
C
C Test the input arguments.
C
INFO = 0
NR = 0
IWRK = MAX( M*N, M*M + 2*N + 4*M + 1 )
DO 10 I = 1, MIN( INDCON, N )
NR = NR + NBLK( I )
IF( I.GT.1 ) THEN
IF( NBLK( I-1 ).LT.NBLK( I ) )
$ INFO = -8
END IF
10 CONTINUE
IF( N.LT.0 ) THEN
INFO = -1
ELSE IF( M.LT.0 ) THEN
INFO = -2
ELSE IF( INDCON.LT.0 .OR. INDCON.GT.N ) THEN
INFO = -3
ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
INFO = -5
ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
INFO = -7
ELSE IF( NR.NE.N ) THEN
INFO = -8
ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
INFO = -12
ELSE IF( LDG.LT.MAX( 1, M ) ) THEN
INFO = -16
ELSE IF( LDWORK.LT.IWRK ) THEN
INFO = -20
END IF
C
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SB01DD', -INFO )
RETURN
END IF
C
C Quick return if possible.
C
IF( MIN( M, N, INDCON ).EQ.0 ) THEN
COUNT = 0
DWORK( 1 ) = ONE
RETURN
END IF
C
MAXWRK = IWRK
TOLDEF = TOL
IF ( TOLDEF.LE.ZERO ) THEN
C
C Use the default tolerance, based on machine precision.
C
TOLDEF = DBLE( N*N )*DLAMCH( 'EPSILON' )
END IF
C
IRMX = 2*N + 1
IWRK = IRMX + M*M
M1 = NBLK( 1 )
COUNT = 1
INDCRT = INDCON
NBLKCR = NBLK( INDCRT )
C
C Compute the Frobenius norm of [ B A ] (used for rank estimation),
C taking into account the structure.
C
NR = M1
NC = 1
SVLMXB = DLANGE( 'Frobenius', M1, M, B, LDB, DWORK )
SVLMXA = ZERO
C
DO 20 I = 1, INDCRT - 1
NR = NR + NBLK( I+1 )
SVLMXA = DLAPY2( SVLMXA,
$ DLANGE( 'Frobenius', NR, NBLK( I ),
$ A( 1, NC ), LDA, DWORK ) )
NC = NC + NBLK( I )
20 CONTINUE
C
SVLMXA = DLAPY2( SVLMXA,
$ DLANGE( 'Frobenius', N, NBLKCR, A( 1, NC ), LDA,
$ DWORK ) )
L = 1
MR = NBLKCR
NR = N - MR + 1
30 CONTINUE
C WHILE( INDCRT.GT.1 )LOOP
IF( INDCRT.GT.1 ) THEN
C
C Assign next eigenvalue/eigenvector.
C
LP1 = L + M1
INDCN1 = INDCRT - 1
MR1 = NBLK( INDCN1 )
NR1 = NR - MR1
COMPLX = WI(L).NE.ZERO
CALL DCOPY( MR, Y( COUNT ), 1, DWORK( NR ), 1 )
COUNT = COUNT + MR
NC = 1
IF( COMPLX ) THEN
CALL DCOPY( MR, Y( COUNT ), 1, DWORK( N+NR ), 1 )
COUNT = COUNT + MR
WI( L+1 ) = WI( L )*WI( L+1 )
NC = 2
END IF
C
C Compute and transform eigenvector.
C
DO 50 IP = 1, INDCRT
IF( IP.NE.INDCRT ) THEN
CALL DLACPY( 'Full', MR, MR1, A( NR, NR1 ), LDA,
$ DWORK( IRMX ), M )
IF( IP.EQ.1 ) THEN
MP1 = MR
NP1 = NR + MP1
ELSE
MP1 = MR + 1
NP1 = NR + MP1
S = DASUM( MP1, DWORK( NR ), 1 )
IF( COMPLX ) S = S + DASUM( MP1, DWORK( N+NR ), 1 )
IF( S.NE.ZERO ) THEN
C
C Scale eigenvector elements.
C
CALL DSCAL( MP1, ONE/S, DWORK( NR ), 1 )
IF( COMPLX ) THEN
CALL DSCAL( MP1, ONE/S, DWORK( N+NR ), 1 )
IF( NP1.LE.N )
$ DWORK( N+NP1 ) = DWORK( N+NP1 ) / S
END IF
END IF
END IF
C
C Compute the right-hand side of the eigenvector equations.
C
CALL DCOPY( MR, DWORK( NR ), 1, DWORK( NR1 ), 1 )
CALL DSCAL( MR, WR( L ), DWORK( NR1 ), 1 )
CALL DGEMV( 'No transpose', MR, MP1, -ONE, A( NR, NR ),
$ LDA, DWORK( NR ), 1, ONE, DWORK( NR1 ), 1 )
IF( COMPLX ) THEN
CALL DAXPY( MR, WI( L+1 ), DWORK( N+NR ), 1,
$ DWORK( NR1 ), 1 )
CALL DCOPY( MR, DWORK( NR ), 1, DWORK( N+NR1 ), 1 )
CALL DAXPY( MR, WR( L+1 ), DWORK( N+NR ), 1,
$ DWORK( N+NR1 ), 1 )
CALL DGEMV( 'No transpose', MR, MP1, -ONE,
$ A( NR, NR ), LDA, DWORK( N+NR ), 1, ONE,
$ DWORK( N+NR1 ), 1 )
IF( NP1.LE.N )
$ CALL DAXPY( MR, -DWORK( N+NP1 ), A( NR, NP1 ), 1,
$ DWORK( N+NR1 ), 1 )
END IF
C
C Solve linear equations for eigenvector elements.
C
CALL MB02QD( 'FreeElements', 'NoPermuting', MR, MR1, NC,
$ TOLDEF, SVLMXA, DWORK( IRMX ), M,
$ DWORK( NR1 ), N, Y( COUNT ), IWORK, RANK,
$ SVAL, DWORK( IWRK ), LDWORK-IWRK+1, INFO )
MAXWRK = MAX( MAXWRK, INT( DWORK( IWRK ) ) + IWRK - 1 )
IF( RANK.LT.MR ) GO TO 80
C
COUNT = COUNT + ( MR1 - MR )*NC
NJ = NR1
ELSE
NJ = NR
END IF
NI = NR + MR - 1
IF( IP.EQ.1 ) THEN
KMR = MR - 1
ELSE
KMR = MR
IF( IP.EQ.2 ) THEN
NI = NI + NBLKCR
ELSE
NI = NI + NBLK( INDCRT-IP+2 ) + 1
IF( COMPLX ) NI = MIN( NI+1, N )
END IF
END IF
C
DO 40 KK = 1, KMR
K = NR + MR - KK
IF( IP.EQ.1 ) K = N - KK
CALL DLARTG( DWORK( K ), DWORK( K+1 ), P, Q, R )
DWORK( K ) = R
DWORK( K+1 ) = ZERO
C
C Transform A.
C
CALL DROT( N-NJ+1, A( K, NJ ), LDA, A( K+1, NJ ), LDA,
$ P, Q )
CALL DROT( NI, A( 1, K ), 1, A( 1, K+1 ), 1, P, Q )
C
IF( K.LT.LP1 ) THEN
C
C Transform B.
C
CALL DROT( M, B( K, 1 ), LDB, B( K+1, 1 ), LDB, P, Q )
END IF
C
C Accumulate transformations.
C
CALL DROT( N, Z( 1, K ), 1, Z( 1, K+1 ), 1, P, Q )
C
IF( COMPLX ) THEN
CALL DROT( 1, DWORK( N+K ), 1, DWORK( N+K+1 ), 1, P,
$ Q )
K = K + 1
IF( K.LT.N ) THEN
CALL DLARTG( DWORK( N+K ), DWORK( N+K+1 ), P, Q,
$ R )
DWORK( N+K ) = R
DWORK( N+K+1 ) = ZERO
C
C Transform A.
C
CALL DROT( N-NJ+1, A( K, NJ ), LDA, A( K+1, NJ ),
$ LDA, P, Q )
CALL DROT( NI, A( 1, K ), 1, A( 1, K+1 ), 1, P, Q )
C
IF( K.LE.LP1 ) THEN
C
C Transform B.
C
CALL DROT( M, B( K, 1 ), LDB, B( K+1, 1 ), LDB,
$ P, Q )
END IF
C
C Accumulate transformations.
C
CALL DROT( N, Z( 1, K ), 1, Z( 1, K+1 ), 1, P, Q )
C
END IF
END IF
40 CONTINUE
C
IF( IP.NE.INDCRT ) THEN
MR = MR1
NR = NR1
IF( IP.NE.INDCN1 ) THEN
INDCN2 = INDCRT - IP - 1
MR1 = NBLK( INDCN2 )
NR1 = NR1 - MR1
END IF
END IF
50 CONTINUE
C
IF( .NOT.COMPLX ) THEN
C
C Find one column of G.
C
CALL DLACPY( 'Full', M1, M, B( L+1, 1 ), LDB, DWORK( IRMX ),
$ M )
CALL DCOPY( M1, A( L+1, L ), 1, G( 1, L ), 1 )
ELSE
C
C Find two columns of G.
C
IF( LP1.LT.N ) THEN
LP1 = LP1 + 1
K = L + 2
ELSE
K = L + 1
END IF
CALL DLACPY( 'Full', M1, M, B( K, 1 ), LDB, DWORK( IRMX ),
$ M )
CALL DLACPY( 'Full', M1, 2, A( K, L ), LDA, G( 1, L ), LDG )
IF( K.EQ.L+1 ) THEN
G( 1, L ) = G( 1, L ) -
$ ( DWORK( N+L+1 ) / DWORK( L ) )*WI( L+1 )
G( 1, L+1 ) = G( 1, L+1 ) - WR(L+1) +
$ ( DWORK( N+L ) / DWORK( L ) )*WI( L+1 )
END IF
END IF
C
CALL MB02QD( 'FreeElements', 'NoPermuting', M1, M, NC, TOLDEF,
$ SVLMXB, DWORK( IRMX ), M, G( 1, L ), LDG,
$ Y( COUNT ), IWORK, RANK, SVAL, DWORK( IWRK ),
$ LDWORK-IWRK+1, INFO )
MAXWRK = MAX( MAXWRK, INT( DWORK( IWRK ) ) + IWRK - 1 )
IF( RANK.LT.M1 ) GO TO 80
C
COUNT = COUNT + ( M - M1 )*NC
CALL DGEMM( 'No transpose', 'No transpose', LP1, NC, M, -ONE,
$ B, LDB, G( 1, L ), LDG, ONE, A( 1, L ), LDA )
L = L + 1
NBLKCR = NBLKCR - 1
IF( NBLKCR.EQ.0 ) THEN
INDCRT = INDCRT - 1
NBLKCR = NBLK( INDCRT )
END IF
IF( COMPLX ) THEN
WI( L ) = -WI( L-1 )
L = L + 1
NBLKCR = NBLKCR - 1
IF( NBLKCR.EQ.0 ) THEN
INDCRT = INDCRT - 1
IF( INDCRT.GT.0 ) NBLKCR = NBLK( INDCRT )
END IF
END IF
MR = NBLKCR
NR = N - MR + 1
GO TO 30
END IF
C END WHILE 30
C
IF( L.LE.N ) THEN
C
C Find the remaining columns of G.
C
C QR decomposition of the free eigenvectors.
C
DO 60 I = 1, MR - 1
IA = L + I - 1
MI = MR - I + 1
CALL DCOPY( MI, Y( COUNT ), 1, DWORK( 1 ), 1 )
COUNT = COUNT + MI
CALL DLARFG( MI, DWORK( 1 ), DWORK( 2 ), 1, R )
DWORK( 1 ) = ONE
C
C Transform A.
C
CALL DLARF( 'Left', MI, MR, DWORK( 1 ), 1, R, A( IA, L ),
$ LDA, DWORK( N+1 ) )
CALL DLARF( 'Right', N, MI, DWORK( 1 ), 1, R, A( 1, IA ),
$ LDA, DWORK( N+1 ) )
C
C Transform B.
C
CALL DLARF( 'Left', MI, M, DWORK( 1 ), 1, R, B( IA, 1 ),
$ LDB, DWORK( N+1 ) )
C
C Accumulate transformations.
C
CALL DLARF( 'Right', N, MI, DWORK( 1 ), 1, R, Z( 1, IA ),
$ LDZ, DWORK( N+1 ) )
60 CONTINUE
C
I = 0
C REPEAT
70 CONTINUE
I = I + 1
IA = L + I - 1
IF( WI( IA ).EQ.ZERO ) THEN
CALL DCOPY( MR, A( IA, L ), LDA, G( I, L ), LDG )
CALL DAXPY( MR-I, -ONE, Y( COUNT ), 1, G( I, L+I ), LDG )
COUNT = COUNT + MR - I
G( I, IA ) = G( I, IA ) - WR( IA )
ELSE
CALL DLACPY( 'Full', 2, MR, A( IA, L ), LDA, G( I, L ),
$ LDG )
CALL DAXPY( MR-I-1, -ONE, Y( COUNT ), 2, G( I, L+I+1 ),
$ LDG )
CALL DAXPY( MR-I-1, -ONE, Y( COUNT+1 ), 2,
$ G( I+1, L+I+1 ), LDG )
COUNT = COUNT + 2*( MR - I - 1 )
G( I, IA ) = G(I, IA ) - WR( IA )
G( I, IA+1 ) = G(I, IA+1 ) - WI( IA )
G( I+1, IA ) = G(I+1, IA ) - WI( IA+1 )
G( I+1, IA+1 ) = G(I+1, IA+1 ) - WR( IA+1 )
I = I + 1
END IF
IF( I.LT.MR ) GO TO 70
C UNTIL I.GE.MR
C
CALL DLACPY( 'Full', MR, M, B( L, 1 ), LDB, DWORK( IRMX ), M )
CALL MB02QD( 'FreeElements', 'NoPermuting', MR, M, MR, TOLDEF,
$ SVLMXB, DWORK( IRMX ), M, G( 1, L ), LDG,
$ Y( COUNT ), IWORK, RANK, SVAL, DWORK( IWRK ),
$ LDWORK-IWRK+1, INFO )
MAXWRK = MAX( MAXWRK, INT( DWORK( IWRK ) ) + IWRK - 1 )
IF( RANK.LT.MR ) GO TO 80
C
COUNT = COUNT + ( M - MR )*MR
CALL DGEMM( 'No transpose', 'No transpose', N, MR, M, -ONE, B,
$ LDB, G( 1, L ), LDG, ONE, A( 1, L ), LDA )
END IF
C
C Transform G:
C G := G * Z'.
C
CALL DGEMM( 'No transpose', 'Transpose', M, N, N, ONE, G, LDG,
$ Z, LDZ, ZERO, DWORK( 1 ), M )
CALL DLACPY( 'Full', M, N, DWORK( 1 ), M, G, LDG )
COUNT = COUNT - 1
C
IF( N.GT.2) THEN
C
C Set the elements of A below the Hessenberg part to zero.
C
CALL DLASET( 'Lower', N-2, N-2, ZERO, ZERO, A( 3, 1 ), LDA )
END IF
DWORK( 1 ) = MAXWRK
RETURN
C
C Exit with INFO = 1 if the pair ( A, B ) is not controllable or
C the free parameters are not set appropriately.
C
80 INFO = 1
RETURN
C *** Last line of SB01DD ***
END