Purpose
To construct and solve a linear algebraic system of order 2*M whose coefficient matrix has zeros below the third subdiagonal, and zero elements on the third subdiagonal with even column indices. Such systems appear when solving discrete-time Sylvester equations using the Hessenberg-Schur method.Specification
SUBROUTINE SB04QU( N, M, IND, A, LDA, B, LDB, C, LDC, D, IPR, $ INFO ) C .. Scalar Arguments .. INTEGER INFO, IND, LDA, LDB, LDC, M, N C .. Array Arguments .. INTEGER IPR(*) DOUBLE PRECISION A(LDA,*), B(LDB,*), C(LDC,*), D(*)Arguments
Input/Output Parameters
N (input) INTEGER The order of the matrix B. N >= 0. M (input) INTEGER The order of the matrix A. M >= 0. IND (input) INTEGER IND and IND - 1 specify the indices of the columns in C to be computed. IND > 1. A (input) DOUBLE PRECISION array, dimension (LDA,M) The leading M-by-M part of this array must contain an upper Hessenberg matrix. LDA INTEGER The leading dimension of array A. LDA >= MAX(1,M). B (input) DOUBLE PRECISION array, dimension (LDB,N) The leading N-by-N part of this array must contain a matrix in real Schur form. LDB INTEGER The leading dimension of array B. LDB >= MAX(1,N). C (input/output) DOUBLE PRECISION array, dimension (LDC,N) On entry, the leading M-by-N part of this array must contain the coefficient matrix C of the equation. On exit, the leading M-by-N part of this array contains the matrix C with columns IND-1 and IND updated. LDC INTEGER The leading dimension of array C. LDC >= MAX(1,M).Workspace
D DOUBLE PRECISION array, dimension (2*M*M+8*M) IPR INTEGER array, dimension (4*M)Error Indicator
INFO INTEGER = 0: successful exit; > 0: if INFO = IND, a singular matrix was encountered.Method
A special linear algebraic system of order 2*M, whose coefficient matrix has zeros below the third subdiagonal and zero elements on the third subdiagonal with even column indices, is constructed and solved. The coefficient matrix is stored compactly, row-wise.References
[1] Golub, G.H., Nash, S. and Van Loan, C.F. A Hessenberg-Schur method for the problem AX + XB = C. IEEE Trans. Auto. Contr., AC-24, pp. 909-913, 1979. [2] Sima, V. Algorithms for Linear-quadratic Optimization. Marcel Dekker, Inc., New York, 1996.Numerical Aspects
None.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None