2
3
4
5
6
7
8
9 INTEGER IA, JA, M, N
10
11
12 INTEGER DESCA( * )
13 COMPLEX 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 COMPLEX ONE, ZERO
128 parameter( one = ( 1.0e+0, 0.0e+0 ),
129 $ zero = ( 0.0e+0, 0.0e+0 ) )
130
131
132 CHARACTER COLBTOP, ROWBTOP
133 INTEGER I, IACOL, IAROW, IB, ICOFF, ICTXT, IIA, IN,
134 $ IPT, IPV, IPW, JJA, JV, K, MYCOL, MYROW, NPCOL,
135 $ NPROW, NQ
136
137
138 INTEGER DESCV( DLEN_ )
139
140
143 $ pb_topset
144
145
146 INTEGER ICEIL, NUMROC
148
149
151
152
153
154
155
156 ictxt = desca( ctxt_ )
157 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
158
160 in =
min(
iceil( ia+m-k, desca( mb_ ) ) * desca( mb_ ), ia+m-1 )
161
162 icoff = mod( ja-1, desca( nb_ ) )
163 CALL infog2l( ia+m-k, ja, desca, nprow, npcol, myrow, mycol,
164 $ iia, jja, iarow, iacol )
165 nq =
numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
166 ipv = 1
167 ipt = ipv + nq * desca( mb_ )
168 ipw = ipt + desca( mb_ ) * desca( mb_ )
169 CALL pb_topget( ictxt, 'Broadcast', 'Rowwise', rowbtop )
170 CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', colbtop )
171 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', ' ' )
172 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', 'I-ring' )
173
174 CALL descset( descv, desca( mb_), n + icoff, desca( mb_ ),
175 $ desca( nb_ ), iarow, iacol, ictxt, desca( mb_ ) )
176
177
178
179 ib = in - ia - m + k + 1
180 jv = 1 + n - k + icoff
181
182
183
184 CALL pclarft(
'Backward',
'Rowwise', n-m+in-ia+1, ib, a, ia+m-k,
185 $ ja, desca, tau, work( ipt ), work( ipw ) )
186
187
188
189 CALL pclacpy(
'All', ib, n-m+in-ia+1, a, ia+m-k, ja, desca,
190 $ work( ipv ), 1, icoff+1, descv )
191 CALL pclaset(
'Upper', ib, ib, zero, one, work( ipv ), 1, jv,
192 $ descv )
193
194
195
196
197 CALL pclaset(
'All', ib, n-k, zero, zero, a, ia+m-k, ja,
198 $ desca )
199 CALL pclaset(
'Lower', ib-1, ib, zero, zero, a, ia+m-k+1,
200 $ ja+n-k, desca )
201
202
203
204 CALL pclarfb(
'Right',
'Conjugate transpose',
'Backward',
205 $ 'Rowwise', in-ia+1, n-m+in-ia+1, ib, work( ipv ), 1,
206 $ icoff+1, descv, work( ipt ), a, ia, ja, desca,
207 $ work( ipw ) )
208
209 descv( rsrc_ ) = mod( descv( rsrc_ ) + 1, nprow )
210
211
212
213 DO 10 i = in+1, ia+m-1, desca( mb_ )
214 ib =
min( ia+m-i, desca( mb_ ) )
215 jv = 1 + n - m + i - ia + icoff
216
217
218
219 CALL pclarft(
'Backward',
'Rowwise', n-m+i+ib-ia, ib, a, i, ja,
220 $ desca, tau, work( ipt ), work( ipw ) )
221
222
223
224 CALL pclacpy(
'All', ib, n-m+i+ib-ia, a, i, ja, desca,
225 $ work( ipv ), 1, icoff+1, descv )
226 CALL pclaset(
'Upper', ib, ib, zero, one, work( ipv ), 1, jv,
227 $ descv )
228
229
230
231
232 CALL pclaset(
'All', ib, n-m+i-ia, zero, zero, a, i, ja,
233 $ desca )
234 CALL pclaset(
'Lower', ib-1, ib, zero, zero, a, i+1,
235 $ ja+n-m+i-ia, desca )
236
237
238
239 CALL pclarfb(
'Right',
'Conjugate transpose',
'Backward',
240 $ 'Rowwise', i+ib-ia, n-m+i+ib-ia, ib, work( ipv ),
241 $ 1, icoff+1, descv, work( ipt ), a, ia, ja, desca,
242 $ work( ipw ) )
243
244 descv( rsrc_ ) = mod( descv( rsrc_ ) + 1, nprow )
245
246 10 CONTINUE
247
248 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', rowbtop )
249 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', colbtop )
250
251 RETURN
252
253
254
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 pclaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
subroutine pclacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
subroutine pclarfb(side, trans, direct, storev, m, n, k, v, iv, jv, descv, t, c, ic, jc, descc, work)
subroutine pclarft(direct, storev, n, k, v, iv, jv, descv, tau, t, work)