ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdtrsen.f
Go to the documentation of this file.
1  SUBROUTINE pdtrsen( JOB, COMPQ, SELECT, PARA, N, T, IT, JT,
2  $ DESCT, Q, IQ, JQ, DESCQ, WR, WI, M, S, SEP, WORK, LWORK,
3  $ IWORK, LIWORK, INFO )
4 *
5 * Contribution from the Department of Computing Science and HPC2N,
6 * Umea University, Sweden
7 *
8 * -- ScaLAPACK computational routine (version 2.0.1) --
9 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
10 * Univ. of Colorado Denver and University of California, Berkeley.
11 * January, 2012
12 *
13  IMPLICIT NONE
14 *
15 * .. Scalar Arguments ..
16  CHARACTER COMPQ, JOB
17  INTEGER INFO, LIWORK, LWORK, M, N,
18  $ it, jt, iq, jq
19  DOUBLE PRECISION S, SEP
20 * ..
21 * .. Array Arguments ..
22  LOGICAL SELECT( N )
23  INTEGER PARA( 6 ), DESCT( * ), DESCQ( * ), IWORK( * )
24  DOUBLE PRECISION Q( * ), T( * ), WI( * ), WORK( * ), WR( * )
25 * ..
26 *
27 * Purpose
28 * =======
29 *
30 * PDTRSEN reorders the real Schur factorization of a real matrix
31 * A = Q*T*Q**T, so that a selected cluster of eigenvalues appears
32 * in the leading diagonal blocks of the upper quasi-triangular matrix
33 * T, and the leading columns of Q form an orthonormal basis of the
34 * corresponding right invariant subspace. The reordering is performed
35 * by PDTRORD.
36 *
37 * Optionally the routine computes the reciprocal condition numbers of
38 * the cluster of eigenvalues and/or the invariant subspace. SCASY
39 * library is needed for condition estimation.
40 *
41 * T must be in Schur form (as returned by PDLAHQR), that is, block
42 * upper triangular with 1-by-1 and 2-by-2 diagonal blocks.
43 *
44 * Notes
45 * =====
46 *
47 * Each global data object is described by an associated description
48 * vector. This vector stores the information required to establish
49 * the mapping between an object element and its corresponding process
50 * and memory location.
51 *
52 * Let A be a generic term for any 2D block cyclicly distributed array.
53 * Such a global array has an associated description vector DESCA.
54 * In the following comments, the character _ should be read as
55 * "of the global array".
56 *
57 * NOTATION STORED IN EXPLANATION
58 * --------------- -------------- --------------------------------------
59 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
60 * DTYPE_A = 1.
61 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
62 * the BLACS process grid A is distribu-
63 * ted over. The context itself is glo-
64 * bal, but the handle (the integer
65 * value) may vary.
66 * M_A (global) DESCA( M_ ) The number of rows in the global
67 * array A.
68 * N_A (global) DESCA( N_ ) The number of columns in the global
69 * array A.
70 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
71 * the rows of the array.
72 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
73 * the columns of the array.
74 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
75 * row of the array A is distributed.
76 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
77 * first column of the array A is
78 * distributed.
79 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
80 * array. LLD_A >= MAX(1,LOCr(M_A)).
81 *
82 * Let K be the number of rows or columns of a distributed matrix,
83 * and assume that its process grid has dimension p x q.
84 * LOCr( K ) denotes the number of elements of K that a process
85 * would receive if K were distributed over the p processes of its
86 * process column.
87 * Similarly, LOCc( K ) denotes the number of elements of K that a
88 * process would receive if K were distributed over the q processes of
89 * its process row.
90 * The values of LOCr() and LOCc() may be determined via a call to the
91 * ScaLAPACK tool function, NUMROC:
92 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
93 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
94 * An upper bound for these quantities may be computed by:
95 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
96 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
97 *
98 * Arguments
99 * =========
100 *
101 * JOB (global input) CHARACTER*1
102 * Specifies whether condition numbers are required for the
103 * cluster of eigenvalues (S) or the invariant subspace (SEP):
104 * = 'N': none;
105 * = 'E': for eigenvalues only (S);
106 * = 'V': for invariant subspace only (SEP);
107 * = 'B': for both eigenvalues and invariant subspace (S and
108 * SEP).
109 *
110 * COMPQ (global input) CHARACTER*1
111 * = 'V': update the matrix Q of Schur vectors;
112 * = 'N': do not update Q.
113 *
114 * SELECT (global input) LOGICAL array, dimension (N)
115 * SELECT specifies the eigenvalues in the selected cluster. To
116 * select a real eigenvalue w(j), SELECT(j) must be set to
117 * .TRUE.. To select a complex conjugate pair of eigenvalues
118 * w(j) and w(j+1), corresponding to a 2-by-2 diagonal block,
119 * either SELECT(j) or SELECT(j+1) or both must be set to
120 * .TRUE.; a complex conjugate pair of eigenvalues must be
121 * either both included in the cluster or both excluded.
122 *
123 * PARA (global input) INTEGER*6
124 * Block parameters (some should be replaced by calls to
125 * PILAENV and others by meaningful default values):
126 * PARA(1) = maximum number of concurrent computational windows
127 * allowed in the algorithm;
128 * 0 < PARA(1) <= min(NPROW,NPCOL) must hold;
129 * PARA(2) = number of eigenvalues in each window;
130 * 0 < PARA(2) < PARA(3) must hold;
131 * PARA(3) = window size; PARA(2) < PARA(3) < DESCT(MB_)
132 * must hold;
133 * PARA(4) = minimal percentage of flops required for
134 * performing matrix-matrix multiplications instead
135 * of pipelined orthogonal transformations;
136 * 0 <= PARA(4) <= 100 must hold;
137 * PARA(5) = width of block column slabs for row-wise
138 * application of pipelined orthogonal
139 * transformations in their factorized form;
140 * 0 < PARA(5) <= DESCT(MB_) must hold.
141 * PARA(6) = the maximum number of eigenvalues moved together
142 * over a process border; in practice, this will be
143 * approximately half of the cross border window size
144 * 0 < PARA(6) <= PARA(2) must hold;
145 *
146 * N (global input) INTEGER
147 * The order of the globally distributed matrix T. N >= 0.
148 *
149 * T (local input/output) DOUBLE PRECISION array,
150 * dimension (LLD_T,LOCc(N)).
151 * On entry, the local pieces of the global distributed
152 * upper quasi-triangular matrix T, in Schur form. On exit, T is
153 * overwritten by the local pieces of the reordered matrix T,
154 * again in Schur form, with the selected eigenvalues in the
155 * globally leading diagonal blocks.
156 *
157 * IT (global input) INTEGER
158 * JT (global input) INTEGER
159 * The row and column index in the global array T indicating the
160 * first column of sub( T ). IT = JT = 1 must hold.
161 *
162 * DESCT (global and local input) INTEGER array of dimension DLEN_.
163 * The array descriptor for the global distributed matrix T.
164 *
165 * Q (local input/output) DOUBLE PRECISION array,
166 * dimension (LLD_Q,LOCc(N)).
167 * On entry, if COMPQ = 'V', the local pieces of the global
168 * distributed matrix Q of Schur vectors.
169 * On exit, if COMPQ = 'V', Q has been postmultiplied by the
170 * global orthogonal transformation matrix which reorders T; the
171 * leading M columns of Q form an orthonormal basis for the
172 * specified invariant subspace.
173 * If COMPQ = 'N', Q is not referenced.
174 *
175 * IQ (global input) INTEGER
176 * JQ (global input) INTEGER
177 * The column index in the global array Q indicating the
178 * first column of sub( Q ). IQ = JQ = 1 must hold.
179 *
180 * DESCQ (global and local input) INTEGER array of dimension DLEN_.
181 * The array descriptor for the global distributed matrix Q.
182 *
183 * WR (global output) DOUBLE PRECISION array, dimension (N)
184 * WI (global output) DOUBLE PRECISION array, dimension (N)
185 * The real and imaginary parts, respectively, of the reordered
186 * eigenvalues of T. The eigenvalues are in principle stored in
187 * the same order as on the diagonal of T, with WR(i) = T(i,i)
188 * and, if T(i:i+1,i:i+1) is a 2-by-2 diagonal block, WI(i) > 0
189 * and WI(i+1) = -WI(i).
190 * Note also that if a complex eigenvalue is sufficiently
191 * ill-conditioned, then its value may differ significantly
192 * from its value before reordering.
193 *
194 * M (global output) INTEGER
195 * The dimension of the specified invariant subspace.
196 * 0 <= M <= N.
197 *
198 * S (global output) DOUBLE PRECISION
199 * If JOB = 'E' or 'B', S is a lower bound on the reciprocal
200 * condition number for the selected cluster of eigenvalues.
201 * S cannot underestimate the true reciprocal condition number
202 * by more than a factor of sqrt(N). If M = 0 or N, S = 1.
203 * If JOB = 'N' or 'V', S is not referenced.
204 *
205 * SEP (global output) DOUBLE PRECISION
206 * If JOB = 'V' or 'B', SEP is the estimated reciprocal
207 * condition number of the specified invariant subspace. If
208 * M = 0 or N, SEP = norm(T).
209 * If JOB = 'N' or 'E', SEP is not referenced.
210 *
211 * WORK (local workspace/output) DOUBLE PRECISION array,
212 * dimension (LWORK)
213 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
214 *
215 * LWORK (local input) INTEGER
216 * The dimension of the array WORK.
217 *
218 * If LWORK = -1, then a workspace query is assumed; the routine
219 * only calculates the optimal size of the WORK array, returns
220 * this value as the first entry of the WORK array, and no error
221 * message related to LWORK is issued by PXERBLA.
222 *
223 * IWORK (local workspace/output) INTEGER array, dimension (LIWORK)
224 *
225 * LIWORK (local input) INTEGER
226 * The dimension of the array IWORK.
227 *
228 * If LIWORK = -1, then a workspace query is assumed; the
229 * routine only calculates the optimal size of the IWORK array,
230 * returns this value as the first entry of the IWORK array, and
231 * no error message related to LIWORK is issued by PXERBLA.
232 *
233 * INFO (global output) INTEGER
234 * = 0: successful exit
235 * < 0: if INFO = -i, the i-th argument had an illegal value.
236 * If the i-th argument is an array and the j-entry had
237 * an illegal value, then INFO = -(i*1000+j), if the i-th
238 * argument is a scalar and had an illegal value, then INFO = -i.
239 * > 0: here we have several possibilites
240 * *) Reordering of T failed because some eigenvalues are too
241 * close to separate (the problem is very ill-conditioned);
242 * T may have been partially reordered, and WR and WI
243 * contain the eigenvalues in the same order as in T.
244 * On exit, INFO = {the index of T where the swap failed}.
245 * *) A 2-by-2 block to be reordered split into two 1-by-1
246 * blocks and the second block failed to swap with an
247 * adjacent block.
248 * On exit, INFO = {the index of T where the swap failed}.
249 * *) If INFO = N+1, there is no valid BLACS context (see the
250 * BLACS documentation for details).
251 * *) If INFO = N+2, the routines used in the calculation of
252 * the condition numbers raised a positive warning flag
253 * (see the documentation for PGESYCTD and PSYCTCON of the
254 * SCASY library).
255 * *) If INFO = N+3, PGESYCTD raised an input error flag;
256 * please report this bug to the authors (see below).
257 * If INFO = N+4, PSYCTCON raised an input error flag;
258 * please report this bug to the authors (see below).
259 * In a future release this subroutine may distinguish between
260 * the case 1 and 2 above.
261 *
262 * Method
263 * ======
264 *
265 * This routine performs parallel eigenvalue reordering in real Schur
266 * form by parallelizing the approach proposed in [3]. The condition
267 * number estimation part is performed by using techniques and code
268 * from SCASY, see http://www.cs.umu.se/research/parallel/scasy.
269 *
270 * Additional requirements
271 * =======================
272 *
273 * The following alignment requirements must hold:
274 * (a) DESCT( MB_ ) = DESCT( NB_ ) = DESCQ( MB_ ) = DESCQ( NB_ )
275 * (b) DESCT( RSRC_ ) = DESCQ( RSRC_ )
276 * (c) DESCT( CSRC_ ) = DESCQ( CSRC_ )
277 *
278 * All matrices must be blocked by a block factor larger than or
279 * equal to two (3). This to simplify reordering across processor
280 * borders in the presence of 2-by-2 blocks.
281 *
282 * Limitations
283 * ===========
284 *
285 * This algorithm cannot work on submatrices of T and Q, i.e.,
286 * IT = JT = IQ = JQ = 1 must hold. This is however no limitation
287 * since PDLAHQR does not compute Schur forms of submatrices anyway.
288 *
289 * References
290 * ==========
291 *
292 * [1] Z. Bai and J. W. Demmel; On swapping diagonal blocks in real
293 * Schur form, Linear Algebra Appl., 186:73--95, 1993. Also as
294 * LAPACK Working Note 54.
295 *
296 * [2] Z. Bai, J. W. Demmel, and A. McKenney; On computing condition
297 * numbers for the nonsymmetric eigenvalue problem, ACM Trans.
298 * Math. Software, 19(2):202--223, 1993. Also as LAPACK Working
299 * Note 13.
300 *
301 * [3] D. Kressner; Block algorithms for reordering standard and
302 * generalized Schur forms, ACM TOMS, 32(4):521-532, 2006.
303 * Also LAPACK Working Note 171.
304 *
305 * [4] R. Granat, B. Kagstrom, and D. Kressner; Parallel eigenvalue
306 * reordering in real Schur form, Concurrency and Computations:
307 * Practice and Experience, 21(9):1225-1250, 2009. Also as
308 * LAPACK Working Note 192.
309 *
310 * Parallel execution recommendations
311 * ==================================
312 *
313 * Use a square grid, if possible, for maximum performance. The block
314 * parameters in PARA should be kept well below the data distribution
315 * block size. In particular, see [3,4] for recommended settings for
316 * these parameters.
317 *
318 * In general, the parallel algorithm strives to perform as much work
319 * as possible without crossing the block borders on the main block
320 * diagonal.
321 *
322 * Contributors
323 * ============
324 *
325 * Implemented by Robert Granat, Dept. of Computing Science and HPC2N,
326 * Umea University, Sweden, March 2007,
327 * in collaboration with Bo Kagstrom and Daniel Kressner.
328 * Modified by Meiyue Shao, October 2011.
329 *
330 * Revisions
331 * =========
332 *
333 * Please send bug-reports to granat@cs.umu.se
334 *
335 * Keywords
336 * ========
337 *
338 * Real Schur form, eigenvalue reordering, Sylvester matrix equation
339 *
340 * =====================================================================
341 * ..
342 * .. Parameters ..
343  CHARACTER TOP
344  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
345  $ lld_, mb_, m_, nb_, n_, rsrc_
346  DOUBLE PRECISION ZERO, ONE
347  PARAMETER ( TOP = '1-Tree',
348  $ block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
349  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
350  $ rsrc_ = 7, csrc_ = 8, lld_ = 9,
351  $ zero = 0.0d+0, one = 1.0d+0 )
352 * ..
353 * .. Local Scalars ..
354  LOGICAL LQUERY, WANTBH, WANTQ, WANTS, WANTSP
355  INTEGER ICOFFT12, ICTXT, IDUM1, IDUM2, IERR, ILOC1,
356  $ ipw1, iter, itt, jloc1, jtt, k, liwmin, lldt,
357  $ lldq, lwmin, mmax, mmin, myrow, mycol, n1, n2,
358  $ nb, noexsy, npcol, nprocs, nprow, space,
359  $ t12rows, t12cols, tcols, tcsrc, trows, trsrc,
360  $ wrk1, iwrk1, wrk2, iwrk2, wrk3, iwrk3
361  DOUBLE PRECISION DPDUM1, ELEM, EST, SCALE, RNORM
362 * .. Local Arrays ..
363  INTEGER DESCT12( DLEN_ ), MBNB2( 2 )
364 * ..
365 * .. External Functions ..
366  LOGICAL LSAME
367  INTEGER NUMROC
368  DOUBLE PRECISION PDLANGE
369  EXTERNAL lsame, numroc, pdlange
370 * ..
371 * .. External Subroutines ..
372  EXTERNAL blacs_gridinfo, chk1mat, descinit,
373  $ igamx2d, infog2l, pdlacpy, pdtrord, pxerbla,
374  $ pchk1mat, pchk2mat
375 * $ , PGESYCTD, PSYCTCON
376 * ..
377 * .. Intrinsic Functions ..
378  INTRINSIC max, min, sqrt
379 * ..
380 * .. Executable Statements ..
381 *
382 * Get grid parameters
383 *
384  ictxt = desct( ctxt_ )
385  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
386  nprocs = nprow*npcol
387 *
388 * Test if grid is O.K., i.e., the context is valid
389 *
390  info = 0
391  IF( nprow.EQ.-1 ) THEN
392  info = n+1
393  END IF
394 *
395 * Check if workspace
396 *
397  lquery = lwork.EQ.-1 .OR. liwork.EQ.-1
398 *
399 * Test dimensions for local sanity
400 *
401  IF( info.EQ.0 ) THEN
402  CALL chk1mat( n, 5, n, 5, it, jt, desct, 9, info )
403  END IF
404  IF( info.EQ.0 ) THEN
405  CALL chk1mat( n, 5, n, 5, iq, jq, descq, 13, info )
406  END IF
407 *
408 * Check the blocking sizes for alignment requirements
409 *
410  IF( info.EQ.0 ) THEN
411  IF( desct( mb_ ).NE.desct( nb_ ) ) info = -(1000*9 + mb_)
412  END IF
413  IF( info.EQ.0 ) THEN
414  IF( descq( mb_ ).NE.descq( nb_ ) ) info = -(1000*13 + mb_)
415  END IF
416  IF( info.EQ.0 ) THEN
417  IF( desct( mb_ ).NE.descq( mb_ ) ) info = -(1000*9 + mb_)
418  END IF
419 *
420 * Check the blocking sizes for minimum sizes
421 *
422  IF( info.EQ.0 ) THEN
423  IF( n.NE.desct( mb_ ) .AND. desct( mb_ ).LT.3 )
424  $ info = -(1000*9 + mb_)
425  IF( n.NE.descq( mb_ ) .AND. descq( mb_ ).LT.3 )
426  $ info = -(1000*13 + mb_)
427  END IF
428 *
429 * Check parameters in PARA
430 *
431  nb = desct( mb_ )
432  IF( info.EQ.0 ) THEN
433  IF( para(1).LT.1 .OR. para(1).GT.min(nprow,npcol) )
434  $ info = -(1000 * 4 + 1)
435  IF( para(2).LT.1 .OR. para(2).GE.para(3) )
436  $ info = -(1000 * 4 + 2)
437  IF( para(3).LT.1 .OR. para(3).GT.nb )
438  $ info = -(1000 * 4 + 3)
439  IF( para(4).LT.0 .OR. para(4).GT.100 )
440  $ info = -(1000 * 4 + 4)
441  IF( para(5).LT.1 .OR. para(5).GT.nb )
442  $ info = -(1000 * 4 + 5)
443  IF( para(6).LT.1 .OR. para(6).GT.para(2) )
444  $ info = -(1000 * 4 + 6)
445  END IF
446 *
447 * Check requirements on IT, JT, IQ and JQ
448 *
449  IF( info.EQ.0 ) THEN
450  IF( it.NE.1 ) info = -7
451  IF( jt.NE.it ) info = -8
452  IF( iq.NE.1 ) info = -11
453  IF( jq.NE.iq ) info = -12
454  END IF
455 *
456 * Test input parameters for global sanity
457 *
458  IF( info.EQ.0 ) THEN
459  CALL pchk1mat( n, 5, n, 5, it, jt, desct, 9, 0, idum1,
460  $ idum2, info )
461  END IF
462  IF( info.EQ.0 ) THEN
463  CALL pchk1mat( n, 5, n, 5, iq, jq, descq, 13, 0, idum1,
464  $ idum2, info )
465  END IF
466  IF( info.EQ.0 ) THEN
467  CALL pchk2mat( n, 5, n, 5, it, jt, desct, 9, n, 5, n, 5,
468  $ iq, jq, descq, 13, 0, idum1, idum2, info )
469  END IF
470 *
471 * Decode and test the input parameters
472 *
473  IF( info.EQ.0 .OR. lquery ) THEN
474  wantbh = lsame( job, 'B' )
475  wants = lsame( job, 'E' ) .OR. wantbh
476  wantsp = lsame( job, 'V' ) .OR. wantbh
477  wantq = lsame( compq, 'V' )
478 *
479  IF( .NOT.lsame( job, 'N' ) .AND. .NOT.wants .AND. .NOT.wantsp )
480  $ THEN
481  info = -1
482  ELSEIF( .NOT.lsame( compq, 'N' ) .AND. .NOT.wantq ) THEN
483  info = -2
484  ELSEIF( n.LT.0 ) THEN
485  info = -4
486  ELSE
487 *
488 * Extract local leading dimension
489 *
490  lldt = desct( lld_ )
491  lldq = descq( lld_ )
492 *
493 * Check the SELECT vector for consistency and set M to the
494 * dimension of the specified invariant subspace.
495 *
496  m = 0
497  DO 10 k = 1, n
498 *
499 * IWORK(1:N) is an integer copy of SELECT.
500 *
501  IF( SELECT(k) ) THEN
502  iwork(k) = 1
503  ELSE
504  iwork(k) = 0
505  END IF
506  IF( k.LT.n ) THEN
507  CALL infog2l( k+1, k, desct, nprow, npcol,
508  $ myrow, mycol, itt, jtt, trsrc, tcsrc )
509  IF( myrow.EQ.trsrc .AND. mycol.EQ.tcsrc ) THEN
510  elem = t( (jtt-1)*lldt + itt )
511  IF( elem.NE.zero ) THEN
512  IF( SELECT(k) .AND. .NOT.SELECT(k+1) ) THEN
513 * INFO = -3
514  iwork(k+1) = 1
515  ELSEIF( .NOT.SELECT(k) .AND. SELECT(k+1) ) THEN
516 * INFO = -3
517  iwork(k) = 1
518  END IF
519  END IF
520  END IF
521  END IF
522  IF( SELECT(k) ) m = m + 1
523  10 CONTINUE
524  mmax = m
525  mmin = m
526  IF( nprocs.GT.1 )
527  $ CALL igamx2d( ictxt, 'All', top, 1, 1, mmax, 1, -1,
528  $ -1, -1, -1, -1 )
529  IF( nprocs.GT.1 )
530  $ CALL igamn2d( ictxt, 'All', top, 1, 1, mmin, 1, -1,
531  $ -1, -1, -1, -1 )
532  IF( mmax.GT.mmin ) THEN
533  m = mmax
534  IF( nprocs.GT.1 )
535  $ CALL igamx2d( ictxt, 'All', top, n, 1, iwork, n,
536  $ -1, -1, -1, -1, -1 )
537  END IF
538 *
539 * Set parameters for deep pipelining in parallel
540 * Sylvester solver.
541 *
542  mbnb2( 1 ) = min( max( para( 3 ), para( 2 )*2 ), nb )
543  mbnb2( 2 ) = mbnb2( 1 )
544 *
545 * Compute needed workspace
546 *
547  n1 = m
548  n2 = n - m
549  IF( wants ) THEN
550 c CALL PGESYCTD( 'Solve', 'Schur', 'Schur', 'Notranspose',
551 c $ 'Notranspose', -1, 'Demand', N1, N2, T, 1, 1, DESCT,
552 c $ T, N1+1, N1+1, DESCT, T, 1, N1+1, DESCT, MBNB2,
553 c $ WORK, -1, IWORK(N+1), -1, NOEXSY, SCALE, IERR )
554  wrk1 = int(work(1))
555  iwrk1 = iwork(n+1)
556  wrk1 = 0
557  iwrk1 = 0
558  ELSE
559  wrk1 = 0
560  iwrk1 = 0
561  END IF
562 *
563  IF( wantsp ) THEN
564 c CALL PSYCTCON( 'Notranspose', 'Notranspose', -1,
565 c $ 'Demand', N1, N2, T, 1, 1, DESCT, T, N1+1, N1+1,
566 c $ DESCT, MBNB2, WORK, -1, IWORK(N+1), -1, EST, ITER,
567 c $ IERR )
568  wrk2 = int(work(1))
569  iwrk2 = iwork(n+1)
570  wrk2 = 0
571  iwrk2 = 0
572  ELSE
573  wrk2 = 0
574  iwrk2 = 0
575  END IF
576 *
577  trows = numroc( n, nb, myrow, desct(rsrc_), nprow )
578  tcols = numroc( n, nb, mycol, desct(csrc_), npcol )
579  wrk3 = n + 7*nb**2 + 2*trows*para( 3 ) + tcols*para( 3 ) +
580  $ max( trows*para( 3 ), tcols*para( 3 ) )
581  iwrk3 = 5*para( 1 ) + para(2)*para(3) -
582  $ para(2) * (para(2) + 1 ) / 2
583 *
584  IF( wantsp ) THEN
585  lwmin = max( 1, max( wrk2, wrk3) )
586  liwmin = max( 1, max( iwrk2, iwrk3 ) )+n
587  ELSE IF( lsame( job, 'N' ) ) THEN
588  lwmin = max( 1, wrk3 )
589  liwmin = iwrk3+n
590  ELSE IF( lsame( job, 'E' ) ) THEN
591  lwmin = max( 1, max( wrk1, wrk3) )
592  liwmin = max( 1, max( iwrk1, iwrk3 ) )+n
593  END IF
594 *
595  IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
596  info = -20
597  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
598  info = -22
599  END IF
600  END IF
601  END IF
602 *
603 * Global maximum on info
604 *
605  IF( nprocs.GT.1 )
606  $ CALL igamx2d( ictxt, 'All', top, 1, 1, info, 1, -1, -1, -1,
607  $ -1, -1 )
608 *
609 * Return if some argument is incorrect
610 *
611  IF( info.NE.0 .AND. .NOT.lquery ) THEN
612  m = 0
613  s = one
614  sep = zero
615  CALL pxerbla( ictxt, 'PDTRSEN', -info )
616  RETURN
617  ELSEIF( lquery ) THEN
618  work( 1 ) = dble(lwmin)
619  iwork( 1 ) = liwmin
620  RETURN
621  END IF
622 *
623 * Quick return if possible.
624 *
625  IF( m.EQ.n .OR. m.EQ.0 ) THEN
626  IF( wants )
627  $ s = one
628  IF( wantsp )
629  $ sep = pdlange( '1', n, n, t, it, jt, desct, work )
630  GO TO 50
631  END IF
632 *
633 * Reorder the eigenvalues.
634 *
635  CALL pdtrord( compq, iwork, para, n, t, it, jt,
636  $ desct, q, iq, jq, descq, wr, wi, m, work, lwork,
637  $ iwork(n+1), liwork-n, info )
638 *
639  IF( wants ) THEN
640 *
641 * Solve Sylvester equation T11*R - R*T2 = scale*T12 for R in
642 * parallel.
643 *
644 * Copy T12 to workspace.
645 *
646  CALL infog2l( 1, n1+1, desct, nprow, npcol, myrow,
647  $ mycol, iloc1, jloc1, trsrc, tcsrc )
648  icofft12 = mod( n1, nb )
649  t12rows = numroc( n1, nb, myrow, trsrc, nprow )
650  t12cols = numroc( n2+icofft12, nb, mycol, tcsrc, npcol )
651  CALL descinit( desct12, n1, n2+icofft12, nb, nb, trsrc,
652  $ tcsrc, ictxt, max(1,t12rows), ierr )
653  CALL pdlacpy( 'All', n1, n2, t, 1, n1+1, desct, work,
654  $ 1, 1+icofft12, desct12 )
655 *
656 * Solve the equation to get the solution in workspace.
657 *
658  space = desct12( lld_ ) * t12cols
659  ipw1 = 1 + space
660 c CALL PGESYCTD( 'Solve', 'Schur', 'Schur', 'Notranspose',
661 c $ 'Notranspose', -1, 'Demand', N1, N2, T, 1, 1, DESCT, T,
662 c $ N1+1, N1+1, DESCT, WORK, 1, 1+ICOFFT12, DESCT12, MBNB2,
663 c $ WORK(IPW1), LWORK-SPACE+1, IWORK(N+1), LIWORK-N, NOEXSY,
664 c $ SCALE, IERR )
665  IF( ierr.LT.0 ) THEN
666  info = n+3
667  ELSE
668  info = n+2
669  END IF
670 *
671 * Estimate the reciprocal of the condition number of the cluster
672 * of eigenvalues.
673 *
674  rnorm = pdlange( 'Frobenius', n1, n2, work, 1, 1+icofft12,
675  $ desct12, dpdum1 )
676  IF( rnorm.EQ.zero ) THEN
677  s = one
678  ELSE
679  s = scale / ( sqrt( scale*scale / rnorm+rnorm )*
680  $ sqrt( rnorm ) )
681  END IF
682  END IF
683 *
684  IF( wantsp ) THEN
685 *
686 * Estimate sep(T11,T21) in parallel.
687 *
688 c CALL PSYCTCON( 'Notranspose', 'Notranspose', -1, 'Demand', N1,
689 c $ N2, T, 1, 1, DESCT, T, N1+1, N1+1, DESCT, MBNB2, WORK,
690 c $ LWORK, IWORK(N+1), LIWORK-N, EST, ITER, IERR )
691  est = est * sqrt(dble(n1*n2))
692  sep = one / est
693  IF( ierr.LT.0 ) THEN
694  info = n+4
695  ELSE
696  info = n+2
697  END IF
698  END IF
699 *
700 * Return to calling program.
701 *
702  50 CONTINUE
703 *
704  RETURN
705 *
706 * End of PDTRSEN
707 *
708  END
709 *
pdtrsen
subroutine pdtrsen(JOB, COMPQ, SELECT, PARA, N, T, IT, JT, DESCT, Q, IQ, JQ, DESCQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdtrsen.f:4
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pchk1mat
subroutine pchk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:3
pchk2mat
subroutine pchk2mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, MB, MBPOS0, NB, NBPOS0, IB, JB, DESCB, DESCBPOS0, NEXTRA, EX, EXPOS, INFO)
Definition: pchkxmat.f:175
pdtrord
subroutine pdtrord(COMPQ, SELECT, PARA, N, T, IT, JT, DESCT, Q, IQ, JQ, DESCQ, WR, WI, M, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pdtrord.f:4
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
chk1mat
subroutine chk1mat(MA, MAPOS0, NA, NAPOS0, IA, JA, DESCA, DESCAPOS0, INFO)
Definition: chk1mat.f:3
pxerbla
subroutine pxerbla(ICTXT, SRNAME, INFO)
Definition: pxerbla.f:2
pdlacpy
subroutine pdlacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pdlacpy.f:3
min
#define min(A, B)
Definition: pcgemr.c:181