ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pshseqr.f
Go to the documentation of this file.
1  SUBROUTINE pshseqr( JOB, COMPZ, N, ILO, IHI, H, DESCH, WR, WI, Z,
2  $ DESCZ, WORK, LWORK, IWORK, LIWORK, INFO )
3 *
4 * Contribution from the Department of Computing Science and HPC2N,
5 * Umea University, Sweden
6 *
7 * -- ScaLAPACK driver routine (version 2.0.1) --
8 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
9 * Univ. of Colorado Denver and University of California, Berkeley.
10 * January, 2012
11 *
12  IMPLICIT NONE
13 *
14 * .. Scalar Arguments ..
15  INTEGER IHI, ILO, INFO, LWORK, LIWORK, N
16  CHARACTER COMPZ, JOB
17 * ..
18 * .. Array Arguments ..
19  INTEGER DESCH( * ) , DESCZ( * ), IWORK( * )
20  REAL H( * ), WI( N ), WORK( * ), WR( N ), Z( * )
21 * ..
22 * Purpose
23 * =======
24 *
25 * PSHSEQR computes the eigenvalues of an upper Hessenberg matrix H
26 * and, optionally, the matrices T and Z from the Schur decomposition
27 * H = Z*T*Z**T, where T is an upper quasi-triangular matrix (the
28 * Schur form), and Z is the orthogonal matrix of Schur vectors.
29 *
30 * Optionally Z may be postmultiplied into an input orthogonal
31 * matrix Q so that this routine can give the Schur factorization
32 * of a matrix A which has been reduced to the Hessenberg form H
33 * by the orthogonal matrix Q: A = Q*H*Q**T = (QZ)*T*(QZ)**T.
34 *
35 * Notes
36 * =====
37 *
38 * Each global data object is described by an associated description
39 * vector. This vector stores the information required to establish
40 * the mapping between an object element and its corresponding process
41 * and memory location.
42 *
43 * Let A be a generic term for any 2D block cyclicly distributed array.
44 * Such a global array has an associated description vector DESCA.
45 * In the following comments, the character _ should be read as
46 * "of the global array".
47 *
48 * NOTATION STORED IN EXPLANATION
49 * --------------- -------------- --------------------------------------
50 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
51 * DTYPE_A = 1.
52 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
53 * the BLACS process grid A is distribu-
54 * ted over. The context itself is glo-
55 * bal, but the handle (the integer
56 * value) may vary.
57 * M_A (global) DESCA( M_ ) The number of rows in the global
58 * array A.
59 * N_A (global) DESCA( N_ ) The number of columns in the global
60 * array A.
61 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
62 * the rows of the array.
63 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
64 * the columns of the array.
65 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
66 * row of the array A is distributed.
67 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
68 * first column of the array A is
69 * distributed.
70 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
71 * array. LLD_A >= MAX(1,LOCr(M_A)).
72 *
73 * Let K be the number of rows or columns of a distributed matrix,
74 * and assume that its process grid has dimension p x q.
75 * LOCr( K ) denotes the number of elements of K that a process
76 * would receive if K were distributed over the p processes of its
77 * process column.
78 * Similarly, LOCc( K ) denotes the number of elements of K that a
79 * process would receive if K were distributed over the q processes of
80 * its process row.
81 * The values of LOCr() and LOCc() may be determined via a call to the
82 * ScaLAPACK tool function, NUMROC:
83 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
84 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
85 * An upper bound for these quantities may be computed by:
86 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
87 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
88 *
89 * Arguments
90 * =========
91 *
92 * JOB (global input) CHARACTER*1
93 * = 'E': compute eigenvalues only;
94 * = 'S': compute eigenvalues and the Schur form T.
95 *
96 * COMPZ (global input) CHARACTER*1
97 * = 'N': no Schur vectors are computed;
98 * = 'I': Z is initialized to the unit matrix and the matrix Z
99 * of Schur vectors of H is returned;
100 * = 'V': Z must contain an orthogonal matrix Q on entry, and
101 * the product Q*Z is returned.
102 *
103 * N (global input) INTEGER
104 * The order of the Hessenberg matrix H (and Z if WANTZ).
105 * N >= 0.
106 *
107 * ILO (global input) INTEGER
108 * IHI (global input) INTEGER
109 * It is assumed that H is already upper triangular in rows
110 * and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
111 * set by a previous call to PSGEBAL, and then passed to PSGEHRD
112 * when the matrix output by PSGEBAL is reduced to Hessenberg
113 * form. Otherwise ILO and IHI should be set to 1 and N
114 * respectively. If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
115 * If N = 0, then ILO = 1 and IHI = 0.
116 *
117 * H (global input/output) REAL array, dimension
118 * (DESCH(LLD_),*)
119 * On entry, the upper Hessenberg matrix H.
120 * On exit, if JOB = 'S', H is upper quasi-triangular in
121 * rows and columns ILO:IHI, with 1-by-1 and 2-by-2 blocks on
122 * the main diagonal. The 2-by-2 diagonal blocks (corresponding
123 * to complex conjugate pairs of eigenvalues) are returned in
124 * standard form, with H(i,i) = H(i+1,i+1) and
125 * H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and JOB = 'E', the
126 * contents of H are unspecified on exit.
127 *
128 * DESCH (global and local input) INTEGER array of dimension DLEN_.
129 * The array descriptor for the distributed matrix H.
130 *
131 * WR (global output) REAL array, dimension (N)
132 * WI (global output) REAL array, dimension (N)
133 * The real and imaginary parts, respectively, of the computed
134 * eigenvalues ILO to IHI are stored in the corresponding
135 * elements of WR and WI. If two eigenvalues are computed as a
136 * complex conjugate pair, they are stored in consecutive
137 * elements of WR and WI, say the i-th and (i+1)th, with
138 * WI(i) > 0 and WI(i+1) < 0. If JOB = 'S', the
139 * eigenvalues are stored in the same order as on the diagonal
140 * of the Schur form returned in H.
141 *
142 * Z (global input/output) REAL array.
143 * If COMPZ = 'V', on entry Z must contain the current
144 * matrix Z of accumulated transformations from, e.g., PSGEHRD,
145 * and on exit Z has been updated; transformations are applied
146 * only to the submatrix Z(ILO:IHI,ILO:IHI).
147 * If COMPZ = 'N', Z is not referenced.
148 * If COMPZ = 'I', on entry Z need not be set and on exit,
149 * if INFO = 0, Z contains the orthogonal matrix Z of the Schur
150 * vectors of H.
151 *
152 * DESCZ (global and local input) INTEGER array of dimension DLEN_.
153 * The array descriptor for the distributed matrix Z.
154 *
155 * WORK (local workspace) REAL array, dimension(LWORK)
156 *
157 * LWORK (local input) INTEGER
158 * The length of the workspace array WORK.
159 *
160 * IWORK (local workspace) INTEGER array, dimension (LIWORK)
161 *
162 * LIWORK (local input) INTEGER
163 * The length of the workspace array IWORK.
164 *
165 * INFO (output) INTEGER
166 * = 0: successful exit
167 * .LT. 0: if INFO = -i, the i-th argument had an illegal
168 * value (see also below for -7777 and -8888).
169 * .GT. 0: if INFO = i, PSHSEQR failed to compute all of
170 * the eigenvalues. Elements 1:ilo-1 and i+1:n of WR
171 * and WI contain those eigenvalues which have been
172 * successfully computed. (Failures are rare.)
173 *
174 * If INFO .GT. 0 and JOB = 'E', then on exit, the
175 * remaining unconverged eigenvalues are the eigen-
176 * values of the upper Hessenberg matrix rows and
177 * columns ILO through INFO of the final, output
178 * value of H.
179 *
180 * If INFO .GT. 0 and JOB = 'S', then on exit
181 *
182 * (*) (initial value of H)*U = U*(final value of H)
183 *
184 * where U is an orthogonal matrix. The final
185 * value of H is upper Hessenberg and quasi-triangular
186 * in rows and columns INFO+1 through IHI.
187 *
188 * If INFO .GT. 0 and COMPZ = 'V', then on exit
189 *
190 * (final value of Z) = (initial value of Z)*U
191 *
192 * where U is the orthogonal matrix in (*) (regard-
193 * less of the value of JOB.)
194 *
195 * If INFO .GT. 0 and COMPZ = 'I', then on exit
196 * (final value of Z) = U
197 * where U is the orthogonal matrix in (*) (regard-
198 * less of the value of JOB.)
199 *
200 * If INFO .GT. 0 and COMPZ = 'N', then Z is not
201 * accessed.
202 *
203 * = -7777: PSLAQR0 failed to converge and PSLAQR1 was called
204 * instead. This could happen. Mostly due to a bug.
205 * Please, send a bug report to the authors.
206 * = -8888: PSLAQR1 failed to converge and PSLAQR0 was called
207 * instead. This should not happen.
208 *
209 * ================================================================
210 * Based on contributions by
211 * Robert Granat, Department of Computing Science and HPC2N,
212 * Umea University, Sweden.
213 * ================================================================
214 *
215 * Restrictions: The block size in H and Z must be square and larger
216 * than or equal to six (6) due to restrictions in PSLAQR1, PSLAQR5
217 * and SLAQR6. Moreover, H and Z need to be distributed identically
218 * with the same context.
219 *
220 * ================================================================
221 * References:
222 * K. Braman, R. Byers, and R. Mathias,
223 * The Multi-Shift QR Algorithm Part I: Maintaining Well Focused
224 * Shifts, and Level 3 Performance.
225 * SIAM J. Matrix Anal. Appl., 23(4):929--947, 2002.
226 *
227 * K. Braman, R. Byers, and R. Mathias,
228 * The Multi-Shift QR Algorithm Part II: Aggressive Early
229 * Deflation.
230 * SIAM J. Matrix Anal. Appl., 23(4):948--973, 2002.
231 *
232 * R. Granat, B. Kagstrom, and D. Kressner,
233 * A Novel Parallel QR Algorithm for Hybrid Distributed Momory HPC
234 * Systems.
235 * SIAM J. Sci. Comput., 32(4):2345--2378, 2010.
236 *
237 * ================================================================
238 * .. Parameters ..
239  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
240  $ lld_, mb_, m_, nb_, n_, rsrc_
241  LOGICAL CRSOVER
242  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
243  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
244  $ rsrc_ = 7, csrc_ = 8, lld_ = 9,
245  $ crsover = .true. )
246  INTEGER NTINY
247  parameter( ntiny = 11 )
248  INTEGER NL
249  parameter( nl = 49 )
250  REAL ZERO, ONE
251  parameter( zero = 0.0e0, one = 1.0e0 )
252 * ..
253 * .. Local Scalars ..
254  INTEGER I, KBOT, NMIN, LLDH, LLDZ, ICTXT, NPROW, NPCOL,
255  $ myrow, mycol, hrows, hcols, ipw, nh, nb,
256  $ ii, jj, hrsrc, hcsrc, nprocs, iloc1, jloc1,
257  $ hrsrc1, hcsrc1, k, iloc2, jloc2, iloc3, jloc3,
258  $ iloc4, jloc4, hrsrc2, hcsrc2, hrsrc3, hcsrc3,
259  $ hrsrc4, hcsrc4, liwkopt
260  LOGICAL INITZ, LQUERY, WANTT, WANTZ, PAIR, BORDER
261  REAL TMP1, TMP2, TMP3, TMP4, DUM1, DUM2, DUM3,
262  $ dum4, elem1, elem2, elem3, elem4,
263  $ cs, sn, elem5, tmp, lwkopt
264 * ..
265 * .. Local Arrays ..
266  INTEGER DESCH2( DLEN_ )
267 * ..
268 * .. External Functions ..
269  INTEGER PILAENVX, NUMROC, ICEIL
270  LOGICAL LSAME
271  EXTERNAL pilaenvx, lsame, numroc, iceil
272 * ..
273 * .. External Subroutines ..
274  EXTERNAL pslacpy, pslaqr1, pslaqr0, pslaset, pxerbla
275 * ..
276 * .. Intrinsic Functions ..
277  INTRINSIC float, max, min
278 * ..
279 * .. Executable Statements ..
280 *
281 * Decode and check the input parameters.
282 *
283  info = 0
284  ictxt = desch( ctxt_ )
285  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
286  nprocs = nprow*npcol
287  IF( nprow.EQ.-1 ) info = -(600+ctxt_)
288  IF( info.EQ.0 ) THEN
289  wantt = lsame( job, 'S' )
290  initz = lsame( compz, 'I' )
291  wantz = initz .OR. lsame( compz, 'V' )
292  lldh = desch( lld_ )
293  lldz = descz( lld_ )
294  nb = desch( mb_ )
295  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
296 *
297  IF( .NOT.lsame( job, 'E' ) .AND. .NOT.wantt ) THEN
298  info = -1
299  ELSE IF( .NOT.lsame( compz, 'N' ) .AND. .NOT.wantz ) THEN
300  info = -2
301  ELSE IF( n.LT.0 ) THEN
302  info = -3
303  ELSE IF( ilo.LT.1 .OR. ilo.GT.max( 1, n ) ) THEN
304  info = -4
305  ELSE IF( ihi.LT.min( ilo, n ) .OR. ihi.GT.n ) THEN
306  info = -5
307  ELSEIF( descz( ctxt_ ).NE.desch( ctxt_ ) ) THEN
308  info = -( 1000+ctxt_ )
309  ELSEIF( desch( mb_ ).NE.desch( nb_ ) ) THEN
310  info = -( 700+nb_ )
311  ELSEIF( descz( mb_ ).NE.descz( nb_ ) ) THEN
312  info = -( 1000+nb_ )
313  ELSEIF( desch( mb_ ).NE.descz( mb_ ) ) THEN
314  info = -( 1000+mb_ )
315  ELSEIF( desch( mb_ ).LT.6 ) THEN
316  info = -( 700+nb_ )
317  ELSEIF( descz( mb_ ).LT.6 ) THEN
318  info = -( 1000+mb_ )
319  ELSE
320  CALL chk1mat( n, 3, n, 3, 1, 1, desch, 7, info )
321  IF( info.EQ.0 )
322  $ CALL chk1mat( n, 3, n, 3, 1, 1, descz, 11, info )
323  IF( info.EQ.0 )
324  $ CALL pchk2mat( n, 3, n, 3, 1, 1, desch, 7, n, 3, n, 3,
325  $ 1, 1, descz, 11, 0, iwork, iwork, info )
326  END IF
327  END IF
328 *
329 * Compute required workspace.
330 *
331  CALL pslaqr1( wantt, wantz, n, ilo, ihi, h, desch, wr, wi,
332  $ ilo, ihi, z, descz, work, -1, iwork, -1, info )
333  lwkopt = work(1)
334  liwkopt = iwork(1)
335  CALL pslaqr0( wantt, wantz, n, ilo, ihi, h, desch, wr, wi,
336  $ ilo, ihi, z, descz, work, -1, iwork, -1, info, 0 )
337  IF( n.LT.nl ) THEN
338  hrows = numroc( nl, nb, myrow, desch(rsrc_), nprow )
339  hcols = numroc( nl, nb, mycol, desch(csrc_), npcol )
340  work(1) = work(1) + float(2*hrows*hcols)
341  END IF
342  lwkopt = max( lwkopt, work(1) )
343  liwkopt = max( liwkopt, iwork(1) )
344  work(1) = lwkopt
345  iwork(1) = liwkopt
346 *
347  IF( .NOT.lquery .AND. lwork.LT.int(lwkopt) ) THEN
348  info = -13
349  ELSEIF( .NOT.lquery .AND. liwork.LT.liwkopt ) THEN
350  info = -15
351  END IF
352 *
353  IF( info.NE.0 ) THEN
354 *
355 * Quick return in case of invalid argument.
356 *
357  CALL pxerbla( ictxt, 'PSHSEQR', -info )
358  RETURN
359 *
360  ELSE IF( n.EQ.0 ) THEN
361 *
362 * Quick return in case N = 0; nothing to do.
363 *
364  RETURN
365 *
366  ELSE IF( lquery ) THEN
367 *
368 * Quick return in case of a workspace query.
369 *
370  RETURN
371 *
372  ELSE
373 *
374 * Copy eigenvalues isolated by PSGEBAL.
375 *
376  DO 10 i = 1, ilo - 1
377  CALL infog2l( i, i, desch, nprow, npcol, myrow, mycol, ii,
378  $ jj, hrsrc, hcsrc )
379  IF( myrow.EQ.hrsrc .AND. mycol.EQ.hcsrc ) THEN
380  wr( i ) = h( (jj-1)*lldh + ii )
381  ELSE
382  wr( i ) = zero
383  END IF
384  wi( i ) = zero
385  10 CONTINUE
386  IF( ilo.GT.1 )
387  $ CALL sgsum2d( ictxt, 'All', '1-Tree', ilo-1, 1, wr, n, -1,
388  $ -1 )
389  DO 20 i = ihi + 1, n
390  CALL infog2l( i, i, desch, nprow, npcol, myrow, mycol, ii,
391  $ jj, hrsrc, hcsrc )
392  IF( myrow.EQ.hrsrc .AND. mycol.EQ.hcsrc ) THEN
393  wr( i ) = h( (jj-1)*lldh + ii )
394  ELSE
395  wr( i ) = zero
396  END IF
397  wi( i ) = zero
398  20 CONTINUE
399  IF( ihi.LT.n )
400  $ CALL sgsum2d( ictxt, 'All', '1-Tree', n-ihi, 1, wr(ihi+1),
401  $ n, -1, -1 )
402 *
403 * Initialize Z, if requested.
404 *
405  IF( initz )
406  $ CALL pslaset( 'A', n, n, zero, one, z, 1, 1, descz )
407 *
408 * Quick return if possible.
409 *
410  nprocs = nprow*npcol
411  IF( ilo.EQ.ihi ) THEN
412  CALL infog2l( ilo, ilo, desch, nprow, npcol, myrow,
413  $ mycol, ii, jj, hrsrc, hcsrc )
414  IF( myrow.EQ.hrsrc .AND. mycol.EQ.hcsrc ) THEN
415  wr( ilo ) = h( (jj-1)*lldh + ii )
416  IF( nprocs.GT.1 )
417  $ CALL sgebs2d( ictxt, 'All', '1-Tree', 1, 1, wr(ilo),
418  $ 1 )
419  ELSE
420  CALL sgebr2d( ictxt, 'All', '1-Tree', 1, 1, wr(ilo),
421  $ 1, hrsrc, hcsrc )
422  END IF
423  wi( ilo ) = zero
424  RETURN
425  END IF
426 *
427 * PSLAQR1/PSLAQR0 crossover point.
428 *
429  nh = ihi-ilo+1
430  nmin = pilaenvx( ictxt, 12, 'PSHSEQR',
431  $ job( : 1 ) // compz( : 1 ), n, ilo, ihi, lwork )
432  nmin = max( ntiny, nmin )
433 *
434 * PSLAQR0 for big matrices; PSLAQR1 for small ones.
435 *
436  IF( (.NOT. crsover .AND. nh.GT.ntiny) .OR. nh.GT.nmin .OR.
437  $ desch(rsrc_).NE.0 .OR. desch(csrc_).NE.0 ) THEN
438  CALL pslaqr0( wantt, wantz, n, ilo, ihi, h, desch, wr, wi,
439  $ ilo, ihi, z, descz, work, lwork, iwork, liwork, info,
440  $ 0 )
441  IF( info.GT.0 .AND. ( desch(rsrc_).NE.0 .OR.
442  $ desch(csrc_).NE.0 ) ) THEN
443 *
444 * A rare PSLAQR0 failure! PSLAQR1 sometimes succeeds
445 * when PSLAQR0 fails.
446 *
447  kbot = info
448  CALL pslaqr1( wantt, wantz, n, ilo, ihi, h, desch, wr,
449  $ wi, ilo, ihi, z, descz, work, lwork, iwork,
450  $ liwork, info )
451  info = -7777
452  END IF
453  ELSE
454 *
455 * Small matrix.
456 *
457  CALL pslaqr1( wantt, wantz, n, ilo, ihi, h, desch, wr, wi,
458  $ ilo, ihi, z, descz, work, lwork, iwork, liwork, info )
459 *
460  IF( info.GT.0 ) THEN
461 *
462 * A rare PSLAQR1 failure! PSLAQR0 sometimes succeeds
463 * when PSLAQR1 fails.
464 *
465  kbot = info
466 *
467  IF( n.GE.nl ) THEN
468 *
469 * Larger matrices have enough subdiagonal scratch
470 * space to call PSLAQR0 directly.
471 *
472  CALL pslaqr0( wantt, wantz, n, ilo, kbot, h, desch,
473  $ wr, wi, ilo, ihi, z, descz, work, lwork,
474  $ iwork, liwork, info, 0 )
475  ELSE
476 *
477 * Tiny matrices don't have enough subdiagonal
478 * scratch space to benefit from PSLAQR0. Hence,
479 * tiny matrices must be copied into a larger
480 * array before calling PSLAQR0.
481 *
482  hrows = numroc( nl, nb, myrow, desch(rsrc_), nprow )
483  hcols = numroc( nl, nb, mycol, desch(csrc_), npcol )
484  CALL descinit( desch2, nl, nl, nb, nb, desch(rsrc_),
485  $ desch(csrc_), ictxt, max(1, hrows), info )
486  CALL pslacpy( 'All', n, n, h, 1, 1, desch, work, 1,
487  $ 1, desch2 )
488  CALL pselset( work, n+1, n, desch2, zero )
489  CALL pslaset( 'All', nl, nl-n, zero, zero, work, 1,
490  $ n+1, desch2 )
491  ipw = 1 + desch2(lld_)*hcols
492  CALL pslaqr0( wantt, wantz, nl, ilo, kbot, work,
493  $ desch2, wr, wi, ilo, ihi, z, descz,
494  $ work(ipw), lwork-ipw+1, iwork,
495  $ liwork, info, 0 )
496  IF( wantt .OR. info.NE.0 )
497  $ CALL pslacpy( 'All', n, n, work, 1, 1, desch2,
498  $ h, 1, 1, desch )
499  END IF
500  info = -8888
501  END IF
502  END IF
503 *
504 * Clear out the trash, if necessary.
505 *
506  IF( ( wantt .OR. info.NE.0 ) .AND. n.GT.2 )
507  $ CALL pslaset( 'L', n-2, n-2, zero, zero, h, 3, 1, desch )
508 *
509 * Force any 2-by-2 blocks to be complex conjugate pairs of
510 * eigenvalues by removing false such blocks.
511 *
512  DO 30 i = ilo, ihi-1
513  CALL pselget( 'All', ' ', tmp3, h, i+1, i, desch )
514  IF( tmp3.NE.0.0e+00 ) THEN
515  CALL pselget( 'All', ' ', tmp1, h, i, i, desch )
516  CALL pselget( 'All', ' ', tmp2, h, i, i+1, desch )
517  CALL pselget( 'All', ' ', tmp4, h, i+1, i+1, desch )
518  CALL slanv2( tmp1, tmp2, tmp3, tmp4, dum1, dum2, dum3,
519  $ dum4, cs, sn )
520  IF( tmp3.EQ.0.0e+00 ) THEN
521  IF( wantt ) THEN
522  IF( i+2.LE.n )
523  $ CALL psrot( n-i-1, h, i, i+2, desch,
524  $ desch(m_), h, i+1, i+2, desch, desch(m_),
525  $ cs, sn, work, lwork, info )
526  CALL psrot( i-1, h, 1, i, desch, 1, h, 1, i+1,
527  $ desch, 1, cs, sn, work, lwork, info )
528  END IF
529  IF( wantz ) THEN
530  CALL psrot( n, z, 1, i, descz, 1, z, 1, i+1, descz,
531  $ 1, cs, sn, work, lwork, info )
532  END IF
533  CALL pselset( h, i, i, desch, tmp1 )
534  CALL pselset( h, i, i+1, desch, tmp2 )
535  CALL pselset( h, i+1, i, desch, tmp3 )
536  CALL pselset( h, i+1, i+1, desch, tmp4 )
537  END IF
538  END IF
539  30 CONTINUE
540 *
541 * Read out eigenvalues: first let all the processes compute the
542 * eigenvalue inside their diagonal blocks in parallel, except for
543 * the eigenvalue located next to a block border. After that,
544 * compute all eigenvalues located next to the block borders.
545 * Finally, do a global summation over WR and WI so that all
546 * processors receive the result.
547 *
548  DO 40 k = ilo, ihi
549  wr( k ) = zero
550  wi( k ) = zero
551  40 CONTINUE
552  nb = desch( mb_ )
553 *
554 * Loop 50: extract eigenvalues from the blocks which are not laid
555 * out across a border of the processor mesh, except for those 1x1
556 * blocks on the border.
557 *
558  pair = .false.
559  DO 50 k = ilo, ihi
560  IF( .NOT. pair ) THEN
561  border = mod( k, nb ).EQ.0 .OR. ( k.NE.1 .AND.
562  $ mod( k, nb ).EQ.1 )
563  IF( .NOT. border ) THEN
564  CALL infog2l( k, k, desch, nprow, npcol, myrow,
565  $ mycol, iloc1, jloc1, hrsrc1, hcsrc1 )
566  IF( myrow.EQ.hrsrc1 .AND. mycol.EQ.hcsrc1 ) THEN
567  elem1 = h((jloc1-1)*lldh+iloc1)
568  IF( k.LT.n ) THEN
569  elem3 = h((jloc1-1)*lldh+iloc1+1)
570  ELSE
571  elem3 = zero
572  END IF
573  IF( elem3.NE.zero ) THEN
574  elem2 = h((jloc1)*lldh+iloc1)
575  elem4 = h((jloc1)*lldh+iloc1+1)
576  CALL slanv2( elem1, elem2, elem3, elem4,
577  $ wr( k ), wi( k ), wr( k+1 ), wi( k+1 ),
578  $ sn, cs )
579  pair = .true.
580  ELSE
581  IF( k.GT.1 ) THEN
582  tmp = h((jloc1-2)*lldh+iloc1)
583  IF( tmp.NE.zero ) THEN
584  elem1 = h((jloc1-2)*lldh+iloc1-1)
585  elem2 = h((jloc1-1)*lldh+iloc1-1)
586  elem3 = h((jloc1-2)*lldh+iloc1)
587  elem4 = h((jloc1-1)*lldh+iloc1)
588  CALL slanv2( elem1, elem2, elem3,
589  $ elem4, wr( k-1 ), wi( k-1 ),
590  $ wr( k ), wi( k ), sn, cs )
591  ELSE
592  wr( k ) = elem1
593  END IF
594  ELSE
595  wr( k ) = elem1
596  END IF
597  END IF
598  END IF
599  END IF
600  ELSE
601  pair = .false.
602  END IF
603  50 CONTINUE
604 *
605 * Loop 60: extract eigenvalues from the blocks which are laid
606 * out across a border of the processor mesh. The processors are
607 * numbered as below:
608 *
609 * 1 | 2
610 * --+--
611 * 3 | 4
612 *
613  DO 60 k = iceil(ilo,nb)*nb, ihi-1, nb
614  CALL infog2l( k, k, desch, nprow, npcol, myrow, mycol,
615  $ iloc1, jloc1, hrsrc1, hcsrc1 )
616  CALL infog2l( k, k+1, desch, nprow, npcol, myrow, mycol,
617  $ iloc2, jloc2, hrsrc2, hcsrc2 )
618  CALL infog2l( k+1, k, desch, nprow, npcol, myrow, mycol,
619  $ iloc3, jloc3, hrsrc3, hcsrc3 )
620  CALL infog2l( k+1, k+1, desch, nprow, npcol, myrow, mycol,
621  $ iloc4, jloc4, hrsrc4, hcsrc4 )
622  IF( myrow.EQ.hrsrc2 .AND. mycol.EQ.hcsrc2 ) THEN
623  elem2 = h((jloc2-1)*lldh+iloc2)
624  IF( hrsrc1.NE.hrsrc2 .OR. hcsrc1.NE.hcsrc2 )
625  $ CALL sgesd2d( ictxt, 1, 1, elem2, 1, hrsrc1, hcsrc1)
626  END IF
627  IF( myrow.EQ.hrsrc3 .AND. mycol.EQ.hcsrc3 ) THEN
628  elem3 = h((jloc3-1)*lldh+iloc3)
629  IF( hrsrc1.NE.hrsrc3 .OR. hcsrc1.NE.hcsrc3 )
630  $ CALL sgesd2d( ictxt, 1, 1, elem3, 1, hrsrc1, hcsrc1)
631  END IF
632  IF( myrow.EQ.hrsrc4 .AND. mycol.EQ.hcsrc4 ) THEN
633  work(1) = h((jloc4-1)*lldh+iloc4)
634  IF( k+1.LT.n ) THEN
635  work(2) = h((jloc4-1)*lldh+iloc4+1)
636  ELSE
637  work(2) = zero
638  END IF
639  IF( hrsrc1.NE.hrsrc4 .OR. hcsrc1.NE.hcsrc4 )
640  $ CALL sgesd2d( ictxt, 2, 1, work, 2, hrsrc1, hcsrc1 )
641  END IF
642  IF( myrow.EQ.hrsrc1 .AND. mycol.EQ.hcsrc1 ) THEN
643  elem1 = h((jloc1-1)*lldh+iloc1)
644  IF( hrsrc1.NE.hrsrc2 .OR. hcsrc1.NE.hcsrc2 )
645  $ CALL sgerv2d( ictxt, 1, 1, elem2, 1, hrsrc2, hcsrc2)
646  IF( hrsrc1.NE.hrsrc3 .OR. hcsrc1.NE.hcsrc3 )
647  $ CALL sgerv2d( ictxt, 1, 1, elem3, 1, hrsrc3, hcsrc3)
648  IF( hrsrc1.NE.hrsrc4 .OR. hcsrc1.NE.hcsrc4 )
649  $ CALL sgerv2d( ictxt, 2, 1, work, 2, hrsrc4, hcsrc4 )
650  elem4 = work(1)
651  elem5 = work(2)
652  IF( elem5.EQ.zero ) THEN
653  IF( wr( k ).EQ.zero .AND. wi( k ).EQ.zero ) THEN
654  CALL slanv2( elem1, elem2, elem3, elem4, wr( k ),
655  $ wi( k ), wr( k+1 ), wi( k+1 ), sn, cs )
656  ELSEIF( wr( k+1 ).EQ.zero .AND. wi( k+1 ).EQ.zero )
657  $ THEN
658  wr( k+1 ) = elem4
659  END IF
660  ELSEIF( wr( k ).EQ.zero .AND. wi( k ).EQ.zero )
661  $ THEN
662  wr( k ) = elem1
663  END IF
664  END IF
665  60 CONTINUE
666 *
667  IF( nprocs.GT.1 ) THEN
668  CALL sgsum2d( ictxt, 'All', ' ', ihi-ilo+1, 1, wr(ilo), n,
669  $ -1, -1 )
670  CALL sgsum2d( ictxt, 'All', ' ', ihi-ilo+1, 1, wi(ilo), n,
671  $ -1, -1 )
672  END IF
673 *
674  END IF
675 *
676  work(1) = lwkopt
677  iwork(1) = liwkopt
678  RETURN
679 *
680 * End of PSHSEQR
681 *
682  END
pslaqr1
recursive subroutine pslaqr1(WANTT, WANTZ, N, ILO, IHI, A, DESCA, WR, WI, ILOZ, IHIZ, Z, DESCZ, WORK, LWORK, IWORK, ILWORK, INFO)
Definition: pslaqr1.f:5
max
#define max(A, B)
Definition: pcgemr.c:180
pslaqr0
recursive subroutine pslaqr0(WANTT, WANTZ, N, ILO, IHI, H, DESCH, WR, WI, ILOZ, IHIZ, Z, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO, RECLEVEL)
Definition: pslaqr0.f:4
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pselset
subroutine pselset(A, IA, JA, DESCA, ALPHA)
Definition: pselset.f:2
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
pslacpy
subroutine pslacpy(UPLO, M, N, A, IA, JA, DESCA, B, IB, JB, DESCB)
Definition: pslacpy.f:3
descinit
subroutine descinit(DESC, M, N, MB, NB, IRSRC, ICSRC, ICTXT, LLD, INFO)
Definition: descinit.f:3
psrot
subroutine psrot(N, X, IX, JX, DESCX, INCX, Y, IY, JY, DESCY, INCY, CS, SN, WORK, LWORK, INFO)
Definition: psrot.f:3
pselget
subroutine pselget(SCOPE, TOP, ALPHA, A, IA, JA, DESCA)
Definition: pselget.f:2
pshseqr
subroutine pshseqr(JOB, COMPZ, N, ILO, IHI, H, DESCH, WR, WI, Z, DESCZ, WORK, LWORK, IWORK, LIWORK, INFO)
Definition: pshseqr.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
pslaset
subroutine pslaset(UPLO, M, N, ALPHA, BETA, A, IA, JA, DESCA)
Definition: psblastst.f:6863
min
#define min(A, B)
Definition: pcgemr.c:181