Purpose
To annihilate one or two entries on the subdiagonal of the Hessenberg matrix A for dealing with zero elements on the diagonal of the triangular matrix B. MB03YA is an auxiliary routine called by SLICOT Library routines MB03XP and MB03YD.Specification
SUBROUTINE MB03YA( WANTT, WANTQ, WANTZ, N, ILO, IHI, ILOQ, IHIQ, $ POS, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO ) C .. Scalar Arguments .. LOGICAL WANTQ, WANTT, WANTZ INTEGER IHI, IHIQ, ILO, ILOQ, INFO, LDA, LDB, LDQ, LDZ, $ N, POS C .. Array Arguments .. DOUBLE PRECISION A(LDA,*), B(LDB,*), Q(LDQ,*), Z(LDZ,*)Arguments
Mode Parameters
WANTT LOGICAL Indicates whether the user wishes to compute the full Schur form or the eigenvalues only, as follows: = .TRUE. : Compute the full Schur form; = .FALSE.: compute the eigenvalues only. WANTQ LOGICAL Indicates whether or not the user wishes to accumulate the matrix Q as follows: = .TRUE. : The matrix Q is updated; = .FALSE.: the matrix Q is not required. WANTZ LOGICAL Indicates whether or not the user wishes to accumulate the matrix Z as follows: = .TRUE. : The matrix Z is updated; = .FALSE.: the matrix Z is not required.Input/Output Parameters
N (input) INTEGER The order of the matrices A and B. N >= 0. ILO (input) INTEGER IHI (input) INTEGER It is assumed that the matrices A and B are already (quasi) upper triangular in rows and columns 1:ILO-1 and IHI+1:N. The routine works primarily with the submatrices in rows and columns ILO to IHI, but applies the transformations to all the rows and columns of the matrices A and B, if WANTT = .TRUE.. 1 <= ILO <= max(1,N); min(ILO,N) <= IHI <= N. ILOQ (input) INTEGER IHIQ (input) INTEGER Specify the rows of Q and Z to which transformations must be applied if WANTQ = .TRUE. and WANTZ = .TRUE., respectively. 1 <= ILOQ <= ILO; IHI <= IHIQ <= N. POS (input) INTEGER The position of the zero element on the diagonal of B. ILO <= POS <= IHI. A (input/output) DOUBLE PRECISION array, dimension (LDA,N) On entry, the leading N-by-N part of this array must contain the upper Hessenberg matrix A. On exit, the leading N-by-N part of this array contains the updated matrix A where A(POS,POS-1) = 0, if POS > ILO, and A(POS+1,POS) = 0, if POS < IHI. LDA INTEGER The leading dimension of the array A. LDA >= MAX(1,N). B (input/output) DOUBLE PRECISION array, dimension (LDB,N) On entry, the leading N-by-N part of this array must contain an upper triangular matrix B with B(POS,POS) = 0. On exit, the leading N-by-N part of this array contains the updated upper triangular matrix B. LDB INTEGER The leading dimension of the array B. LDB >= MAX(1,N). Q (input/output) DOUBLE PRECISION array, dimension (LDQ,N) On entry, if WANTQ = .TRUE., then the leading N-by-N part of this array must contain the current matrix Q of transformations accumulated by MB03XP. On exit, if WANTQ = .TRUE., then the leading N-by-N part of this array contains the matrix Q updated in the submatrix Q(ILOQ:IHIQ,ILO:IHI). If WANTQ = .FALSE., Q is not referenced. LDQ INTEGER The leading dimension of the array Q. LDQ >= 1. If WANTQ = .TRUE., LDQ >= MAX(1,N). Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) On entry, if WANTZ = .TRUE., then the leading N-by-N part of this array must contain the current matrix Z of transformations accumulated by MB03XP. On exit, if WANTZ = .TRUE., then the leading N-by-N part of this array contains the matrix Z updated in the submatrix Z(ILOQ:IHIQ,ILO:IHI). If WANTZ = .FALSE., Z is not referenced. LDZ INTEGER The leading dimension of the array Z. LDZ >= 1. If WANTZ = .TRUE., LDZ >= MAX(1,N).Error Indicator
INFO INTEGER = 0: successful exit; < 0: if INFO = -i, the i-th argument had an illegal value.Method
The method is illustrated by Wilkinson diagrams for N = 5, POS = 3: [ x x x x x ] [ x x x x x ] [ x x x x x ] [ o x x x x ] A = [ o x x x x ], B = [ o o o x x ]. [ o o x x x ] [ o o o x x ] [ o o o x x ] [ o o o o x ] First, a QR factorization is applied to A(1:3,1:3) and the resulting nonzero in the updated matrix B is immediately annihilated by a Givens rotation acting on columns 1 and 2: [ x x x x x ] [ x x x x x ] [ x x x x x ] [ o x x x x ] A = [ o o x x x ], B = [ o o o x x ]. [ o o x x x ] [ o o o x x ] [ o o o x x ] [ o o o o x ] Secondly, an RQ factorization is applied to A(4:5,4:5) and the resulting nonzero in the updated matrix B is immediately annihilated by a Givens rotation acting on rows 4 and 5: [ x x x x x ] [ x x x x x ] [ x x x x x ] [ o x x x x ] A = [ o o x x x ], B = [ o o o x x ]. [ o o o x x ] [ o o o x x ] [ o o o x x ] [ o o o o x ]References
[1] Bojanczyk, A.W., Golub, G.H., and Van Dooren, P. The periodic Schur decomposition: Algorithms and applications. Proc. of the SPIE Conference (F.T. Luk, Ed.), 1770, pp. 31-42, 1992.Numerical Aspects
The algorithm requires O(N**2) floating point operations and is backward stable.Further Comments
NoneExample
Program Text
NoneProgram Data
NoneProgram Results
None