PSTRORD(3) ScaLAPACK routine of NEC Numeric Library Collection PSTRORD(3)
NAME
PSTRORD - reorders the real Schur factorization of a real matrix A =
Q*T*Q**T, so that a selected cluster of eigenvalues appears in the
leading diagonal blocks of the upper quasi-triangular matrix T, and the
leading columns of Q form an orthonormal basis of the corresponding
right invariant subspace
SYNOPSIS
SUBROUTINE PSTRORD( COMPQ, SELECT, PARA, N, T, IT, JT, DESCT, Q, IQ,
JQ, DESCQ, WR, WI, M, WORK, LWORK, IWORK, LIWORK, INFO )
CHARACTER COMPQ
INTEGER INFO, LIWORK, LWORK, M, N, IT, JT, IQ, JQ
INTEGER SELECT( * )
INTEGER PARA( 6 ), DESCT( * ), DESCQ( * ), IWORK( * )
REAL Q( * ), T( * ), WI( * ), WORK( * ), WR( * )
PURPOSE
PSTRORD reorders the real Schur factorization of a real matrix A =
Q*T*Q**T, so that a selected cluster of eigenvalues appears in the
leading diagonal blocks of the upper quasi-triangular matrix T, and the
leading columns of Q form an orthonormal basis of the corresponding
right invariant subspace.
T must be in Schur form (as returned by PSLAHQR), that is, block upper
triangular with 1-by-1 and 2-by-2 diagonal blocks.
This subroutine uses a delay and accumulate procedure for performing
the off-diagonal updates (see references for details).
Notes
=====
Each global data object is described by an associated description vec-
tor. This vector stores the information required to establish the map-
ping between an object element and its corresponding process and memory
location.
Let A be a generic term for any 2D block cyclicly distributed array.
Such a global array has an associated description vector DESCA. In the
following comments, the character _ should be read as "of the global
array".
NOTATION STORED IN EXPLANATION
--------------- -------------- --------------------------------------
DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
DTYPE_A = 1.
CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
the BLACS process grid A is distribu-
ted over. The context itself is glo-
bal, but the handle (the integer
value) may vary.
M_A (global) DESCA( M_ ) The number of rows in the global
array A.
N_A (global) DESCA( N_ ) The number of columns in the global
array A.
MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
the rows of the array.
NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
the columns of the array.
RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
row of the array A is distributed.
CSRC_A (global) DESCA( CSRC_ ) The process column over which the
first column of the array A is
distributed.
LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
array. LLD_A >= MAX(1,LOCr(M_A)).
Let K be the number of rows or columns of a distributed matrix, and
assume that its process grid has dimension p x q.
LOCr( K ) denotes the number of elements of K that a process would
receive if K were distributed over the p processes of its process col-
umn.
Similarly, LOCc( K ) denotes the number of elements of K that a process
would receive if K were distributed over the q processes of its process
row.
The values of LOCr() and LOCc() may be determined via a call to the
ScaLAPACK tool function, NUMROC:
LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
An upper bound for these quantities may be computed by:
LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
ARGUMENTS
COMPQ (global input) CHARACTER*1
= 'V': update the matrix Q of Schur vectors;
= 'N': do not update Q.
SELECT (global input/output) INTEGER array, dimension (N)
SELECT specifies the eigenvalues in the selected cluster. To
select a real eigenvalue w(j), SELECT(j) must be set to 1.
To select a complex conjugate pair of eigenvalues w(j) and
w(j+1), corresponding to a 2-by-2 diagonal block, either
SELECT(j) or SELECT(j+1) or both must be set to 1;
a complex conjugate pair of eigenvalues must be either both
included in the cluster or both excluded.
On output, the (partial) reordering is displayed.
PARA (global input) INTEGER*6
Block parameters (some should be replaced by calls to PILAENV
and others by meaningful default values):
PARA(1) = maximum number of concurrent computational windows
allowed in the algorithm;
0 < PARA(1) <= min(NPROW,NPCOL) must hold;
PARA(2) = number of eigenvalues in each window;
0 < PARA(2) < PARA(3) must hold;
PARA(3) = window size; PARA(2) < PARA(3) < DESCT(MB_)
must hold;
PARA(4) = minimal percentage of flops required for
performing matrix-matrix multiplications instead
of pipelined orthogonal transformations;
0 <= PARA(4) <= 100 must hold;
PARA(5) = width of block column slabs for row-wise
application of pipelined orthogonal
transformations in their factorized form;
0 < PARA(5) <= DESCT(MB_) must hold.
PARA(6) = the maximum number of eigenvalues moved together
over a process border; in practice, this will be
approximately half of the cross border window size
0 < PARA(6) <= PARA(2) must hold;
N (global input) INTEGER
The order of the globally distributed matrix T. N >= 0.
T (local input/output) REAL array,
dimension (LLD_T,LOCc(N)).
On entry, the local pieces of the global distributed upper
quasi-triangular matrix T, in Schur form. On exit, T is over-
written by the local pieces of the reordered matrix T, again in
Schur form, with the selected eigenvalues in the globally lead-
ing diagonal blocks.
IT (global input) INTEGER
JT (global input) INTEGER
The row and column index in the global array T indicating the
first column of sub( T ). IT = JT = 1 must hold.
DESCT (global and local input) INTEGER array of dimension DLEN_.
The array descriptor for the global distributed matrix T.
Q (local input/output) REAL array, dimension (LLD_Q,LOCc(N)).
On entry, if COMPQ = 'V', the local pieces of the global dis-
tributed matrix Q of Schur vectors.
On exit, if COMPQ = 'V', Q has been postmultiplied by the
global orthogonal transformation matrix which reorders T; the
leading M columns of Q form an orthonormal basis for the speci-
fied invariant subspace.
If COMPQ = 'N', Q is not referenced.
IQ (global input) INTEGER
JQ (global input) INTEGER
The column index in the global array Q indicating the first
column of sub( Q ). IQ = JQ = 1 must hold.
DESCQ (global and local input) INTEGER array of dimension DLEN_.
The array descriptor for the global distributed matrix Q.
WR (global output) REAL array, dimension (N)
WI (global output) REAL array, dimension (N)
The real and imaginary parts, respectively, of the reordered
eigenvalues of T. The eigenvalues are in principle stored in
the same order as on the diagonal of T, with WR(i) = T(i,i)
and, if T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0
and WI(i+1) = -WI(i).
Note also that if a complex eigenvalue is sufficiently ill-con-
ditioned, then its value may differ significantly from its
value before reordering.
M (global output) INTEGER
The dimension of the specified invariant subspace.
0 <= M <= N.
WORK (local workspace/output) REAL array,
dimension (LWORK)
On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
LWORK (local input) INTEGER
The dimension of the array WORK.
If LWORK = -1, then a workspace query is assumed; the routine
only calculates the optimal size of the WORK array, returns
this value as the first entry of the WORK array, and no error
message related to LWORK is issued by PXERBLA.
IWORK (local workspace/output) INTEGER array, dimension (LIWORK)
LIWORK (local input) INTEGER
The dimension of the array IWORK.
If LIWORK = -1, then a workspace query is assumed; the routine
only calculates the optimal size of the IWORK array, returns
this value as the first entry of the IWORK array, and no error
message related to LIWORK is issued by PXERBLA.
INFO (global output) INTEGER
= 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value.
If the i-th argument is an array and the j-entry had an illegal
value, then INFO = -(i*1000+j), if the i-th argument is a
scalar and had an illegal value, then INFO = -i.
> 0: here we have several possibilites
*) Reordering of T failed because some eigenvalues are too
close to separate (the problem is very ill-conditioned);
T may have been partially reordered, and WR and WI
contain the eigenvalues in the same order as in T.
On exit, INFO = {the index of T where the swap failed}.
*) A 2-by-2 block to be reordered split into two 1-by-1
blocks and the second block failed to swap with an
adjacent block.
On exit, INFO = {the index of T where the swap failed}.
*) If INFO = N+1, there is no valid BLACS context (see the
BLACS documentation for details).
In a future release this subroutine may distinguish between the
case 1 and 2 above.
Additional requirements
=======================
The following alignment requirements must hold:
(a) DESCT( MB_ ) = DESCT( NB_ ) = DESCQ( MB_ ) = DESCQ( NB_ )
(b) DESCT( RSRC_ ) = DESCQ( RSRC_ )
(c) DESCT( CSRC_ ) = DESCQ( CSRC_ )
All matrices must be blocked by a block factor larger than or equal to
two (3). This is to simplify reordering across processor borders in the
presence of 2-by-2 blocks.
Limitations
===========
This algorithm cannot work on submatrices of T and Q, i.e., IT = JT =
IQ = JQ = 1 must hold. This is however no limitation since PDLAHQR does
not compute Schur forms of submatrices anyway.
Parallel execution recommendations
==================================
Use a square grid, if possible, for maximum performance. The block
parameters in PARA should be kept well below the data distribution
block size.
In general, the parallel algorithm strives to perform as much work as
possible without crossing the block borders on the main block diagonal.
Keywords
========
Real Schur form, eigenvalue reordering
ScaLAPACK routine 31 October 2017 PSTRORD(3)