ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdlawil.f
Go to the documentation of this file.
1  SUBROUTINE pdlawil( 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 * May 1, 1997
7 *
8 * .. Scalar Arguments ..
9  INTEGER II, JJ, M
10  DOUBLE PRECISION H33, H43H34, H44
11 * ..
12 * .. Array Arguments ..
13  INTEGER DESCA( * )
14  DOUBLE PRECISION A( * ), V( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * PDLAWIL 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) DOUBLE PRECISION 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) DOUBLE PRECISION
102 * These three values are for the double shift QR iteration.
103 * Unchanged on exit.
104 *
105 * V (global output) DOUBLE PRECISION array of size 3.
106 * Contains the transform on ouput.
107 *
108 * Implemented by: G. Henry, November 17, 1996
109 *
110 * =====================================================================
111 *
112 * .. Parameters ..
113  INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
114  $ LLD_, MB_, M_, NB_, N_, RSRC_
115  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
116  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
117  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
118 * ..
119 * .. Local Scalars ..
120  INTEGER CONTXT, DOWN, HBL, ICOL, IROW, JSRC, LDA, LEFT,
121  $ MODKM1, MYCOL, MYROW, NPCOL, NPROW, NUM, RIGHT,
122  $ RSRC, UP
123  DOUBLE PRECISION H11, H12, H21, H22, H33S, H44S, S, V1, V2, V3
124 * ..
125 * .. Local Arrays ..
126  DOUBLE PRECISION BUF( 4 )
127 * ..
128 * .. External Subroutines ..
129  EXTERNAL blacs_gridinfo, dgerv2d, dgesd2d, infog2l
130 * ..
131 * .. Intrinsic Functions ..
132  INTRINSIC abs, mod
133 * ..
134 * .. Executable Statements ..
135 *
136  hbl = desca( mb_ )
137  contxt = desca( ctxt_ )
138  lda = desca( lld_ )
139  CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
140  left = mod( mycol+npcol-1, npcol )
141  right = mod( mycol+1, npcol )
142  up = mod( myrow+nprow-1, nprow )
143  down = mod( myrow+1, nprow )
144  num = nprow*npcol
145 *
146 * On node (II,JJ) collect all DIA,SUP,SUB info from M, M+1
147 *
148  modkm1 = mod( m+1, hbl )
149  IF( modkm1.EQ.0 ) THEN
150  IF( ( myrow.EQ.ii ) .AND. ( right.EQ.jj ) .AND.
151  $ ( npcol.GT.1 ) ) THEN
152  CALL infog2l( m+2, m+1, desca, nprow, npcol, myrow, mycol,
153  $ irow, icol, rsrc, jsrc )
154  buf( 1 ) = a( ( icol-1 )*lda+irow )
155  CALL dgesd2d( contxt, 1, 1, buf, 1, ii, jj )
156  END IF
157  IF( ( down.EQ.ii ) .AND. ( right.EQ.jj ) .AND. ( num.GT.1 ) )
158  $ THEN
159  CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol, irow,
160  $ icol, rsrc, jsrc )
161  buf( 1 ) = a( ( icol-1 )*lda+irow )
162  buf( 2 ) = a( ( icol-1 )*lda+irow+1 )
163  buf( 3 ) = a( icol*lda+irow )
164  buf( 4 ) = a( icol*lda+irow+1 )
165  CALL dgesd2d( contxt, 4, 1, buf, 4, ii, jj )
166  END IF
167  IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
168  CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
169  $ irow, icol, rsrc, jsrc )
170  IF( npcol.GT.1 ) THEN
171  CALL dgerv2d( contxt, 1, 1, v3, 1, myrow, left )
172  ELSE
173  v3 = a( ( icol-2 )*lda+irow )
174  END IF
175  IF( num.GT.1 ) THEN
176  CALL dgerv2d( contxt, 4, 1, buf, 4, up, left )
177  h11 = buf( 1 )
178  h21 = buf( 2 )
179  h12 = buf( 3 )
180  h22 = buf( 4 )
181  ELSE
182  h11 = a( ( icol-3 )*lda+irow-2 )
183  h21 = a( ( icol-3 )*lda+irow-1 )
184  h12 = a( ( icol-2 )*lda+irow-2 )
185  h22 = a( ( icol-2 )*lda+irow-1 )
186  END IF
187  END IF
188  END IF
189  IF( modkm1.EQ.1 ) THEN
190  IF( ( down.EQ.ii ) .AND. ( right.EQ.jj ) .AND. ( num.GT.1 ) )
191  $ THEN
192  CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol, irow,
193  $ icol, rsrc, jsrc )
194  CALL dgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
195  $ jj )
196  END IF
197  IF( ( down.EQ.ii ) .AND. ( mycol.EQ.jj ) .AND. ( nprow.GT.1 ) )
198  $ THEN
199  CALL infog2l( m, m+1, desca, nprow, npcol, myrow, mycol,
200  $ irow, icol, rsrc, jsrc )
201  CALL dgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
202  $ jj )
203  END IF
204  IF( ( myrow.EQ.ii ) .AND. ( right.EQ.jj ) .AND.
205  $ ( npcol.GT.1 ) ) THEN
206  CALL infog2l( m+1, m, desca, nprow, npcol, myrow, mycol,
207  $ irow, icol, rsrc, jsrc )
208  CALL dgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
209  $ jj )
210  END IF
211  IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
212  CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
213  $ irow, icol, rsrc, jsrc )
214  IF( num.GT.1 ) THEN
215  CALL dgerv2d( contxt, 1, 1, h11, 1, up, left )
216  ELSE
217  h11 = a( ( icol-3 )*lda+irow-2 )
218  END IF
219  IF( nprow.GT.1 ) THEN
220  CALL dgerv2d( contxt, 1, 1, h12, 1, up, mycol )
221  ELSE
222  h12 = a( ( icol-2 )*lda+irow-2 )
223  END IF
224  IF( npcol.GT.1 ) THEN
225  CALL dgerv2d( contxt, 1, 1, h21, 1, myrow, left )
226  ELSE
227  h21 = a( ( icol-3 )*lda+irow-1 )
228  END IF
229  h22 = a( ( icol-2 )*lda+irow-1 )
230  v3 = a( ( icol-2 )*lda+irow )
231  END IF
232  END IF
233  IF( ( myrow.NE.ii ) .OR. ( mycol.NE.jj ) )
234  $ RETURN
235 *
236  IF( modkm1.GT.1 ) THEN
237  CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
238  $ irow, icol, rsrc, jsrc )
239  h11 = a( ( icol-3 )*lda+irow-2 )
240  h21 = a( ( icol-3 )*lda+irow-1 )
241  h12 = a( ( icol-2 )*lda+irow-2 )
242  h22 = a( ( icol-2 )*lda+irow-1 )
243  v3 = a( ( icol-2 )*lda+irow )
244  END IF
245 *
246  h44s = h44 - h11
247  h33s = h33 - h11
248  v1 = ( h33s*h44s-h43h34 ) / h21 + h12
249  v2 = h22 - h11 - h33s - h44s
250  s = abs( v1 ) + abs( v2 ) + abs( v3 )
251  v1 = v1 / s
252  v2 = v2 / s
253  v3 = v3 / s
254  v( 1 ) = v1
255  v( 2 ) = v2
256  v( 3 ) = v3
257 *
258  RETURN
259 *
260 * End of PDLAWIL
261 *
262  END
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pdlawil
subroutine pdlawil(II, JJ, M, A, DESCA, H44, H33, H43H34, V)
Definition: pdlawil.f:2