2
3
4
5
6
7
8
9 INTEGER II, JJ, M
10 REAL H33, H43H34, H44
11
12
13 INTEGER DESCA( * )
14 REAL 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 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
120 INTEGER CONTXT, DOWN, HBL, ICOL, IROW, JSRC, LDA, LEFT,
121 $ MODKM1, MYCOL, MYROW, NPCOL, NPROW, NUM, RIGHT,
122 $ RSRC, UP
123 REAL H22, H33S, H44S, S, V1, V2
124
125
126 REAL BUF( 4 )
127 REAL H11( 1 )
128 REAL H12( 1 )
129 REAL H21( 1 )
130 REAL V3( 1 )
131
132
133 EXTERNAL blacs_gridinfo, sgerv2d, sgesd2d,
infog2l
134
135
136 INTRINSIC abs, mod
137
138
139
140 hbl = desca( mb_ )
141 contxt = desca( ctxt_ )
142 lda = desca( lld_ )
143 CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
144 left = mod( mycol+npcol-1, npcol )
145 right = mod( mycol+1, npcol )
146 up = mod( myrow+nprow-1, nprow )
147 down = mod( myrow+1, nprow )
148 num = nprow*npcol
149
150
151
152 modkm1 = mod( m+1, hbl )
153 IF( modkm1.EQ.0 ) THEN
154 IF( ( myrow.EQ.ii ) .AND. ( right.EQ.jj ) .AND.
155 $ ( npcol.GT.1 ) ) THEN
156 CALL infog2l( m+2, m+1, desca, nprow, npcol, myrow, mycol,
157 $ irow, icol, rsrc, jsrc )
158 buf( 1 ) = a( ( icol-1 )*lda+irow )
159 CALL sgesd2d( contxt, 1, 1, buf, 1, ii, jj )
160 END IF
161 IF( ( down.EQ.ii ) .AND. ( right.EQ.jj ) .AND. ( num.GT.1 ) )
162 $ THEN
163 CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol, irow,
164 $ icol, rsrc, jsrc )
165 buf( 1 ) = a( ( icol-1 )*lda+irow )
166 buf( 2 ) = a( ( icol-1 )*lda+irow+1 )
167 buf( 3 ) = a( icol*lda+irow )
168 buf( 4 ) = a( icol*lda+irow+1 )
169 CALL sgesd2d( contxt, 4, 1, buf, 4, ii, jj )
170 END IF
171 IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
172 CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
173 $ irow, icol, rsrc, jsrc )
174 IF( npcol.GT.1 ) THEN
175 CALL sgerv2d( contxt, 1, 1, v3, 1, myrow, left )
176 ELSE
177 v3( 1 ) = a( ( icol-2 )*lda+irow )
178 END IF
179 IF( num.GT.1 ) THEN
180 CALL sgerv2d( contxt, 4, 1, buf, 4, up, left )
181 h11( 1 ) = buf( 1 )
182 h21( 1 ) = buf( 2 )
183 h12( 1 ) = buf( 3 )
184 h22 = buf( 4 )
185 ELSE
186 h11( 1 ) = a( ( icol-3 )*lda+irow-2 )
187 h21( 1 ) = a( ( icol-3 )*lda+irow-1 )
188 h12( 1 ) = a( ( icol-2 )*lda+irow-2 )
189 h22 = a( ( icol-2 )*lda+irow-1 )
190 END IF
191 END IF
192 END IF
193 IF( modkm1.EQ.1 ) THEN
194 IF( ( down.EQ.ii ) .AND. ( right.EQ.jj ) .AND. ( num.GT.1 ) )
195 $ THEN
196 CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol, irow,
197 $ icol, rsrc, jsrc )
198 CALL sgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
199 $ jj )
200 END IF
201 IF( ( down.EQ.ii ) .AND. ( mycol.EQ.jj ) .AND. ( nprow.GT.1 ) )
202 $ THEN
203 CALL infog2l( m, m+1, desca, nprow, npcol, myrow, mycol,
204 $ irow, icol, rsrc, jsrc )
205 CALL sgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
206 $ jj )
207 END IF
208 IF( ( myrow.EQ.ii ) .AND. ( right.EQ.jj ) .AND.
209 $ ( npcol.GT.1 ) ) THEN
210 CALL infog2l( m+1, m, desca, nprow, npcol, myrow, mycol,
211 $ irow, icol, rsrc, jsrc )
212 CALL sgesd2d( contxt, 1, 1, a( ( icol-1 )*lda+irow ), 1, ii,
213 $ jj )
214 END IF
215 IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
216 CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
217 $ irow, icol, rsrc, jsrc )
218 IF( num.GT.1 ) THEN
219 CALL sgerv2d( contxt, 1, 1, h11, 1, up, left )
220 ELSE
221 h11( 1 ) = a( ( icol-3 )*lda+irow-2 )
222 END IF
223 IF( nprow.GT.1 ) THEN
224 CALL sgerv2d( contxt, 1, 1, h12, 1, up, mycol )
225 ELSE
226 h12( 1 ) = a( ( icol-2 )*lda+irow-2 )
227 END IF
228 IF( npcol.GT.1 ) THEN
229 CALL sgerv2d( contxt, 1, 1, h21, 1, myrow, left )
230 ELSE
231 h21( 1 ) = a( ( icol-3 )*lda+irow-1 )
232 END IF
233 h22 = a( ( icol-2 )*lda+irow-1 )
234 v3( 1 ) = a( ( icol-2 )*lda+irow )
235 END IF
236 END IF
237 IF( ( myrow.NE.ii ) .OR. ( mycol.NE.jj ) )
238 $ RETURN
239
240 IF( modkm1.GT.1 ) THEN
241 CALL infog2l( m+2, m+2, desca, nprow, npcol, myrow, mycol,
242 $ irow, icol, rsrc, jsrc )
243 h11( 1 ) = a( ( icol-3 )*lda+irow-2 )
244 h21( 1 ) = a( ( icol-3 )*lda+irow-1 )
245 h12( 1 ) = a( ( icol-2 )*lda+irow-2 )
246 h22 = a( ( icol-2 )*lda+irow-1 )
247 v3( 1 ) = a( ( icol-2 )*lda+irow )
248 END IF
249
250 h44s = h44 - h11( 1 )
251 h33s = h33 - h11( 1 )
252 v1 = ( h33s*h44s-h43h34 ) / h21( 1 ) + h12( 1 )
253 v2 = h22 - h11( 1 ) - h33s - h44s
254 s = abs( v1 ) + abs( v2 ) + abs( v3( 1 ) )
255 v1 = v1 / s
256 v2 = v2 / s
257 v3( 1 ) = v3( 1 ) / s
258 v( 1 ) = v1
259 v( 2 ) = v2
260 v( 3 ) = v3( 1 )
261
262 RETURN
263
264
265
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)