OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pdlaqr4.f
Go to the documentation of this file.
1 SUBROUTINE pdlaqr4( WANTT, WANTZ, N, ILO, IHI, A, DESCA, WR, WI,
2 $ ILOZ, IHIZ, Z, DESCZ, T, LDT, V, LDV, WORK,
3 $ LWORK, INFO )
4*
5* Contribution from the Department of Computing Science and HPC2N,
6* Umea University, Sweden
7*
8* -- ScaLAPACK routine (version 2.0.2) --
9* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver
10* May 1 2012
11*
12 IMPLICIT NONE
13*
14* .. Scalar Arguments ..
15 LOGICAL WANTT, WANTZ
16 INTEGER IHI, IHIZ, ILO, ILOZ, INFO, LDT, LDV, LWORK, N
17* ..
18* .. Array Arguments ..
19 INTEGER DESCA( * ), DESCZ( * )
20 DOUBLE PRECISION A( * ), T( LDT, * ), V( LDV, * ), WI( * ),
21 $ work( * ), wr( * ), z( * )
22* ..
23*
24* Purpose
25* =======
26*
27* PDLAQR4 is an auxiliary routine used to find the Schur decomposition
28* and or eigenvalues of a matrix already in Hessenberg form from cols
29* ILO to IHI. This routine requires that the active block is small
30* enough, i.e. IHI-ILO+1 .LE. LDT, so that it can be solved by LAPACK.
31* Normally, it is called by PDLAQR1. All the inputs are assumed to be
32* valid without checking.
33*
34* Notes
35* =====
36*
37* Each global data object is described by an associated description
38* vector. This vector stores the information required to establish
39* the mapping between an object element and its corresponding process
40* and memory location.
41*
42* Let A be a generic term for any 2D block cyclicly distributed array.
43* Such a global array has an associated description vector DESCA.
44* In the following comments, the character _ should be read as
45* "of the global array".
46*
47* NOTATION STORED IN EXPLANATION
48* --------------- -------------- --------------------------------------
49* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
50* DTYPE_A = 1.
51* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
52* the BLACS process grid A is distribu-
53* ted over. The context itself is glo-
54* bal, but the handle (the integer
55* value) may vary.
56* M_A (global) DESCA( M_ ) The number of rows in the global
57* array A.
58* N_A (global) DESCA( N_ ) The number of columns in the global
59* array A.
60* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
61* the rows of the array.
62* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
63* the columns of the array.
64* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
65* row of the array A is distributed.
66* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
67* first column of the array A is
68* distributed.
69* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
70* array. LLD_A >= MAX(1,LOCr(M_A)).
71*
72* Let K be the number of rows or columns of a distributed matrix,
73* and assume that its process grid has dimension p x q.
74* LOCr( K ) denotes the number of elements of K that a process
75* would receive if K were distributed over the p processes of its
76* process column.
77* Similarly, LOCc( K ) denotes the number of elements of K that a
78* process would receive if K were distributed over the q processes of
79* its process row.
80* The values of LOCr() and LOCc() may be determined via a call to the
81* ScaLAPACK tool function, NUMROC:
82* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
83* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
84* An upper bound for these quantities may be computed by:
85* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
86* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
87*
88* Arguments
89* =========
90*
91* WANTT (global input) LOGICAL
92* = .TRUE. : the full Schur form T is required;
93* = .FALSE.: only eigenvalues are required.
94*
95* WANTZ (global input) LOGICAL
96* = .TRUE. : the matrix of Schur vectors Z is required;
97* = .FALSE.: Schur vectors are not required.
98*
99* N (global input) INTEGER
100* The order of the Hessenberg matrix A (and Z if WANTZ).
101* N >= 0.
102*
103* ILO (global input) INTEGER
104* IHI (global input) INTEGER
105* It is assumed that A is already upper quasi-triangular in
106* rows and columns IHI+1:N, and that A(ILO,ILO-1) = 0 (unless
107* ILO = 1). PDLAQR4 works primarily with the Hessenberg
108* submatrix in rows and columns ILO to IHI, but applies
109* transformations to all of H if WANTT is .TRUE..
110* 1 <= ILO <= max(1,IHI); IHI <= N.
111*
112* A (global input/output) DOUBLE PRECISION array, dimension
113* (DESCA(LLD_),*)
114* On entry, the upper Hessenberg matrix A.
115* On exit, if WANTT is .TRUE., A is upper quasi-triangular in
116* rows and columns ILO:IHI, with any 2-by-2 or larger diagonal
117* blocks not yet in standard form. If WANTT is .FALSE., the
118* contents of A are unspecified on exit.
119*
120* DESCA (global and local input) INTEGER array of dimension DLEN_.
121* The array descriptor for the distributed matrix A.
122*
123* WR (global replicated output) DOUBLE PRECISION array,
124* dimension (N)
125* WI (global replicated output) DOUBLE PRECISION array,
126* dimension (N)
127* The real and imaginary parts, respectively, of the computed
128* eigenvalues ILO to IHI are stored in the corresponding
129* elements of WR and WI. If two eigenvalues are computed as a
130* complex conjugate pair, they are stored in consecutive
131* elements of WR and WI, say the i-th and (i+1)th, with
132* WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
133* eigenvalues are stored in the same order as on the diagonal
134* of the Schur form returned in A. A may be returned with
135* larger diagonal blocks until the next release.
136*
137* ILOZ (global input) INTEGER
138* IHIZ (global input) INTEGER
139* Specify the rows of Z to which transformations must be
140* applied if WANTZ is .TRUE..
141* 1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
142*
143* Z (global input/output) DOUBLE PRECISION array.
144* If WANTZ is .TRUE., on entry Z must contain the current
145* matrix Z of transformations accumulated by PDHSEQR, and on
146* exit Z has been updated; transformations are applied only to
147* the submatrix Z(ILOZ:IHIZ,ILO:IHI).
148* If WANTZ is .FALSE., Z is not referenced.
149*
150* DESCZ (global and local input) INTEGER array of dimension DLEN_.
151* The array descriptor for the distributed matrix Z.
152*
153* T (local workspace) DOUBLE PRECISION array, dimension LDT*NW.
154*
155* LDT (local input) INTEGER
156* The leading dimension of the array T.
157* LDT >= IHI-ILO+1.
158*
159* V (local workspace) DOUBLE PRECISION array, dimension LDV*NW.
160*
161* LDV (local input) INTEGER
162* The leading dimension of the array V.
163* LDV >= IHI-ILO+1.
164*
165* WORK (local workspace) DOUBLE PRECISION array, dimension LWORK.
166*
167* LWORK (local input) INTEGER
168* The dimension of the work array WORK.
169* LWORK >= IHI-ILO+1.
170* WORK(LWORK) is a local array and LWORK is assumed big enough.
171* Typically LWORK >= 4*LDS*LDS if this routine is called by
172* PDLAQR1. (LDS = 385, see PDLAQR1)
173*
174* INFO (global output) INTEGER
175* < 0: parameter number -INFO incorrect or inconsistent;
176* = 0: successful exit;
177* > 0: PDLAQR4 failed to compute all the eigenvalues ILO to IHI
178* in a total of 30*(IHI-ILO+1) iterations; if INFO = i,
179* elements i+1:ihi of WR and WI contain those eigenvalues
180* which have been successfully computed.
181*
182* ================================================================
183* Implemented by
184* Meiyue Shao, Department of Computing Science and HPC2N,
185* Umea University, Sweden
186*
187* ================================================================
188* References:
189* B. Kagstrom, D. Kressner, and M. Shao,
190* On Aggressive Early Deflation in Parallel Variants of the QR
191* Algorithm.
192* Para 2010, to appear.
193*
194* ================================================================
195* .. Parameters ..
196 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
197 $ LLD_, MB_, M_, NB_, N_, RSRC_
198 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
199 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
200 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
201 DOUBLE PRECISION ZERO, ONE
202 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
203* ..
204* .. Local Scalars ..
205 INTEGER CONTXT, HBL, I, I1, I2, IAFIRST, ICOL, ICOL1,
206 $ ICOL2, II, IROW, IROW1, IROW2, ITMP1, ITMP2,
207 $ ierr, j, jafirst, jj, k, l, lda, ldz, lldtmp,
208 $ mycol, myrow, node, npcol, nprow, nh, nmin, nz,
209 $ hstep, vstep, kkrow, kkcol, kln, ltop, left,
210 $ right, up, down, d1, d2
211* ..
212* .. Local Arrays ..
213 INTEGER DESCT( 9 ), DESCV( 9 ), DESCWH( 9 ),
214 $ DESCWV( 9 )
215* ..
216* .. External Functions ..
217 INTEGER NUMROC, ILAENV
218 EXTERNAL NUMROC, ILAENV
219* ..
220* .. External Subroutines ..
221 EXTERNAL blacs_gridinfo, infog2l, dlaset,
222 $ dlahqr, dlaqr4, descinit, pdgemm, pdgemr2d,
223 $ dgemm, dlamov, dgesd2d, dgerv2d,
224 $ dgebs2d, dgebr2d, igebs2d, igebr2d
225* ..
226* .. Intrinsic Functions ..
227 INTRINSIC max, min, mod
228* ..
229* .. Executable Statements ..
230*
231 info = 0
232*
233 nh = ihi - ilo + 1
234 nz = ihiz - iloz + 1
235 IF( n.EQ.0 .OR. nh.EQ.0 )
236 $ RETURN
237*
238* NODE (IAFIRST,JAFIRST) OWNS A(1,1)
239*
240 hbl = desca( mb_ )
241 contxt = desca( ctxt_ )
242 lda = desca( lld_ )
243 iafirst = desca( rsrc_ )
244 jafirst = desca( csrc_ )
245 ldz = descz( lld_ )
246 CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
247 node = myrow*npcol + mycol
248 left = mod( mycol+npcol-1, npcol )
249 right = mod( mycol+1, npcol )
250 up = mod( myrow+nprow-1, nprow )
251 down = mod( myrow+1, nprow )
252*
253* I1 and I2 are the indices of the first row and last column of A
254* to which transformations must be applied.
255*
256 i = ihi
257 l = ilo
258 IF( wantt ) THEN
259 i1 = 1
260 i2 = n
261 ltop = 1
262 ELSE
263 i1 = l
264 i2 = i
265 ltop = l
266 END IF
267*
268* Copy the diagonal block to local and call LAPACK.
269*
270 CALL infog2l( ilo, ilo, desca, nprow, npcol, myrow, mycol,
271 $ irow, icol, ii, jj )
272 IF ( myrow .EQ. ii ) THEN
273 CALL descinit( desct, nh, nh, nh, nh, ii, jj, contxt,
274 $ ldt, ierr )
275 CALL descinit( descv, nh, nh, nh, nh, ii, jj, contxt,
276 $ ldv, ierr )
277 ELSE
278 CALL descinit( desct, nh, nh, nh, nh, ii, jj, contxt,
279 $ 1, ierr )
280 CALL descinit( descv, nh, nh, nh, nh, ii, jj, contxt,
281 $ 1, ierr )
282 END IF
283 CALL pdgemr2d( nh, nh, a, ilo, ilo, desca, t, 1, 1, desct,
284 $ contxt )
285 IF ( myrow .EQ. ii .AND. mycol .EQ. jj ) THEN
286 CALL dlaset( 'All', nh, nh, zero, one, v, ldv )
287 nmin = ilaenv( 12, 'DLAQR3', 'SV', nh, 1, nh, lwork )
288 IF( nh .GT. nmin ) THEN
289 CALL dlaqr4( .true., .true., nh, 1, nh, t, ldt, wr( ilo ),
290 $ wi( ilo ), 1, nh, v, ldv, work, lwork, info )
291* Clean up the scratch used by DLAQR4.
292 CALL dlaset( 'l', NH-2, NH-2, ZERO, ZERO, T( 3, 1 ), LDT )
293 ELSE
294 CALL DLAHQR( .TRUE., .TRUE., NH, 1, NH, T, LDT, WR( ILO ),
295 $ WI( ILO ), 1, NH, V, LDV, INFO )
296 END IF
297 CALL DGEBS2D( CONTXT, 'all', ' ', NH, NH, V, LDV )
298 CALL IGEBS2D( CONTXT, 'all', ' ', 1, 1, INFO, 1 )
299 ELSE
300 CALL DGEBR2D( CONTXT, 'all', ' ', NH, NH, V, LDV, II, JJ )
301 CALL IGEBR2D( CONTXT, 'all', ' ', 1, 1, INFO, 1, II, JJ )
302 END IF
303.NE. IF( INFO 0 ) INFO = INFO+ILO-1
304*
305* Copy the local matrix back to the diagonal block.
306*
307 CALL PDGEMR2D( NH, NH, T, 1, 1, DESCT, A, ILO, ILO, DESCA,
308 $ CONTXT )
309*
310* Update T and Z.
311*
312.LE. IF( MOD( ILO-1, HBL )+NH HBL ) THEN
313*
314* Simplest case: the diagonal block is located on one processor.
315* Call DGEMM directly to perform the update.
316*
317 HSTEP = LWORK / NH
318 VSTEP = HSTEP
319*
320 IF( WANTT ) THEN
321*
322* Update horizontal slab in A.
323*
324 CALL INFOG2L( ILO, I+1, DESCA, NPROW, NPCOL, MYROW,
325 $ MYCOL, IROW, ICOL, II, JJ )
326.EQ. IF( MYROW II ) THEN
327 ICOL1 = NUMROC( N, HBL, MYCOL, JAFIRST, NPCOL )
328 DO 10 KKCOL = ICOL, ICOL1, HSTEP
329 KLN = MIN( HSTEP, ICOL1-KKCOL+1 )
330 CALL DGEMM( 't', 'n', NH, KLN, NH, ONE, V,
331 $ LDV, A( IROW+(KKCOL-1)*LDA ), LDA, ZERO, WORK,
332 $ NH )
333 CALL DLAMOV( 'a', NH, KLN, WORK, NH,
334 $ A( IROW+(KKCOL-1)*LDA ), LDA )
335 10 CONTINUE
336 END IF
337*
338* Update vertical slab in A.
339*
340 CALL INFOG2L( LTOP, ILO, DESCA, NPROW, NPCOL, MYROW,
341 $ MYCOL, IROW, ICOL, II, JJ )
342.EQ. IF( MYCOL JJ ) THEN
343 CALL INFOG2L( ILO-1, ILO, DESCA, NPROW, NPCOL,
344 $ MYROW, MYCOL, IROW1, ICOL1, ITMP1, ITMP2 )
345.NE. IF( MYROW ITMP1 ) IROW1 = IROW1-1
346 DO 20 KKROW = IROW, IROW1, VSTEP
347 KLN = MIN( VSTEP, IROW1-KKROW+1 )
348 CALL DGEMM( 'n', 'n', KLN, NH, NH, ONE,
349 $ A( KKROW+(ICOL-1)*LDA ), LDA, V, LDV, ZERO,
350 $ WORK, KLN )
351 CALL DLAMOV( 'a', KLN, NH, WORK, KLN,
352 $ A( KKROW+(ICOL-1)*LDA ), LDA )
353 20 CONTINUE
354 END IF
355 END IF
356*
357* Update vertical slab in Z.
358*
359 IF( WANTZ ) THEN
360 CALL INFOG2L( ILOZ, ILO, DESCZ, NPROW, NPCOL, MYROW,
361 $ MYCOL, IROW, ICOL, II, JJ )
362.EQ. IF( MYCOL JJ ) THEN
363 CALL INFOG2L( IHIZ, ILO, DESCZ, NPROW, NPCOL,
364 $ MYROW, MYCOL, IROW1, ICOL1, ITMP1, ITMP2 )
365.NE. IF( MYROW ITMP1 ) IROW1 = IROW1-1
366 DO 30 KKROW = IROW, IROW1, VSTEP
367 KLN = MIN( VSTEP, IROW1-KKROW+1 )
368 CALL DGEMM( 'n', 'n', KLN, NH, NH, ONE,
369 $ Z( KKROW+(ICOL-1)*LDZ ), LDZ, V, LDV, ZERO,
370 $ WORK, KLN )
371 CALL DLAMOV( 'a', KLN, NH, WORK, KLN,
372 $ Z( KKROW+(ICOL-1)*LDZ ), LDZ )
373 30 CONTINUE
374 END IF
375 END IF
376*
377.LE. ELSE IF( MOD( ILO-1, HBL )+NH 2*HBL ) THEN
378*
379* More complicated case: the diagonal block lay on a 2x2
380* processor mesh.
381* Call DGEMM locally and communicate by pair.
382*
383 D1 = HBL - MOD( ILO-1, HBL )
384 D2 = NH - D1
385 HSTEP = LWORK / NH
386 VSTEP = HSTEP
387*
388 IF( WANTT ) THEN
389*
390* Update horizontal slab in A.
391*
392 CALL INFOG2L( ILO, I+1, DESCA, NPROW, NPCOL, MYROW,
393 $ MYCOL, IROW, ICOL, II, JJ )
394.EQ. IF( MYROW UP ) THEN
395.EQ. IF( MYROW II ) THEN
396 ICOL1 = NUMROC( N, HBL, MYCOL, JAFIRST, NPCOL )
397 DO 40 KKCOL = ICOL, ICOL1, HSTEP
398 KLN = MIN( HSTEP, ICOL1-KKCOL+1 )
399 CALL DGEMM( 't', 'n', NH, KLN, NH, ONE, V,
400 $ NH, A( IROW+(KKCOL-1)*LDA ), LDA, ZERO,
401 $ WORK, NH )
402 CALL DLAMOV( 'a', NH, KLN, WORK, NH,
403 $ A( IROW+(KKCOL-1)*LDA ), LDA )
404 40 CONTINUE
405 END IF
406 ELSE
407.EQ. IF( MYROW II ) THEN
408 ICOL1 = NUMROC( N, HBL, MYCOL, JAFIRST, NPCOL )
409 DO 50 KKCOL = ICOL, ICOL1, HSTEP
410 KLN = MIN( HSTEP, ICOL1-KKCOL+1 )
411 CALL DGEMM( 't', 'n', D2, KLN, D1, ONE,
412 $ V( 1, D1+1 ), LDV, A( IROW+(KKCOL-1)*LDA ),
413 $ LDA, ZERO, WORK( D1+1 ), NH )
414 CALL DGESD2D( CONTXT, D2, KLN, WORK( D1+1 ),
415 $ NH, DOWN, MYCOL )
416 CALL DGERV2D( CONTXT, D1, KLN, WORK, NH, DOWN,
417 $ MYCOL )
418 CALL DGEMM( 't', 'N', d1, kln, d1, one,
419 $ v, ldv, a( irow+(kkcol-1)*lda ), lda, one,
420 $ work, nh )
421 CALL dlamov( 'A', d1, kln, work, nh,
422 $ a( irow+(kkcol-1)*lda ), lda )
423 50 CONTINUE
424 ELSE IF( up .EQ. ii ) THEN
425 icol1 = numroc( n, hbl, mycol, jafirst, npcol )
426 DO 60 kkcol = icol, icol1, hstep
427 kln = min( hstep, icol1-kkcol+1 )
428 CALL dgemm( 'T', 'N', d1, kln, d2, one,
429 $ v( d1+1, 1 ), ldv, a( irow+(kkcol-1)*lda ),
430 $ lda, zero, work, nh )
431 CALL dgesd2d( contxt, d1, kln, work, nh, up,
432 $ mycol )
433 CALL dgerv2d( contxt, d2, kln, work( d1+1 ),
434 $ nh, up, mycol )
435 CALL dgemm( 'T', 'N', d2, kln, d2, one,
436 $ v( d1+1, d1+1 ), ldv,
437 $ a( irow+(kkcol-1)*lda ), lda, one,
438 $ work( d1+1 ), nh )
439 CALL dlamov( 'A', d2, kln, work( d1+1 ), nh,
440 $ a( irow+(kkcol-1)*lda ), lda )
441 60 CONTINUE
442 END IF
443 END IF
444*
445* Update vertical slab in A.
446*
447 CALL infog2l( ltop, ilo, desca, nprow, npcol, myrow,
448 $ mycol, irow, icol, ii, jj )
449 IF( mycol .EQ. left ) THEN
450 IF( mycol .EQ. jj ) THEN
451 CALL infog2l( ilo-1, ilo, desca, nprow, npcol,
452 $ myrow, mycol, irow1, icol1, itmp1, itmp2 )
453 IF( myrow .NE. itmp1 ) irow1 = irow1-1
454 DO 70 kkrow = irow, irow1, vstep
455 kln = min( vstep, irow1-kkrow+1 )
456 CALL dgemm( 'N', 'N', kln, nh, nh, one,
457 $ a( kkrow+(icol-1)*lda ), lda, v, ldv,
458 $ zero, work, kln )
459 CALL dlamov( 'A', kln, nh, work, kln,
460 $ a( kkrow+(icol-1)*lda ), lda )
461 70 CONTINUE
462 END IF
463 ELSE
464 IF( mycol .EQ. jj ) THEN
465 CALL infog2l( ilo-1, ilo, desca, nprow, npcol,
466 $ myrow, mycol, irow1, icol1, itmp1, itmp2 )
467 IF( myrow .NE. itmp1 ) irow1 = irow1-1
468 DO 80 kkrow = irow, irow1, vstep
469 kln = min( vstep, irow1-kkrow+1 )
470 CALL dgemm( 'N', 'N', kln, d2, d1, one,
471 $ a( kkrow+(icol-1)*lda ), lda, v( 1, d1+1 ),
472 $ ldv, zero, work( 1+d1*kln ), kln )
473 CALL dgesd2d( contxt, kln, d2, work( 1+d1*kln ),
474 $ kln, myrow, right )
475 CALL dgerv2d( contxt, kln, d1, work, kln, myrow,
476 $ right )
477 CALL dgemm( 'N', 'N', kln, d1, d1, one,
478 $ a( kkrow+(icol-1)*lda ), lda, v, ldv, one,
479 $ work, kln )
480 CALL dlamov( 'A', kln, d1, work, kln,
481 $ a( kkrow+(icol-1)*lda ), lda )
482 80 CONTINUE
483 ELSE IF ( left .EQ. jj ) THEN
484 CALL infog2l( ilo-1, ilo, desca, nprow, npcol,
485 $ myrow, mycol, irow1, icol1, itmp1, itmp2 )
486 IF( myrow .NE. itmp1 ) irow1 = irow1-1
487 DO 90 kkrow = irow, irow1, vstep
488 kln = min( vstep, irow1-kkrow+1 )
489 CALL dgemm( 'N', 'N', kln, d1, d2, one,
490 $ a( kkrow+(icol-1)*lda ), lda, v( d1+1, 1 ),
491 $ ldv, zero, work, kln )
492 CALL dgesd2d( contxt, kln, d1, work, kln, myrow,
493 $ left )
494 CALL dgerv2d( contxt, kln, d2, work( 1+d1*kln ),
495 $ kln, myrow, left )
496 CALL dgemm( 'N', 'N', kln, d2, d2, one,
497 $ a( kkrow+(icol-1)*lda ), lda, v( d1+1, d1+1 ),
498 $ ldv, one, work( 1+d1*kln ), kln )
499 CALL dlamov( 'A', kln, d2, work( 1+d1*kln ), kln,
500 $ a( kkrow+(icol-1)*lda ), lda )
501 90 CONTINUE
502 END IF
503 END IF
504 END IF
505*
506* Update vertical slab in Z.
507*
508 IF( wantz ) THEN
509 CALL infog2l( iloz, ilo, descz, nprow, npcol, myrow,
510 $ mycol, irow, icol, ii, jj )
511 IF( mycol .EQ. left ) THEN
512 IF( mycol .EQ. jj ) THEN
513 CALL infog2l( ihiz, ilo, descz, nprow, npcol,
514 $ myrow, mycol, irow1, icol1, itmp1, itmp2 )
515 IF( myrow .NE. itmp1 ) irow1 = irow1-1
516 DO 100 kkrow = irow, irow1, vstep
517 kln = min( vstep, irow1-kkrow+1 )
518 CALL dgemm( 'N', 'N', kln, nh, nh, one,
519 $ z( kkrow+(icol-1)*ldz ), ldz, v, ldv, zero,
520 $ work, kln )
521 CALL dlamov( 'A', kln, nh, work, kln,
522 $ z( kkrow+(icol-1)*ldz ), ldz )
523 100 CONTINUE
524 END IF
525 ELSE
526 IF( mycol .EQ. jj ) THEN
527 CALL infog2l( ihiz, ilo, descz, nprow, npcol,
528 $ myrow, mycol, irow1, icol1, itmp1, itmp2 )
529 IF( myrow .NE. itmp1 ) irow1 = irow1-1
530 DO 110 kkrow = irow, irow1, vstep
531 kln = min( vstep, irow1-kkrow+1 )
532 CALL dgemm( 'N', 'N', kln, d2, d1, one,
533 $ z( kkrow+(icol-1)*ldz ), ldz, v( 1, d1+1 ),
534 $ ldv, zero, work( 1+d1*kln ), kln )
535 CALL dgesd2d( contxt, kln, d2, work( 1+d1*kln ),
536 $ kln, myrow, right )
537 CALL dgerv2d( contxt, kln, d1, work, kln, myrow,
538 $ right )
539 CALL dgemm( 'N', 'N', kln, d1, d1, one,
540 $ z( kkrow+(icol-1)*ldz ), ldz, v, ldv, one,
541 $ work, kln )
542 CALL dlamov( 'A', kln, d1, work, kln,
543 $ z( kkrow+(icol-1)*ldz ), ldz )
544 110 CONTINUE
545 ELSE IF( left .EQ. jj ) THEN
546 CALL infog2l( ihiz, ilo, descz, nprow, npcol,
547 $ myrow, mycol, irow1, icol1, itmp1, itmp2 )
548 IF( myrow .NE. itmp1 ) irow1 = irow1-1
549 DO 120 kkrow = irow, irow1, vstep
550 kln = min( vstep, irow1-kkrow+1 )
551 CALL dgemm( 'N', 'N', kln, d1, d2, one,
552 $ z( kkrow+(icol-1)*ldz ), ldz, v( d1+1, 1 ),
553 $ ldv, zero, work, kln )
554 CALL dgesd2d( contxt, kln, d1, work, kln, myrow,
555 $ left )
556 CALL dgerv2d( contxt, kln, d2, work( 1+d1*kln ),
557 $ kln, myrow, left )
558 CALL dgemm( 'N', 'N', kln, d2, d2, one,
559 $ z( kkrow+(icol-1)*ldz ), ldz,
560 $ v( d1+1, d1+1 ), ldv, one, work( 1+d1*kln ),
561 $ kln )
562 CALL dlamov( 'A', kln, d2, work( 1+d1*kln ),
563 $ kln, z( kkrow+(icol-1)*ldz ), ldz )
564 120 CONTINUE
565 END IF
566 END IF
567 END IF
568*
569 ELSE
570*
571* Most complicated case: the diagonal block lay across the border
572* of the processor mesh.
573* Treat V as a distributed matrix and call PDGEMM.
574*
575 hstep = lwork / nh * npcol
576 vstep = lwork / nh * nprow
577 lldtmp = numroc( nh, nh, myrow, 0, nprow )
578 lldtmp = max( 1, lldtmp )
579 CALL descinit( descv, nh, nh, nh, nh, 0, 0, contxt,
580 $ lldtmp, ierr )
581 CALL descinit( descwh, nh, hstep, nh, lwork / nh, 0, 0,
582 $ contxt, lldtmp, ierr )
583*
584 IF( wantt ) THEN
585*
586* Update horizontal slab in A.
587*
588 DO 130 kkcol = i+1, n, hstep
589 kln = min( hstep, n-kkcol+1 )
590 CALL pdgemm( 'T', 'N', nh, kln, nh, one, v, 1, 1,
591 $ descv, a, ilo, kkcol, desca, zero, work, 1, 1,
592 $ descwh )
593 CALL pdgemr2d( nh, kln, work, 1, 1, descwh, a,
594 $ ilo, kkcol, desca, contxt )
595 130 CONTINUE
596*
597* Update vertical slab in A.
598*
599 DO 140 kkrow = ltop, ilo-1, vstep
600 kln = min( vstep, ilo-kkrow )
601 lldtmp = numroc( kln, lwork / nh, myrow, 0, nprow )
602 lldtmp = max( 1, lldtmp )
603 CALL descinit( descwv, kln, nh, lwork / nh, nh, 0, 0,
604 $ contxt, lldtmp, ierr )
605 CALL pdgemm( 'N', 'N', kln, nh, nh, one, a, kkrow,
606 $ ilo, desca, v, 1, 1, descv, zero, work, 1, 1,
607 $ descwv )
608 CALL pdgemr2d( kln, nh, work, 1, 1, descwv, a, kkrow,
609 $ ilo, desca, contxt )
610 140 CONTINUE
611 END IF
612*
613* Update vertical slab in Z.
614*
615 IF( wantz ) THEN
616 DO 150 kkrow = iloz, ihiz, vstep
617 kln = min( vstep, ihiz-kkrow+1 )
618 lldtmp = numroc( kln, lwork / nh, myrow, 0, nprow )
619 lldtmp = max( 1, lldtmp )
620 CALL descinit( descwv, kln, nh, lwork / nh, nh, 0, 0,
621 $ contxt, lldtmp, ierr )
622 CALL pdgemm( 'N', 'N', kln, nh, nh, one, z, kkrow,
623 $ ilo, descz, v, 1, 1, descv, zero, work, 1, 1,
624 $ descwv )
625 CALL pdgemr2d( kln, nh, work, 1, 1, descwv, z,
626 $ kkrow, ilo, descz, contxt )
627 150 CONTINUE
628 END IF
629 END IF
630*
631* END OF PDLAQR4
632*
633 END
dlaset
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
dlahqr
subroutine dlahqr(wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, iloz, ihiz, z, ldz, info)
DLAHQR computes the eigenvalues and Schur factorization of an upper Hessenberg matrix,...
Definition dlahqr.f:207
dlaqr4
subroutine dlaqr4(wantt, wantz, n, ilo, ihi, h, ldh, wr, wi, iloz, ihiz, z, ldz, work, lwork, info)
DLAQR4 computes the eigenvalues of a Hessenberg matrix, and optionally the matrices from the Schur de...
Definition dlaqr4.f:263
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:187
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dgebs2d
subroutine dgebs2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1082
dgebr2d
subroutine dgebr2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1123
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition mpi.f:937
descinit
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
Definition mpi.f:777
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
pdlaqr4
subroutine pdlaqr4(wantt, wantz, n, ilo, ihi, a, desca, wr, wi, iloz, ihiz, z, descz, t, ldt, v, ldv, work, lwork, info)
Definition pdlaqr4.f:4