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 ONE, ZERO
128 parameter( one = 1.0e+0, 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, JV, K, MYCOL, MYROW, NPCOL,
134 $ 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
159 in =
min(
iceil( ia+m-k, desca( mb_ ) ) * desca( mb_ ), ia+m-1 )
160
161 icoff = mod( ja-1, desca( nb_ ) )
162 CALL infog2l( ia+m-k, ja, desca, nprow, npcol, myrow, mycol,
163 $ iia, jja, iarow, iacol )
164 nq =
numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
165 ipv = 1
166 ipt = ipv + nq * desca( mb_ )
167 ipw = ipt + desca( mb_ ) * desca( mb_ )
168 CALL pb_topget( ictxt, 'Broadcast', 'Rowwise', rowbtop )
169 CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', colbtop )
170 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', ' ' )
171 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', 'I-ring' )
172
173 CALL descset( descv, desca( mb_), n + icoff, desca( mb_ ),
174 $ desca( nb_ ), iarow, iacol, ictxt, desca( mb_ ) )
175
176
177
178 ib = in - ia - m + k + 1
179 jv = 1 + n - k + icoff
180
181
182
183 CALL pslarft(
'Backward',
'Rowwise', n-m+in-ia+1, ib, a, ia+m-k,
184 $ ja, desca, tau, work( ipt ), work( ipw ) )
185
186
187
188 CALL pslacpy(
'All', ib, n-m+in-ia+1, a, ia+m-k, ja, desca,
189 $ work( ipv ), 1, icoff+1, descv )
190 CALL pslaset(
'Upper', ib, ib, zero, one, work( ipv ), 1, jv,
191 $ descv )
192
193
194
195
196 CALL pslaset(
'All', ib, n-k, zero, zero, a, ia+m-k, ja,
197 $ desca )
198 CALL pslaset(
'Lower', ib-1, ib, zero, zero, a, ia+m-k+1,
199 $ ja+n-k, desca )
200
201
202
203 CALL pslarfb(
'Right',
'Transpose',
'Backward',
'Rowwise',
204 $ in-ia+1, n-m+in-ia+1, ib, work( ipv ), 1, icoff+1,
205 $ descv, work( ipt ), a, ia, ja, desca, work( ipw ) )
206
207 descv( rsrc_ ) = mod( descv( rsrc_ ) + 1, nprow )
208
209
210
211 DO 10 i = in+1, ia+m-1, desca( mb_ )
212 ib =
min( ia+m-i, desca( mb_ ) )
213 jv = 1 + n - m + i - ia + icoff
214
215
216
217 CALL pslarft(
'Backward',
'Rowwise', n-m+i+ib-ia, ib, a, i, ja,
218 $ desca, tau, work( ipt ), work( ipw ) )
219
220
221
222 CALL pslacpy(
'All', ib, n-m+i+ib-ia, a, i, ja, desca,
223 $ work( ipv ), 1, icoff+1, descv )
224 CALL pslaset(
'Upper', ib, ib, zero, one, work( ipv ), 1, jv,
225 $ descv )
226
227
228
229
230 CALL pslaset(
'All', ib, n-m+i-ia, zero, zero, a, i, ja,
231 $ desca )
232 CALL pslaset(
'Lower', ib-1, ib, zero, zero, a, i+1,
233 $ ja+n-m+i-ia, desca )
234
235
236
237 CALL pslarfb(
'Right',
'Transpose',
'Backward',
'Rowwise',
238 $ i+ib-ia, n-m+i+ib-ia, ib, work( ipv ), 1,
239 $ icoff+1, descv, work( ipt ), a, ia, ja, desca,
240 $ work( ipw ) )
241
242 descv( rsrc_ ) = mod( descv( rsrc_ ) + 1, nprow )
243
244 10 CONTINUE
245
246 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', rowbtop )
247 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', colbtop )
248
249 RETURN
250
251
252
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 pslarfb(side, trans, direct, storev, m, n, k, v, iv, jv, descv, t, c, ic, jc, descc, work)
subroutine pslarft(direct, storev, n, k, v, iv, jv, descv, tau, t, work)