2
3
4
5
6
7
8
9 INTEGER II, JJ, M
10 COMPLEX H33, H43H34, H44
11
12
13 INTEGER DESCA( * )
14 COMPLEX A( * ), V( * )
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
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
123 INTEGER CONTXT, DOWN, HBL, ICOL, IROW, JSRC, LDA, LEFT,
124 $ MODKM1, MYCOL, MYROW, NPCOL, NPROW, NUM, RIGHT,
125 $ RSRC, UP
126 REAL S
127 COMPLEX CDUM, H22, H33S, H44S, V1, V2
128
129
130 COMPLEX BUF( 4 ), V3( 1 ), H11( 1 ), H12( 1 ), H21( 1 )
131
132
133 EXTERNAL blacs_gridinfo,
infog2l, cgerv2d, cgesd2d
134
135
136 INTRINSIC abs, real, aimag, mod
137
138
139 REAL CABS1
140
141
142 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
143
144
145
146 hbl = desca( mb_ )
147 contxt = desca( ctxt_ )
148 lda = desca( lld_ )
149 CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
150 left = mod( mycol+npcol-1, npcol )
151 right = mod( mycol+1, npcol )
152 up = mod( myrow+nprow-1, nprow )
153 down = mod( myrow+1, nprow )
154 num = nprow*npcol
155
156
157
158 modkm1 = mod( m+1, hbl )
159 IF( modkm1.EQ.0 ) THEN
160 IF( ( myrow.EQ.ii ) .AND. ( right.EQ.jj ) .AND.
161 $ ( npcol.GT.1 ) ) THEN
162 CALL infog2l( m+2, m+1, desca, nprow, npcol, myrow, mycol,
163 $ irow, icol, rsrc, jsrc )
164 buf( 1 ) = a( ( icol-1 )*lda+irow )
165 CALL cgesd2d( contxt, 1, 1, buf, 1, ii, jj )
166 END IF
167 IF( ( down.EQ.ii ) .AND. ( right.EQ.jj ) .AND. ( num.GT.1 ) )
168 $ THEN
169 CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol, irow,
170 $ icol, rsrc, jsrc )
171 buf( 1 ) = a( ( icol-1 )*lda+irow )
172 buf( 2 ) = a( ( icol-1 )*lda+irow+1 )
173 buf( 3 ) = a( icol*lda+irow )
174 buf( 4 ) = a( icol*lda+irow+1 )
175 CALL cgesd2d( contxt, 4, 1, buf, 4, ii, jj )
176 END IF
177 IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
178 CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
179 $ irow, icol, rsrc, jsrc )
180 IF( npcol.GT.1 ) THEN
181 CALL cgerv2d( contxt, 1, 1, v3, 1, myrow, left )
182 ELSE
183 v3( 1 ) = a( ( icol-2 )*lda+irow )
184 END IF
185 IF( num.GT.1 ) THEN
186 CALL cgerv2d( contxt, 4, 1, buf, 4, up, left )
187 h11( 1 ) = buf( 1 )
188 h21( 1 ) = buf( 2 )
189 h12( 1 ) = buf( 3 )
190 h22 = buf( 4 )
191 ELSE
192 h11( 1 ) = a( ( icol-3 )*lda+irow-2 )
193 h21( 1 ) = a( ( icol-3 )*lda+irow-1 )
194 h12( 1 ) = a( ( icol-2 )*lda+irow-2 )
195 h22 = a( ( icol-2 )*lda+irow-1 )
196 END IF
197 END IF
198 END IF
199 IF( modkm1.EQ.1 ) THEN
200 IF( ( down.EQ.ii ) .AND. ( right.EQ.jj ) .AND. ( num.GT.1 ) )
201 $ THEN
202 CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol, irow,
203 $ icol, rsrc, jsrc )
204 CALL cgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
205 $ jj )
206 END IF
207 IF( ( down.EQ.ii ) .AND. ( mycol.EQ.jj ) .AND. ( nprow.GT.1 ) )
208 $ THEN
209 CALL infog2l( m, m+1, desca, nprow, npcol, myrow, mycol,
210 $ irow, icol, rsrc, jsrc )
211 CALL cgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
212 $ jj )
213 END IF
214 IF( ( myrow.EQ.ii ) .AND. ( right.EQ.jj ) .AND.
215 $ ( npcol.GT.1 ) ) THEN
216 CALL infog2l( m+1, m, desca, nprow, npcol, myrow, mycol,
217 $ irow, icol, rsrc, jsrc )
218 CALL cgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
219 $ jj )
220 END IF
221 IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
222 CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
223 $ irow, icol, rsrc, jsrc )
224 IF( num.GT.1 ) THEN
225 CALL cgerv2d( contxt, 1, 1, h11( 1 ), 1, up, left )
226 ELSE
227 h11( 1 ) = a( ( icol-3 )*lda+irow-2 )
228 END IF
229 IF( nprow.GT.1 ) THEN
230 CALL cgerv2d( contxt, 1, 1, h12, 1, up, mycol )
231 ELSE
232 h12( 1 ) = a( ( icol-2 )*lda+irow-2 )
233 END IF
234 IF( npcol.GT.1 ) THEN
235 CALL cgerv2d( contxt, 1, 1, h21( 1 ), 1, myrow, left )
236 ELSE
237 h21( 1 ) = a( ( icol-3 )*lda+irow-1 )
238 END IF
239 h22 = a( ( icol-2 )*lda+irow-1 )
240 v3( 1 ) = a( ( icol-2 )*lda+irow )
241 END IF
242 END IF
243 IF( ( myrow.NE.ii ) .OR. ( mycol.NE.jj ) )
244 $ RETURN
245
246 IF( modkm1.GT.1 ) THEN
247 CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
248 $ irow, icol, rsrc, jsrc )
249 h11( 1 ) = a( ( icol-3 )*lda+irow-2 )
250 h21( 1 ) = a( ( icol-3 )*lda+irow-1 )
251 h12( 1 ) = a( ( icol-2 )*lda+irow-2 )
252 h22 = a( ( icol-2 )*lda+irow-1 )
253 v3( 1 ) = a( ( icol-2 )*lda+irow )
254 END IF
255
256 h44s = h44 - h11( 1 )
257 h33s = h33 - h11( 1 )
258 v1 = ( h33s*h44s-h43h34 ) / h21( 1 ) + h12( 1 )
259 v2 = h22 - h11( 1 ) - h33s - h44s
260 s = cabs1( v1 ) + cabs1( v2 ) + cabs1( v3( 1 ) )
261 v1 = v1 / s
262 v2 = v2 / s
263 v3( 1 ) = v3( 1 ) / s
264 v( 1 ) = v1
265 v( 2 ) = v2
266 v( 3 ) = v3( 1 )
267
268 RETURN
269
270
271
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)