stepan
/
dynare
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
dynare/mex/sources/libslicot/SB03SY.f

452 lines
14 KiB
Fortran
Raw Blame History

SUBROUTINE SB03SY( JOB, TRANA, LYAPUN, N, T, LDT, U, LDU, XA,
$ LDXA, SEPD, THNORM, 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 estimate the "separation" between the matrices op(A) and
C op(A)',
C
C sepd(op(A),op(A)') = min norm(op(A)'*X*op(A) - X)/norm(X)
C = 1 / norm(inv(Omega))
C
C and/or the 1-norm of Theta, where op(A) = A or A' (A**T), and
C Omega and Theta are linear operators associated to the real
C discrete-time Lyapunov matrix equation
C
C op(A)'*X*op(A) - X = C,
C
C defined by
C
C Omega(W) = op(A)'*W*op(A) - W,
C Theta(W) = inv(Omega(op(W)'*X*op(A) + op(A)'*X*op(W))).
C
C The 1-norm condition estimators are used.
C
C ARGUMENTS
C
C Mode Parameters
C
C JOB CHARACTER*1
C Specifies the computation to be performed, as follows:
C = 'S': Compute the separation only;
C = 'T': Compute the norm of Theta only;
C = 'B': Compute both the separation and the norm of Theta.
C
C TRANA CHARACTER*1
C Specifies the form of op(A) to be used, as follows:
C = 'N': op(A) = A (No transpose);
C = 'T': op(A) = A**T (Transpose);
C = 'C': op(A) = A**T (Conjugate transpose = Transpose).
C
C LYAPUN CHARACTER*1
C Specifies whether or not the original Lyapunov equations
C should be solved, as follows:
C = 'O': Solve the original Lyapunov equations, updating
C the right-hand sides and solutions with the
C matrix U, e.g., X <-- U'*X*U;
C = 'R': Solve reduced Lyapunov equations only, without
C updating the right-hand sides and solutions.
C
C Input/Output Parameters
C
C N (input) INTEGER
C The order of the matrices A and X. N >= 0.
C
C T (input) DOUBLE PRECISION array, dimension (LDT,N)
C The leading N-by-N upper Hessenberg part of this array
C must contain the upper quasi-triangular matrix T in Schur
C canonical form from a Schur factorization of A.
C
C LDT INTEGER
C The leading dimension of array T. LDT >= MAX(1,N).
C
C U (input) DOUBLE PRECISION array, dimension (LDU,N)
C The leading N-by-N part of this array must contain the
C orthogonal matrix U from a real Schur factorization of A.
C If LYAPUN = 'R', the array U is not referenced.
C
C LDU INTEGER
C The leading dimension of array U.
C LDU >= 1, if LYAPUN = 'R';
C LDU >= MAX(1,N), if LYAPUN = 'O'.
C
C XA (input) DOUBLE PRECISION array, dimension (LDXA,N)
C The leading N-by-N part of this array must contain the
C matrix product X*op(A), if LYAPUN = 'O', or U'*X*U*op(T),
C if LYAPUN = 'R', in the Lyapunov equation.
C If JOB = 'S', the array XA is not referenced.
C
C LDXA INTEGER
C The leading dimension of array XA.
C LDXA >= 1, if JOB = 'S';
C LDXA >= MAX(1,N), if JOB = 'T' or 'B'.
C
C SEPD (output) DOUBLE PRECISION
C If JOB = 'S' or JOB = 'B', and INFO >= 0, SEPD contains
C the estimated quantity sepd(op(A),op(A)').
C If JOB = 'T' or N = 0, SEPD is not referenced.
C
C THNORM (output) DOUBLE PRECISION
C If JOB = 'T' or JOB = 'B', and INFO >= 0, THNORM contains
C the estimated 1-norm of operator Theta.
C If JOB = 'S' or N = 0, THNORM is not referenced.
C
C Workspace
C
C IWORK INTEGER array, dimension (N*N)
C
C DWORK DOUBLE PRECISION array, dimension (LDWORK)
C
C LDWORK INTEGER
C The length of the array DWORK.
C LDWORK >= 0, if N = 0;
C LDWORK >= MAX(3,2*N*N), if N > 0.
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 = N+1: if T has (almost) reciprocal eigenvalues;
C perturbed values were used to solve Lyapunov
C equations (but the matrix T is unchanged).
C
C METHOD
C
C SEPD is defined as
C
C sepd( op(A), op(A)' ) = sigma_min( K )
C
C where sigma_min(K) is the smallest singular value of the
C N*N-by-N*N matrix
C
C K = kprod( op(A)', op(A)' ) - I(N**2).
C
C I(N**2) is an N*N-by-N*N identity matrix, and kprod denotes the
C Kronecker product. The routine estimates sigma_min(K) by the
C reciprocal of an estimate of the 1-norm of inverse(K), computed as
C suggested in [1]. This involves the solution of several discrete-
C time Lyapunov equations, either direct or transposed. The true
C reciprocal 1-norm of inverse(K) cannot differ from sigma_min(K) by
C more than a factor of N.
C The 1-norm of Theta is estimated similarly.
C
C REFERENCES
C
C [1] Higham, N.J.
C FORTRAN codes for estimating the one-norm of a real or
C complex matrix, with applications to condition estimation.
C ACM Trans. Math. Softw., 14, pp. 381-396, 1988.
C
C NUMERICAL ASPECTS
C 3
C The algorithm requires 0(N ) operations.
C
C FURTHER COMMENTS
C
C When SEPD is zero, the routine returns immediately, with THNORM
C (if requested) not set. In this case, the equation is singular.
C The option LYAPUN = 'R' may occasionally produce slightly worse
C or better estimates, and it is much faster than the option 'O'.
C
C CONTRIBUTOR
C
C V. Sima, Research Institute for Informatics, Bucharest, Romania,
C Oct. 1998. Partly based on DDLSVX (and then SB03SD) by P. Petkov,
C Tech. University of Sofia, March 1998 (and December 1998).
C
C REVISIONS
C
C February 6, 1999, V. Sima, Katholieke Univ. Leuven, Belgium.
C V. Sima, Research Institute for Informatics, Bucharest, Oct. 2004.
C
C KEYWORDS
C
C Lyapunov equation, orthogonal transformation, real Schur form.
C
C ******************************************************************
C
C .. Parameters ..
DOUBLE PRECISION ZERO, ONE, HALF
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 )
C ..
C .. Scalar Arguments ..
CHARACTER JOB, LYAPUN, TRANA
INTEGER INFO, LDT, LDU, LDWORK, LDXA, N
DOUBLE PRECISION SEPD, THNORM
C ..
C .. Array Arguments ..
INTEGER IWORK( * )
DOUBLE PRECISION DWORK( * ), T( LDT, * ), U( LDU, * ),
$ XA( LDXA, * )
C ..
C .. Local Scalars ..
LOGICAL NOTRNA, UPDATE, WANTS, WANTT
CHARACTER TRANAT, UPLO
INTEGER INFO2, ITMP, KASE, NN
DOUBLE PRECISION BIGNUM, EST, SCALE
C ..
C .. External Functions ..
LOGICAL LSAME
DOUBLE PRECISION DLAMCH, DLANSY
EXTERNAL DLAMCH, DLANSY, LSAME
C ..
C .. External Subroutines ..
EXTERNAL DLACON, DLACPY, DSCAL, DSYR2K, MA02ED, MB01RU,
$ SB03MX, XERBLA
C ..
C .. Intrinsic Functions ..
INTRINSIC MAX
C ..
C .. Executable Statements ..
C
C Decode and Test input parameters.
C
WANTS = LSAME( JOB, 'S' )
WANTT = LSAME( JOB, 'T' )
NOTRNA = LSAME( TRANA, 'N' )
UPDATE = LSAME( LYAPUN, 'O' )
C
NN = N*N
INFO = 0
IF( .NOT. ( WANTS .OR. WANTT .OR. LSAME( JOB, 'B' ) ) ) THEN
INFO = -1
ELSE IF( .NOT.( NOTRNA .OR. LSAME( TRANA, 'T' ) .OR.
$ LSAME( TRANA, 'C' ) ) ) THEN
INFO = -2
ELSE IF( .NOT.( UPDATE .OR. LSAME( LYAPUN, 'R' ) ) ) THEN
INFO = -3
ELSE IF( N.LT.0 ) THEN
INFO = -4
ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
INFO = -6
ELSE IF( LDU.LT.1 .OR. ( UPDATE .AND. LDU.LT.N ) ) THEN
INFO = -8
ELSE IF( LDXA.LT.1 .OR. ( .NOT.WANTS .AND. LDXA.LT.N ) ) THEN
INFO = -10
ELSE IF( LDWORK.LT.0 .OR.
$ ( LDWORK.LT.MAX( 3, 2*NN ) .AND. N.GT.0 ) ) THEN
INFO = -15
END IF
C
IF( INFO.NE.0 ) THEN
CALL XERBLA( 'SB03SY', -INFO )
RETURN
END IF
C
C Quick return if possible.
C
IF( N.EQ.0 )
$ RETURN
C
ITMP = NN + 1
C
IF( NOTRNA ) THEN
TRANAT = 'T'
ELSE
TRANAT = 'N'
END IF
C
IF( .NOT.WANTT ) THEN
C
C Estimate sepd(op(A),op(A)').
C Workspace: max(3,2*N*N).
C
KASE = 0
C
C REPEAT
10 CONTINUE
CALL DLACON( NN, DWORK( ITMP ), DWORK, IWORK, EST, KASE )
IF( KASE.NE.0 ) THEN
C
C Select the triangular part of symmetric matrix to be used.
C
IF( DLANSY( '1-norm', 'Upper', N, DWORK, N, DWORK( ITMP ) )
$ .GE.
$ DLANSY( '1-norm', 'Lower', N, DWORK, N, DWORK( ITMP ) )
$ ) THEN
UPLO = 'U'
ELSE
UPLO = 'L'
END IF
C
IF( UPDATE ) THEN
C
C Transform the right-hand side: RHS := U'*RHS*U.
C
CALL MB01RU( UPLO, 'Transpose', N, N, ZERO, ONE, DWORK,
$ N, U, LDU, DWORK, N, DWORK( ITMP ), NN,
$ INFO2 )
CALL DSCAL( N, HALF, DWORK, N+1 )
END IF
CALL MA02ED( UPLO, N, DWORK, N )
C
IF( KASE.EQ.1 ) THEN
C
C Solve op(T)'*Y*op(T) - Y = scale*RHS.
C
CALL SB03MX( TRANA, N, T, LDT, DWORK, N, SCALE,
$ DWORK( ITMP ), INFO2 )
ELSE
C
C Solve op(T)*W*op(T)' - W = scale*RHS.
C
CALL SB03MX( TRANAT, N, T, LDT, DWORK, N, SCALE,
$ DWORK( ITMP ), INFO2 )
END IF
C
IF( INFO2.GT.0 )
$ INFO = N + 1
C
IF( UPDATE ) THEN
C
C Transform back to obtain the solution: Z := U*Z*U', with
C Z = Y or Z = W.
C
CALL MB01RU( UPLO, 'No transpose', N, N, ZERO, ONE,
$ DWORK, N, U, LDU, DWORK, N, DWORK( ITMP ),
$ NN, INFO2 )
CALL DSCAL( N, HALF, DWORK, N+1 )
C
C Fill in the remaining triangle of the symmetric matrix.
C
CALL MA02ED( UPLO, N, DWORK, N )
END IF
C
GO TO 10
END IF
C UNTIL KASE = 0
C
IF( EST.GT.SCALE ) THEN
SEPD = SCALE / EST
ELSE
BIGNUM = ONE / DLAMCH( 'Safe minimum' )
IF( SCALE.LT.EST*BIGNUM ) THEN
SEPD = SCALE / EST
ELSE
SEPD = BIGNUM
END IF
END IF
C
C Return if the equation is singular.
C
IF( SEPD.EQ.ZERO )
$ RETURN
END IF
C
IF( .NOT.WANTS ) THEN
C
C Estimate norm(Theta).
C Workspace: max(3,2*N*N).
C
KASE = 0
C
C REPEAT
20 CONTINUE
CALL DLACON( NN, DWORK( ITMP ), DWORK, IWORK, EST, KASE )
IF( KASE.NE.0 ) THEN
C
C Select the triangular part of symmetric matrix to be used.
C
IF( DLANSY( '1-norm', 'Upper', N, DWORK, N, DWORK( ITMP ) )
$ .GE.
$ DLANSY( '1-norm', 'Lower', N, DWORK, N, DWORK( ITMP ) )
$ ) THEN
UPLO = 'U'
ELSE
UPLO = 'L'
END IF
C
C Fill in the remaining triangle of the symmetric matrix.
C
CALL MA02ED( UPLO, N, DWORK, N )
C
C Compute RHS = op(W)'*X*op(A) + op(A)'*X*op(W).
C
CALL DSYR2K( UPLO, TRANAT, N, N, ONE, DWORK, N, XA, LDXA,
$ ZERO, DWORK( ITMP ), N )
CALL DLACPY( UPLO, N, N, DWORK( ITMP ), N, DWORK, N )
C
IF( UPDATE ) THEN
C
C Transform the right-hand side: RHS := U'*RHS*U.
C
CALL MB01RU( UPLO, 'Transpose', N, N, ZERO, ONE, DWORK,
$ N, U, LDU, DWORK, N, DWORK( ITMP ), NN,
$ INFO2 )
CALL DSCAL( N, HALF, DWORK, N+1 )
END IF
CALL MA02ED( UPLO, N, DWORK, N )
C
IF( KASE.EQ.1 ) THEN
C
C Solve op(T)'*Y*op(T) - Y = scale*RHS.
C
CALL SB03MX( TRANA, N, T, LDT, DWORK, N, SCALE,
$ DWORK( ITMP ), INFO2 )
ELSE
C
C Solve op(T)*W*op(T)' - W = scale*RHS.
C
CALL SB03MX( TRANAT, N, T, LDT, DWORK, N, SCALE,
$ DWORK( ITMP ), INFO2 )
END IF
C
IF( INFO2.GT.0 )
$ INFO = N + 1
C
IF( UPDATE ) THEN
C
C Transform back to obtain the solution: Z := U*Z*U', with
C Z = Y or Z = W.
C
CALL MB01RU( UPLO, 'No transpose', N, N, ZERO, ONE,
$ DWORK, N, U, LDU, DWORK, N, DWORK( ITMP ),
$ NN, INFO2 )
CALL DSCAL( N, HALF, DWORK, N+1 )
C
C Fill in the remaining triangle of the symmetric matrix.
C
CALL MA02ED( UPLO, N, DWORK, N )
END IF
C
GO TO 20
END IF
C UNTIL KASE = 0
C
IF( EST.LT.SCALE ) THEN
THNORM = EST / SCALE
ELSE
BIGNUM = ONE / DLAMCH( 'Safe minimum' )
IF( EST.LT.SCALE*BIGNUM ) THEN
THNORM = EST / SCALE
ELSE
THNORM = BIGNUM
END IF
END IF
END IF
C
RETURN
C *** Last line of SB03SY ***
END
Reference in New Issue View Git Blame Copy Permalink