2
3
4
5
6
7
8
9 INTEGER IA, JA, M, N
10
11
12 INTEGER DESCA( * )
13 REAL A( * ), TAU( * ), WORK( * )
14
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
117
118
119
120
121
122 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
123 $ LLD_, MB_, M_, NB_, N_, RSRC_
124 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
125 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
126 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
127 REAL ZERO
128 parameter( zero = 0.0e+0 )
129
130
131 CHARACTER COLBTOP, ROWBTOP
132 INTEGER I, IACOL, IAROW, IB, ICOFF, ICTXT, IIA, IN,
133 $ IPT, IPV, IPW, JJA, JM1, JV, L, MYCOL, MYROW,
134 $ NPCOL, NPROW, NQ
135
136
137 INTEGER DESCV( DLEN_ )
138
139
142 $ pb_topset
143
144
145 INTEGER ICEIL, NUMROC
147
148
150
151
152
153
154
155 ictxt = desca( ctxt_ )
156 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
157
158
159
160 IF( n.LT.m )
161 $ RETURN
162
163 l = n - m
164 jm1 = ja +
min( m+1, n ) - 1
165 in =
min(
iceil( ia, desca( mb_ ) ) * desca( mb_ ), ia+m-1 )
166 icoff = mod( ja-1, desca( nb_ ) )
167 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
168 $ iarow, iacol )
169 nq =
numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
170 ipv = 1
171 ipt = ipv + nq * desca( mb_ )
172 ipw = ipt + desca( mb_ ) * desca( mb_ )
173 CALL pb_topget( ictxt, 'Broadcast', 'Rowwise', rowbtop )
174 CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', colbtop )
175 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', ' ' )
176 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', 'I-ring' )
177
178 CALL descset( descv, desca( mb_ ), n + icoff, desca( mb_ ),
179 $ desca( nb_ ), iarow, iacol, ictxt, desca( mb_ ) )
180
181
182
183 ib = in - ia + 1
184 jv = icoff + jm1 - ja + 1
185
186
187
188 CALL pslarzt(
'Backward',
'Rowwise', l, ib, a, ia, jm1, desca,
189 $ tau, work( ipt ), work( ipw ) )
190
191
192
193 CALL pslacpy(
'All', ib, l, a, ia, jm1, desca, work( ipv ), 1,
194 $ jv, descv )
195
196
197
198 CALL pslacpy(
'Lower', ib-1, ib-1, a, ia+1, ja, desca,
199 $ work( ipv ), 1, icoff+1, descv )
200
201
202
203 CALL pslaset(
'All', ib, l, zero, zero, a, ia, jm1, desca )
204 CALL pslaset(
'Lower', ib-1, ib-1, zero, zero, a, ia+1, ja,
205 $ desca )
206
207
208
209 CALL pslarzb(
'Right',
'Transpose',
'Backward',
'Rowwise',
210 $ in-ia+1, n, ib, l, work( ipv ), 1, jv, descv,
211 $ work( ipt ), a, ia, ja, desca, work( ipw ) )
212
213
214
215 CALL pslacpy(
'Lower', ib-1, ib-1, work( ipv ), 1, icoff+1, descv,
216 $ a, ia+1, ja, desca )
217
218 descv( rsrc_ ) = mod( descv( rsrc_ ) + 1, nprow )
219
220
221
222 DO 10 i = in+1, ia+m-1, desca( mb_ )
223 ib =
min( ia+m-i, desca( mb_ ) )
224
225
226
227 CALL pslarzt(
'Backward',
'Rowwise', l, ib, a, i, jm1, desca,
228 $ tau, work( ipt ), work( ipw ) )
229
230
231
232 CALL pslacpy(
'All', ib, l, a, i, jm1, desca, work( ipv ), 1,
233 $ jv, descv )
234
235
236
237 CALL pslacpy(
'Lower', ib-1, ib-1, a, i+1, ja+i-ia, desca,
238 $ work( ipv ), 1, icoff+1+i-ia, descv )
239
240
241
242 CALL pslaset(
'All', ib, l, zero, zero, a, i, jm1, desca )
243 CALL pslaset(
'Lower', ib-1, ib-1, zero, zero, a, i+1, ja+i-ia,
244 $ desca )
245
246
247
248 CALL pslarzb(
'Right',
'Transpose',
'Backward',
'Rowwise',
249 $ i+ib-ia, n-i+ia, ib, l, work( ipv ), 1, jv,
250 $ descv, work( ipt ), a, ia, ja+i-ia, desca,
251 $ work( ipw ) )
252
253 CALL pslacpy(
'Lower', ib-1, ib-1, work( ipv ), 1,
254 $ icoff+1+i-ia, descv, a, i+1, ja+i-ia, desca )
255
256 descv( rsrc_ ) = mod( descv( rsrc_ ) + 1, nprow )
257
258 10 CONTINUE
259
260 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', rowbtop )
261 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', colbtop )
262
263 RETURN
264
265
266
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
integer function iceil(inum, idenom)
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
integer function numroc(n, nb, iproc, isrcproc, nprocs)
subroutine pslaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
subroutine pslacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
subroutine pslarzb(side, trans, direct, storev, m, n, k, l, v, iv, jv, descv, t, c, ic, jc, descc, work)
subroutine pslarzt(direct, storev, n, k, v, iv, jv, descv, tau, t, work)