ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzlawil.f
Go to the documentation of this file.
1  SUBROUTINE pzlawil( II, JJ, M, A, DESCA, H44, H33, H43H34, V )
2 *
3 * -- ScaLAPACK routine (version 1.7) --
4 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5 * and University of California, Berkeley.
6 * July 31, 2001
7 *
8 * .. Scalar Arguments ..
9  INTEGER II, JJ, M
10  COMPLEX*16 H33, H43H34, H44
11 * ..
12 * .. Array Arguments ..
13  INTEGER DESCA( * )
14  COMPLEX*16 A( * ), V( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * PZLAWIL gets the transform given by H44,H33, & H43H34 into V
21 * starting at row M.
22 *
23 * Notes
24 * =====
25 *
26 * Each global data object is described by an associated description
27 * vector. This vector stores the information required to establish
28 * the mapping between an object element and its corresponding process
29 * and memory location.
30 *
31 * Let A be a generic term for any 2D block cyclicly distributed array.
32 * Such a global array has an associated description vector DESCA.
33 * In the following comments, the character _ should be read as
34 * "of the global array".
35 *
36 * NOTATION STORED IN EXPLANATION
37 * --------------- -------------- --------------------------------------
38 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
39 * DTYPE_A = 1.
40 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
41 * the BLACS process grid A is distribu-
42 * ted over. The context itself is glo-
43 * bal, but the handle (the integer
44 * value) may vary.
45 * M_A (global) DESCA( M_ ) The number of rows in the global
46 * array A.
47 * N_A (global) DESCA( N_ ) The number of columns in the global
48 * array A.
49 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
50 * the rows of the array.
51 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
52 * the columns of the array.
53 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
54 * row of the array A is distributed.
55 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
56 * first column of the array A is
57 * distributed.
58 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
59 * array. LLD_A >= MAX(1,LOCr(M_A)).
60 *
61 * Let K be the number of rows or columns of a distributed matrix,
62 * and assume that its process grid has dimension p x q.
63 * LOCr( K ) denotes the number of elements of K that a process
64 * would receive if K were distributed over the p processes of its
65 * process column.
66 * Similarly, LOCc( K ) denotes the number of elements of K that a
67 * process would receive if K were distributed over the q processes of
68 * its process row.
69 * The values of LOCr() and LOCc() may be determined via a call to the
70 * ScaLAPACK tool function, NUMROC:
71 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
72 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
73 * An upper bound for these quantities may be computed by:
74 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
75 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
76 *
77 * Arguments
78 * =========
79 *
80 * II (global input) INTEGER
81 * Row owner of H(M+2,M+2)
82 *
83 * JJ (global input) INTEGER
84 * Column owner of H(M+2,M+2)
85 *
86 * M (global input) INTEGER
87 * On entry, this is where the transform starts (row M.)
88 * Unchanged on exit.
89 *
90 * A (global input) COMPLEX*16 array, dimension
91 * (DESCA(LLD_),*)
92 * On entry, the Hessenberg matrix.
93 * Unchanged on exit.
94 *
95 * DESCA (global and local input) INTEGER array of dimension DLEN_.
96 * The array descriptor for the distributed matrix A.
97 * Unchanged on exit.
98 *
99 * H44
100 * H33
101 * H43H34 (global input) COMPLEX*16
102 * These three values are for the double shift QR iteration.
103 * Unchanged on exit.
104 *
105 * V (global output) COMPLEX*16 array of size 3.
106 * Contains the transform on ouput.
107 *
108 * Further Details
109 * ===============
110 *
111 * Implemented by: M. Fahey, May 28, 1999
112 *
113 * =====================================================================
114 *
115 * .. Parameters ..
116  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
117  $ LLD_, MB_, M_, NB_, N_, RSRC_
118  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
119  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
120  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
121 * ..
122 * .. Local Scalars ..
123  INTEGER CONTXT, DOWN, HBL, ICOL, IROW, JSRC, LDA, LEFT,
124  $ MODKM1, MYCOL, MYROW, NPCOL, NPROW, NUM, RIGHT,
125  $ RSRC, UP
126  DOUBLE PRECISION S
127  COMPLEX*16 CDUM, H11, H12, H21, H22, H33S, H44S, V1, V2,
128  $ V3
129 * ..
130 * .. Local Arrays ..
131  COMPLEX*16 BUF( 4 )
132 * ..
133 * .. External Subroutines ..
134  EXTERNAL blacs_gridinfo, infog2l, zgerv2d, zgesd2d
135 * ..
136 * .. Intrinsic Functions ..
137  INTRINSIC abs, dble, dimag, mod
138 * ..
139 * .. Statement Functions ..
140  DOUBLE PRECISION CABS1
141 * ..
142 * .. Statement Function definitions ..
143  cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
144 * ..
145 * .. Executable Statements ..
146 *
147  hbl = desca( mb_ )
148  contxt = desca( ctxt_ )
149  lda = desca( lld_ )
150  CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
151  left = mod( mycol+npcol-1, npcol )
152  right = mod( mycol+1, npcol )
153  up = mod( myrow+nprow-1, nprow )
154  down = mod( myrow+1, nprow )
155  num = nprow*npcol
156 *
157 * On node (II,JJ) collect all DIA,SUP,SUB info from M, M+1
158 *
159  modkm1 = mod( m+1, hbl )
160  IF( modkm1.EQ.0 ) THEN
161  IF( ( myrow.EQ.ii ) .AND. ( right.EQ.jj ) .AND.
162  $ ( npcol.GT.1 ) ) THEN
163  CALL infog2l( m+2, m+1, desca, nprow, npcol, myrow, mycol,
164  $ irow, icol, rsrc, jsrc )
165  buf( 1 ) = a( ( icol-1 )*lda+irow )
166  CALL zgesd2d( contxt, 1, 1, buf, 1, ii, jj )
167  END IF
168  IF( ( down.EQ.ii ) .AND. ( right.EQ.jj ) .AND. ( num.GT.1 ) )
169  $ THEN
170  CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol, irow,
171  $ icol, rsrc, jsrc )
172  buf( 1 ) = a( ( icol-1 )*lda+irow )
173  buf( 2 ) = a( ( icol-1 )*lda+irow+1 )
174  buf( 3 ) = a( icol*lda+irow )
175  buf( 4 ) = a( icol*lda+irow+1 )
176  CALL zgesd2d( contxt, 4, 1, buf, 4, ii, jj )
177  END IF
178  IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
179  CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
180  $ irow, icol, rsrc, jsrc )
181  IF( npcol.GT.1 ) THEN
182  CALL zgerv2d( contxt, 1, 1, v3, 1, myrow, left )
183  ELSE
184  v3 = a( ( icol-2 )*lda+irow )
185  END IF
186  IF( num.GT.1 ) THEN
187  CALL zgerv2d( contxt, 4, 1, buf, 4, up, left )
188  h11 = buf( 1 )
189  h21 = buf( 2 )
190  h12 = buf( 3 )
191  h22 = buf( 4 )
192  ELSE
193  h11 = a( ( icol-3 )*lda+irow-2 )
194  h21 = a( ( icol-3 )*lda+irow-1 )
195  h12 = a( ( icol-2 )*lda+irow-2 )
196  h22 = a( ( icol-2 )*lda+irow-1 )
197  END IF
198  END IF
199  END IF
200  IF( modkm1.EQ.1 ) THEN
201  IF( ( down.EQ.ii ) .AND. ( right.EQ.jj ) .AND. ( num.GT.1 ) )
202  $ THEN
203  CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol, irow,
204  $ icol, rsrc, jsrc )
205  CALL zgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
206  $ jj )
207  END IF
208  IF( ( down.EQ.ii ) .AND. ( mycol.EQ.jj ) .AND. ( nprow.GT.1 ) )
209  $ THEN
210  CALL infog2l( m, m+1, desca, nprow, npcol, myrow, mycol,
211  $ irow, icol, rsrc, jsrc )
212  CALL zgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
213  $ jj )
214  END IF
215  IF( ( myrow.EQ.ii ) .AND. ( right.EQ.jj ) .AND.
216  $ ( npcol.GT.1 ) ) THEN
217  CALL infog2l( m+1, m, desca, nprow, npcol, myrow, mycol,
218  $ irow, icol, rsrc, jsrc )
219  CALL zgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
220  $ jj )
221  END IF
222  IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
223  CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
224  $ irow, icol, rsrc, jsrc )
225  IF( num.GT.1 ) THEN
226  CALL zgerv2d( contxt, 1, 1, h11, 1, up, left )
227  ELSE
228  h11 = a( ( icol-3 )*lda+irow-2 )
229  END IF
230  IF( nprow.GT.1 ) THEN
231  CALL zgerv2d( contxt, 1, 1, h12, 1, up, mycol )
232  ELSE
233  h12 = a( ( icol-2 )*lda+irow-2 )
234  END IF
235  IF( npcol.GT.1 ) THEN
236  CALL zgerv2d( contxt, 1, 1, h21, 1, myrow, left )
237  ELSE
238  h21 = a( ( icol-3 )*lda+irow-1 )
239  END IF
240  h22 = a( ( icol-2 )*lda+irow-1 )
241  v3 = a( ( icol-2 )*lda+irow )
242  END IF
243  END IF
244  IF( ( myrow.NE.ii ) .OR. ( mycol.NE.jj ) )
245  $ RETURN
246 *
247  IF( modkm1.GT.1 ) THEN
248  CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
249  $ irow, icol, rsrc, jsrc )
250  h11 = a( ( icol-3 )*lda+irow-2 )
251  h21 = a( ( icol-3 )*lda+irow-1 )
252  h12 = a( ( icol-2 )*lda+irow-2 )
253  h22 = a( ( icol-2 )*lda+irow-1 )
254  v3 = a( ( icol-2 )*lda+irow )
255  END IF
256 *
257  h44s = h44 - h11
258  h33s = h33 - h11
259  v1 = ( h33s*h44s-h43h34 ) / h21 + h12
260  v2 = h22 - h11 - h33s - h44s
261  s = cabs1( v1 ) + cabs1( v2 ) + cabs1( v3 )
262  v1 = v1 / s
263  v2 = v2 / s
264  v3 = v3 / s
265  v( 1 ) = v1
266  v( 2 ) = v2
267  v( 3 ) = v3
268 *
269  RETURN
270 *
271 * End of PZLAWIL
272 *
273  END
pzlawil
subroutine pzlawil(II, JJ, M, A, DESCA, H44, H33, H43H34, V)
Definition: pzlawil.f:2
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3