Technische Universit¨at Berlin
Institut f¨ur Mathematik
STCSSP: A FORTRAN 77 routine to
compute a structured staircase form for a
(skew-)symmetric/(skew-)symmetric
matrix pencil.
T. Br¨ull and V. Mehrmann
Preprint 31-2007
Preprint-Reihe des Instituts f¨ur Mathematik
Technische Universit¨at Berlin
Report 31-2007 March 2008
STCSSP: A FORTRAN 77 routine to compute a structured
staircase form for a (skew-) symmetric / (skew-) symmetric
matrix pencil. ∗
Tobias Br¨ull †Volker Mehrmann †
March 19, 2008
1 Purpose
Let (N, H)∈Rn,n ×Rn,n be a matrix pencil that satisfies one of the conditions
N=NTand H=HT, (1a)
N=NTand H=−HT, (1b)
N=−NTand H=HT, (2a)
N=−NTand H=−HT. (2b)
Then we call (N, H)(skew-) symmetric / (skew-) symmetric. Likewise, one could say, that a
pencil (N, H) is (skew-) symmetric / (skew-) symmetric if and only if
N=opN(N) and H=opH(H),(3)
where opN(N) = NTor opN(N) = −NTand opH(H) = HTor opH(H) = −HT. Pencils
that satisfy (2a) are called even, see [1]. STCSSP computes a structured staircase form for a
real (skew-) symmetric / (skew-) symmetric matrix pencil. The staircase form is achieved by
an orthogonal transformation of the form
(UTNU, UTHU) = (Nnew, Hnew),(4)
with U∈Rn,n orthogonal. Obviously, we have (Nnew, Hnew) = (opN(Nnew), opH(Hnew))
and thus the pencil in staircase form (Nnew, Hnew) is (skew-) symmetric / (skew-) symmetric
again.
For a real symmetric matrix Awe know that all eigenvalues are real. Thus, one can count
the positive, negative, and zero eigenvalues. We call the triple (π, ν, ξ) the inertia index of A,
if Ahas πpositive eigenvalues, νnegative eigenvalues, and ξzero eigenvalues.
∗The work has been supported by DFG research grant no. Me790/16-1.
†Institut f¨ur Mathematik, MA 4-5, Technische Universit¨at Berlin, Straße des 17. Juni 136, D-10623 Berlin,
Fed. Rep. Germany. E-mail: {bruell,mehrmann}@math.tu-berlin.de.
1
2 The theory
The algorithm implemented is based on the following theorem, which is a slight generalization
of [1, Theorem 3.1] from even matrix pencils to (skew-) symmetric / (skew-) symmetric matrix
pencils.
Theorem 2.1. (Skew-) symmetric / (skew-) symmetric staircase form. With the
operators defined in (3) consider the (skew-) symmetric / (skew-) symmetric matrix pencil
(N, H) = (opN(N), opH(H)), where N, H ∈Rn,n. There exists a real orthogonal matrix
U∈Rn,n, such that
UTNU =
n1
.
.
.
.
.
.
nm
l
qm
.
.
.
q2
q1
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
N11 . . . . . . N1,m N1,m+1 N1,m+2 . . . N1,2m0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.Nm−1,m+2 .
.
.
opN(N1,m)· · · · · · Nm,m Nm,m+1 0
opN(N1,m+1). . . . . . opN(Nm,m+1)Nm+1,m+1
opN(N1,m+2)· · · opN(Nm−1,m+2) 0
.
.
.
.
.
.
.
.
.
opN(N1,2m).
.
.
0
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
(5)
UTHU =
n1
.
.
.
.
.
.
nm
l
qm
.
.
.
.
.
.
q1
2
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
6
4
H11 · · · · · · H1,m H1,m+1 H1,m+2 . . . . . . H1,2m+1
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
opH(H1,m). . . . . . Hm,m Hm,m+1 Hm,m+2
opH(H1,m+1). . . . . . opH(Hm,m+1)Hm+1,m+1
opH(H1,m+2). . . . . . opH(Hm,m+2)
.
.
.
.
.
.
.
.
.
.
.
.
opH(H1,2m+1)
3
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
7
5
,
where for i= 1,...,m we have Nii =opN(Nii),Hii =opH(Hii). Further, we know that
q1≥n1≥q2≥n2≥...≥qm≥nm,
Nj,2m+1−j∈Rnj,qj+1 ,1≤j≤m−1,
Nm+1,m+1 =∆ 0
0 0 ,∆ = opN(∆) ∈Rp,p,
Hj,2m+2−j=Γj0∈Rnj,qj,Γj∈Rnj,nj,1≤j≤m,
Hm+1,m+1 =Σ11 Σ12
Σ21 Σ22 ,Σ11 ∈Rp,p,Σ22 ∈Rl−p,l−p,
Hm+1,m+1 =opH(Hm+1,m+1),
and the blocks Σ22 and ∆and Γj,j= 1,...,m are nonsingular.
Note, that what is called pin Theorem 2.1 is corresponding to 2pin [1, Theorem 3.1]. The
form (5) is by far not unique. Even the quantities q1, n1, q2, n2,... are not unique, but they
become unique when using the following algorithm to compute the form (5) which is a slight
2
generalization of [1, Algorithm 1] from even matrix pencils to (skew-) symmetric / (skew-)
symmetric matrix pencils and represents a constructive proof of Theorem 2.1.
Algorithm 2.2. Staircase algorithm for (skew-) symmetric / (skew-) symmetric matrix pen-
cils.
With the operators defined in (3) consider the (skew-) symmetric / (skew-) symmetric matrix
pencil (N, H) = (opN(N), opH(H)), where N, H ∈Rn,n. Then this algorithm computes an
orthogonal matrix U∈Rn,n such that UTNU,UTHU are in the form (5). In addition, for
each of the matrices Nand H, which is real symmetric, the algorithm produces a unique
sequence of inertia indices, i.e., if only Nor His real symmetric one sequence of inertia
indices is generated and if both Nand Hare real symmetric two sequences of inertia indices
are generated.
Set flag = 0,m=n0=q0=r0= 0,l=n,
N=N22 =N, H=H, U =I.
DO WHILE flag = 0
Perform a rank revealing factorization of N22 ∈Rl−rm,l−rm,
N22 =U1∆ 0
0 0 UT
1,
with ∆ = opN(∆) ∈Rp,p nonsingular. If the matrix Nis real symmetric, also
store the inertia indices of ∆as (πN
m+1, νN
m+1,0). Set
N1=U10
0IrmT
NU10
0Irm=∆ 0
0 0 ,
H1=U10
0IrmT
HU10
0Irm=ˆ
H11 ˆ
H12
opH(ˆ
H12)ˆ
H22 ,
partitioned analogously, with ˆ
H11 =opH(ˆ
H11), and ˆ
H22 =opH(ˆ
H22).
(Here ˆ
H22 ∈Rl−p,l−p).
IF p=lTHEN
Set flag = 1 and
U=
In1+...+nm0 0
0U10
0 0 Iq1+...+qm
.(6)
ELSE
Set m=m+ 1.
Perform a rank revealing decomposition of ˆ
H22,
ˆ
H22 =U2Σ 0
0 0 UT
2,
where Σ = opH(Σ) ∈Rµ,µ is nonsingular. If the matrix His real symmetric,
also store the inertia indices of Σas (πm, νm,0). Set rm=µ=πm+νm.
3
Set
N2=Ip0
0U2T
N1Ip0
0U2=
∆ 0 0
0 0 0
0 0 0
,
H2=Ip0
0U2T
H1Ip0
0U2=
˜
H11 ˜
H12 ˜
H13
opH(˜
H12) Σ 0
opH(˜
H13) 0 0
,
partitioned analogously.
IF µ=l−pTHEN
Set flag = 1 and
ˆ
U=U10
0Irm−1 Ip0
0U2,
U=
In1+...+nm−10 0
0ˆ
U0
0 0 Iq1+...+qm−1
ELSE
Perform a rank revealing factorization or SVD
˜
H13 =U3Γm0
0 0 VT
3,
where Γm∈Rτ,τ is nonsingular.
Set nm=τ,qm=l−p−µand
N3=
U30 0
0Iµ0
0 0 V3
T
N2
U30 0
0Iµ0
0 0 V3
=
N11 N12 0 0 0
opN(N12)N22 0 0 0
00 0 0 0
0 0 0 0 0
00 0 0 0
,
H3=
U30 0
0Iµ0
0 0 V3
T
H2
U30 0
0Iµ0
0 0 V3
=
H11 H12 H13 Γm0
opH(H12)H22 H23 0 0
opH(H13)opH(H23) Σ 0 0
opH(Γm) 0 0 0 0
0 0 0 0 0
,
ˆ
U=U10
0Irm−1 Ip0
0U2
U30 0
0Iµ0
0 0 V3
,
U=
In1+...+nm−10 0
0ˆ
U0
0 0 Iq1+...+qm−1
.
4
Set
N=
p−τ µ
p−τN22 0
µ0 0 ,H=
p−τ µ
p−τH22 H23
µ opH(H23) Σ ∈Rl,l,
and l=p−τ+µ.
END IF
END IF
Form H=UTHU,N=UTNU, and U=UU.
END WHILE
Algorithm 2.2 is implemented in the function STCSSP.
Remark 1. When we compute the form (5) for an even pencil with the help of Algorithm 2.2
we can determine the structured Kronecker canonical form of Thompson [3], for the original
even pencil (N, H)with the help of the values computed by the algorithm, see [1, Theorem
3.3].
Let us end this section with two graphs that give an overview of the subroutines that
come with STCSSP.
Fig.1: STCSSP and its subroutines. An arrow pointing from Ato Bmeans that routine A
may call subroutine Bat some point.
5
DBDSQX slight adaption of the LAPACK routine DBDSQR,
which computes a singular value decomposition
DLANSK skew-symmetric modification of the LAPACK routine DLANSY
DLASKE skew-symmetric modification of the LAPACK routine DLATRD
DOSKRD compute the eigenvectors and eigenvalues of a real
skew-symmetric tridiagonal matrix with the help of
DBDSQX
DOSTQR slight adaption of the LAPACK routine DSTEQR
DSKMV skew-symmetric modification of the BLAS routine DSYMV
DSKR2 skew-symmetric modification of the BLAS routine DSYR2
DSKR2K skew-symmetric modification of the BLAS routine DSYR2K
DSKSZ split off the odd dimension from a
skew-symmetric matrix
DSKTD2 skew-symmetric modification of the LAPACK routine DSYTD2
DSKTRD skew-symmetric modification of the LAPACK routine DSYTRD
MB01SS adaption of MB01SW (see below) for skew-symmetric matrices
MB01ST adaption of MB01SU (see below) for skew-symmetric matrices
MB01SU a slight adaption of the SLICOT routine MB01RU (uses BLAS3)
MB01SW a copy of the SLICOT routine MB01RW (uses BLAS2)
OTRSM apply (a part of) a block-diagonal Orthogonal
equivalence TRansformation to a (skew-)Symmetric
Matrix
OTRSZ apply (a part of) a block-diagonal Orthogonal
equivalence TRansformation to a (skew-)Symmetric
matrix, which contains some Zero blocks
SFSKEW returns the real Schur Form of a
real SKEW-symmetric matrix
SFSTUP once SFSKEW or SFSYM has finished its
computations, this routine SeTsUP a matrix
therefrom
SFSYM2 computes an ordered real Schur Form of a
real SYMmetric matrix
SFSYM returns an ordered real Schur Form of a
real SYMmetric matrix
Tab.1: The functionality implemented in the subroutines.
6
3 Specification
SUBROUTINE STCSSP( SYMN, SYMH, UPLON, UPLOH, COMPZ, NDIM, N, LDN,
$ H, LDH, U, LDU, M, PINVEC, NUNVEC, PIVEC,
$ NUVEC, NVEC, QVEC, P, L, TOL, DWORK, LDWORK,
$ INFO )
C .. Scalar Arguments ..
CHARACTER*1 SYMN, SYMH, UPLON, UPLOH, COMPZ
INTEGER NDIM, LDN, LDH, LDU, INFO, LDWORK, M, P, L
DOUBLE PRECISION TOL
C .. Array Arguments ..
DOUBLE PRECISION N( LDN, * ), H( LDH, * ), U( LDU, * ),
$ DWORK( * )
INTEGER PINVEC( * ), NUNVEC( * ), PIVEC( * ),
$ NUVEC( * ), NVEC( * ), QVEC( * )
4 Argument List
4.1 Mode Parameters
SYMN, SYMH -CHARACTER*1
These parameters define which of the cases (1a), (1b), (2a), or (2b) is considered.
Each of the two parameters may thereby be set to either ’S’, which means that
the corresponding matrix is assumed to be symmetric, or to ’N’, which means that
the corresponding matrix is not assumed to be symmetric but it is assumed to be
skew-symmetric. SYMN refers to the matrix Nand SYMH refers to the matrix H,
as in Algorithm 2.2. Other values are not allowed and will result in an erroneous
termination of the algorithm. For example, to compute the structured staircase
form of an even pencil, set SYMN = ’N’ and SYMH = ’S’.
With the operators defined in (3) one could also say that SYMN = ’N’ (or SYMH =
’N’, resp.) means opN(N) = −NT(or opH(H) = −HT, resp.) and SYMN = ’S’
(or SYMH = ’S’, resp.) means opN(N) = NT(or opH(H) = HT, resp.).
Since the diagonal of a real skew-symmetric matrices is always zero, it does not
have to be stored. This is why the diagonal elements in the array N(or H, resp.)
are not considered in the computation, once SYMN = ’N’ (or SYMH = ’N’, resp.).
UPLON, UPLOH -CHARACTER*1
These parameters tell the algorithm how the matrices Nand H(as in Algorithm
2.2) are given on input and also how the staircase form is to be stored, once the
algorithm has finished successfully. Each of the two parameters may thereby be
set to either ’U’, which means that the upper triangular part of the array N(or
H, resp.) is assumed to hold the upper triangular part of the matrix N(or H,
resp.), or to ’L’, which means that the lower triangular part of the array N(or H,
resp.) is assumed to hold the lower triangular part of the matrix N(or H, resp.).
UPLON refers to the matrix Nand UPLOH refers to the matrix H, as in Algorithm
2.2. Other values are not allowed and will result in an erroneous termination of
7
the algorithm. The unused part of the arrays may be used as workspace, but the
values in this unused part are not considered. Note, that if one of the parameters
SYMN or SYMH is set to ’N’ then only the strict upper or lower, respectively, part of
the corresponding array is used for computation (although the rest may be used
as workspace). These parameters also control how the staircase form is returned
in the case of an successful exit. The same regions in the arrays Nand Hthat
matter on entry contain the matrices that are transformed to staircase form, on
exit.
COMPZ -CHARACTER*1
This parameter controls, whether the orthogonal transformation Uas in Algorithm
2.2 shall also be computed or not. It is faster not to compute the transformation.
Setting this parameter to ’N’ does not compute the orthogonal transformation.
In this case the array Uis not referenced, and can be supplied as a dummy array
(i.e. set LDU = 1 and declare this array to be U(1,1) in the calling program).
Setting this parameter to ’V’ causes the algorithm to compute the orthogonal
transformation into the array U.
4.2 Input/Output Parameters
NDIM -(input) INTEGER
The order of the matrices N and H, and thus the dimension of the problem. NDIM
has to be greater or equal to 0, otherwise the algorithm will return in failure.
N-(input/output) DOUBLE PRECISION, array (LDN, NDIM)
On entry, the leading NDIM-by-NDIM part of the array has to contain the matrix
N(as in Algorithm 2.2), according to the parameters SYMN and UPLON.
On successful exit, the leading NDIM-by-NDIM part of this array contains a part of
the staircase form of the matrix N, according to Algorithm 2.2. The parameters
SYMN and UPLON also determine which part of this array contains the actual data. If
UPLON = ’U’, then only the upper triangular part contains the upper triangular
part of the staircase form and if UPLON = ’L’ only the lower triangular part
contains the lower triangular part of the staircase form. Also, if SYMN = ’N’,
then the diagonal in the array Nhas no meaning, on exit.
LDN -(input) INTEGER
The leading dimension of the array N. The parameter LDN has to be greater or
equal to MAX(1,NDIM), otherwise the algorithm will return in failure.
H-(input/output) DOUBLE PRECISION, array (LDH, NDIM)
On entry, the leading NDIM-by-NDIM part of the array has to contain the matrix
H(as in Algorithm 2.2), according to the parameters SYMH and UPLOH.
On successful exit, the leading NDIM-by-NDIM part of this array contains a part of
the staircase form of the matrix H, according to Algorithm 2.2. The parameters
SYMH and UPLOH also determine which part of this array contains the actual data. If
8
UPLOH = ’U’, then only the upper triangular part contains the upper triangular
part of the staircase form and if UPLOH = ’L’ only the lower triangular part
contains the lower triangular part of the staircase form. Also, if SYMH = ’N’,
then the diagonal in the array Hhas no meaning, on exit.
LDH -(input) INTEGER
The leading dimension of the array H. The parameter LDH has to be greater or
equal to MAX(1,NDIM), otherwise the algorithm will return in failure.
U-(output) DOUBLE PRECISION, array (LDU, NDIM)
If COMPZ = ’V’ and the algorithm terminated successfully, then the leading NDIM-
by-NDIM part of this array contains the orthogonal transformation matrix Uas
computed by Algorithm 2.2, which has been used to reduce the original matrices
Nand Hto structured staircase form (5). If COMPZ = ’N’, then Uis not referenced
and can be supplied as a dummy array (i.e. set LDU = 1 and declare this array to
be U(1,1) in the calling program).
LDU -(input) INTEGER
The leading dimension of the array U.LDU has to be greater or equal to 1 in any
case and also to be greater or equal to NDIM, if COMPZ = ’V’.
M-(output) INTEGER
The number of reduction steps that where necessary to reveal the Kronecker struc-
ture. Mis always greater or equal to 0 and is corresponds to min Algorithm 2.2.
PINVEC, NUNVEC -(output) INTEGER, array (NDIM+1)
On exit, with INFO = 0 and SYMN = ’S’, the first M+1 entries of the arrays contain
the inertia indices corresponding to the Nmatrix, analogously to the πN
i’s and
νN
i’s in Algorithm 2.2. PINVEC(M+1) and NUNVEC(M+1) may not be used. In this
case they are both zero. This happens when the algorithm does not exit because
a submatrix of Nwith full rank was discovered (i.e. the algorithm did not exit
from (6) ). On exit, with INFO = 0 and SYMN = ’N’, the first M+1 entries of the
array contain only zeros.
PIVEC, NUVEC -(output) INTEGER, array (NDIM)
On exit, with INFO = 0 and SYMH = ’S’, the first Mentries of the arrays contain
the inertia indices corresponding to the Hmatrix, analogously to the πi’s and νi’s
in Algorithm 2.2. On exit, with INFO = 0 and SYMH = ’N’, the first Mentries of
the array contain only zeros.
NVEC, QVEC -(output) INTEGER, array (NDIM)
On exit, with INFO = 0, the first Mentries of the arrays contain the ni’s and qi’s as
in Theorem 2.1 and Algorithm 2.2. Thus, these arrays describe the block structure
of the structured staircase form.
P-(output) INTEGER
9
On exit, with INFO = 0, this integer contains the number of the finite eigenvalues
of the pencil and thus the number of finite eigenvalues of the regular, index 1 part
of the pencil. We have 0≤P≤L, on successful exit. This parameter corresponds
to the value pin Algorithm 2.2 and Theorem 2.1.
L-(output) INTEGER
On exit, with INFO = 0, this integer contains the size of the regular, index 1 part of
the pencil. We have P≤L≤NDIM on successful exit. This parameter corresponds
to the value lin Algorithm 2.2 and Theorem 2.1.
4.3 Tolerance
TOL -(input) DOUBLE PRECISION
Absolute value, below which an eigenvalue or singular value shall be considered
zero. If TOL <= 0.0 is given on entry, the tolerance is automatically chosen to be
NDIM·eps, where eps is the machine precision, as returned by the LAPACK routine
DLAMCH.
In a later version of STCSSP we will implement the rank decision algorithm de-
scribed in [2, Section 5] which involves an additional input parameter GAP. This
rank decision algorithm also uses different tolerances for rank decisions in Nand
Hdepending on the Frobenius norms of Nand H.
4.4 Workspace
DWORK -DOUBLE PRECISION, array (LDWORK)
On exit, if INFO = 0,DWORK(1) returns the optimal value of LDWORK.
LDWORK -INTEGER
The length of the array DWORK that may be used by the algorithm as workspace.
LDWORK has to be greater or equal to NDIM*NDIM + 3*NDIM + 3, otherwise the
algorithm will issue an error message (i.e., the algorithm will return with INFO =
-22).
4.5 Error Indicator
INFO -INTEGER
INFO = 0: Successful exit.
INFO < 0: If INFO = -i, the i-th argument had an illegal value.
INFO = 1: Calculating the (skew-)symmetric (which one was computed de-
pends on the parameter SYMN) Schur-form of a part of the matrix
N failed.
INFO = 2: Calculating the (skew-)symmetric (which one was computed de-
pends on the parameter SYMH) Schur-form of a part of the matrix
H failed.
INFO = 3: The LAPACK routine DGESVD returned with an info value greater than
zero, i.e., DGESVD was not able to compute a SVD.
10
References
[1] R. Byers, V. Mehrmann, and H. Xu. A structured staircase algorithm for skew-
symmetric/symmetric pencils. Electr. Trans. Num. Anal., 26:1–33, 2007.
[2] J. W. Demmel and B. K˚agstr¨om. The generalized Schur decomposition of an arbitrary
pencil A - λB: Robust software with error bounds and applications. Part II: Software and
applications. ACM Trans. Math. Software, 19(2):175–201, 1993.
[3] R. C. Thompson. Pencils of complex and real symmetric and skew matrices. Linear
Algebra Appl., 147:323–371, 1991.
11
5 Example
We consider the following even matrix pencil which is taken from [1].
αN −βH =α
0 1 0 0 0
−1 0 0 0 0
0 0 0 0 0
0 0 0 0 1
0 0 0 −1 0
−β
00100
01000
10000
00010
00004
(7)
5.1 Program Text
PROGRAM STCSSP_example
implicit none
C
C DEMO : Demonstration program for STCSSP
C
*
* .. Parameters ..
INTEGER NIN , NOUT
PARAMETER ( NIN = 5, NOUT = 6 )
INTEGER NMAX
PARAMETER ( NMAX = 20 )
INTEGER LDN , LDH , LDU
PARAMETER ( LDN = NMAX , LDH = NMAX , LDU = NMAX )
INTEGER LDWORK
PARAMETER ( LDWORK = NMAX * NMAX + 3* NMAX + 3 )
* .. Local Scalars ..
CHARACTER *1 SYMN , SYMH , UPLON , UPLOH , COMPZ
INTEGER NDIM , M, P, L
INTEGER PINVEC ( NMAX +1) , NUNVEC (NMAX +1) ,
1 PIVEC(NMAX), NUVEC(NMAX)
INTEGER NVEC( NMAX ), QVEC ( NMAX)
DOUBLE PRECISION N(NMAX ,NMAX ), H(NMAX , NMAX ), U(NMAX , NMAX)
DOUBLE PRECISION TOL1
DOUBLE PRECISION DWORK ( LDWORK )
INTEGER INFO
INTEGER I, J
* .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
* .. External Subroutines ..
EXTERNAL STCSSP
*
WRITE ( NOUT , FMT = 99999 )
* Skip the heading in the data file and read the data .
READ ( NIN , FMT = ’() ’ )
READ ( NIN , FMT = * ) SYMN , SYMH , UPLON , UPLOH , COMPZ
READ ( NIN , FMT = * ) NDIM
READ ( NIN , FMT = * ) TOL1
IF ( NDIM.LT .0 .OR. NDIM .GT.NMAX ) THEN
WRITE ( NOUT , FMT = 99991 ) N
ELSE
READ ( NIN , FMT = * ) ( ( N(I,J), J = 1,NDIM ), I = 1, NDIM )
READ ( NIN , FMT = * ) ( ( H(I,J), J = 1,NDIM ), I = 1, NDIM )
C
12
C Run the structured staircase algorithm
C
CALL STCSSP ( SYMN , SYMH , UPLON , UPLOH , COMPZ , NDIM , N, LDN ,
$ H, LDH , U, LDU , M, PINVEC , NUNVEC , PIVEC ,
$ NUVEC , NVEC , QVEC , P, L, TOL1 , DWORK , LDWORK ,
$ INFO )
C
C Complete the matrices
C
CALL SUPMAT ( SYMN , UPLON , N , LDN , NDIM )
CALL SUPMAT ( SYMH , UPLOH , H , LDH , NDIM )
C
C Show output to console
C
IF ( INFO .NE .0 ) THEN
WRITE ( NOUT , FMT = 99998 ) INFO
ELSE
WRITE ( NOUT , FMT = 99997 ) P, L, M
WRITE ( NOUT , FMT = 99996 )
DO 20 I = 1, NDIM
WRITE ( NOUT , FMT = 99993 ) ( N(I,J), J = 1, NDIM )
20 CONTINUE
WRITE ( NOUT , FMT = 99995 )
DO 30 I = 1, NDIM
WRITE ( NOUT , FMT = 99993 ) ( H(I,J), J = 1, NDIM )
30 CONTINUE
IF( LSAME (COMPZ , ’V’) ) THEN
WRITE ( NOUT , FMT = 99994 )
DO 40 I = 1, NDIM
WRITE ( NOUT , FMT = 99993 ) ( U(I,J), J = 1, NDIM )
40 CONTINUE
END IF
WRITE ( NOUT , FMT = 99990 )
DO 50 I = 1,M+1
WRITE ( NOUT , FMT = 99980 ) I, PINVEC (I), NUNVEC (I)
50 CONTINUE
WRITE ( NOUT , FMT = 99989 )
DO 60 I = 1,M
WRITE ( NOUT , FMT = 99980 ) I, PIVEC (I), NUVEC (I)
60 CONTINUE
WRITE ( NOUT , FMT = 99988 )
DO 70 I = 1,M
WRITE ( NOUT , FMT = 99980 ) I, NVEC(I), QVEC(I)
70 CONTINUE
END IF
END IF
STOP
C
99999 FORMAT (’ STCSSP EXAMPLE PROGRAM ’ ,/1X,
1 ’Computing staircase form of the pencil alpha N - beta H’)
99998 FORMAT (’ STCSPP returned in failure with INFO = ’,I2)
99997 FORMAT (/ ’ P = ’,I3 ,’ ; L = ’,I3 ,’ ; M = ’,I3)
99996 FORMAT (/ ’ The staircase form of matrix N is ’)
99995 FORMAT (/ ’ The staircase form of matrix H is ’)
99994 FORMAT (/ ’ The orthogonal transformation U is ’)
99993 FORMAT (20(1X,F8.4))
99992 FORMAT (/’ ’ )
13
99991 FORMAT (/ ’ Dimension of the problem is out of range .’,/’ N = ’,I5)
99990 FORMAT (/ ’ The characteristic values for N are ’ ,/’ i’ ,6X,
1 ’PI(i )’ ,6X, ’NU (i) ’)
99989 FORMAT (/ ’ The characteristic values for H are ’ ,/’ i’ ,6X,
1 ’PI(i )’ ,6X, ’NU (i) ’)
99988 FORMAT (/ ’ The dimension of the blocks are ’ ,/’ i’ ,6X,
1 ’N(i) ’ ,6X ,’Q(i) ’)
99980 FORMAT (I4 ,3X ,I7 ,3X ,I7 )
END PROGRAM STCSSP_example
C
C Subroutine to complete the returned matrix in all positions
C
SUBROUTINE SUPMAT ( SYM , UPLO , A, LDA , N )
implicit none
C .. Parameters ..
DOUBLE PRECISION ZERO , ONE
PARAMETER ( ZERO = 0.0 D0 , ONE = 1.0 D0 )
C .. Arguments ..
CHARACTER *1 SYM , UPLO
DOUBLE PRECISION A(LDA ,*), FACT
INTEGER LDA , N, I, J
C .. External Functions ..
LOGICAL LSAME
EXTERNAL LSAME
C
C Handle the strict upper / lower triangular part
C
IF( LSAME ( SYM , ’S’ ) ) THEN
FACT = ONE
ELSE
FACT = -ONE
END IF
IF( LSAME ( UPLO , ’U’ ) ) THEN
DO 200 I =1 ,N -1
DO 100 J=I+1,N
A(J,I) = FACT * A(I,J)
100 CONTINUE
200 CONTINUE
ELSE
DO 400 I =1 ,N -1
DO 300 J=I+1,N
A(I,J) = FACT * A(J,I)
300 CONTINUE
400 CONTINUE
END IF
C
C Handle the diagonal
C
IF( LSAME ( SYM , ’N’ ) ) THEN
DO 500 I=1, N
A(I ,I) = ZERO
500 CONTINUE
END IF
C
END SUBROUTINE SUPMAT
5.2 Program Data
14
STCSSP Example Program Data
NSULV
5
1e -12
-7 1 0 0 0
-7 -7 0 0 0
-7 -7 -7 0 0
-7 -7 -7 -7 1
-7 -7 -7 -7 -7
0 -7 -7 -7 -7
0 1 -7 -7 -7
1 0 0 -7 -7
0 0 0 1 -7
00004
5.3 Program Results
STCSSP EXAMPLE PROGRAM
Computing staircase form of the pencil alpha N - beta H
P = 2 ; L = 3 ; M = 2
The staircase form of matrix N is
0.0000 0.0000 0.0000 1.0000 0.0000
0.0000 0.0000 -1.0000 0.0000 0.0000
0.0000 1.0000 0.0000 0.0000 0.0000
-1.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000
The staircase form of matrix H is
0.0000 0.0000 0.0000 0.0000 1.0000
0.0000 1.0000 0.0000 0.0000 0.0000
0.0000 0.0000 4.0000 0.0000 0.0000
0.0000 0.0000 0.0000 1.0000 0.0000
1.0000 0.0000 0.0000 0.0000 0.0000
The orthogonal transformation U is
-1.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 -1.0000 0.0000
0.0000 0.0000 0.0000 0.0000 -1.0000
0.0000 -1.0000 0.0000 0.0000 0.0000
0.0000 0.0000 1.0000 0.0000 0.0000
The characteristic values for N are
i PI(i) NU (i)
1 0 0
2 0 0
3 0 0
The characteristic values for H are
i PI(i) NU (i)
1 0 0
2 1 0
The dimension of the blocks are
i N(i) Q(i)
1 1 1
2 0 0
15
