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, IL, IN,
134 $ IPT, IPV, IPW, J, JJA, JV, K, MYCOL, MYROW,
135 $ NPCOL, NPROW, NQ
136
137
138 INTEGER DESCV( DLEN_ )
139
140
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, desca( mb_ ) ) * desca( mb_ ), ia+k-1 )
160 il =
max( ( (ia+k-2) / desca( mb_ ) ) * desca( mb_ ) + 1, ia )
161
162 icoff = mod( ja-1, desca( nb_ ) )
163 CALL infog2l( il, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
164 $ 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', 'D-ring' )
173
174 CALL descset( descv, desca( mb_ ), n + icoff, desca( mb_ ),
175 $ desca( nb_ ), iarow, iacol, ictxt, desca( mb_ ) )
176
177 DO 10 i = il, in+1, -desca( mb_ )
178 ib =
min( ia+k-i, desca( mb_ ) )
179 j = ja + i - ia
180 jv = 1 + i - ia + icoff
181
182
183
184 CALL pclarft(
'Forward',
'Rowwise', n-j+ja, ib, a, i, j, desca,
185 $ tau, work( ipt ), work( ipw ) )
186
187
188
189 CALL pclacpy(
'Upper', ib, n-j+ja, a, i, j, desca, work( ipv ),
190 $ 1, jv, descv )
191 CALL pclaset(
'Lower', ib, n-j+ja, zero, one, work( ipv ), 1,
192 $ jv, descv )
193
194
195
196
197 CALL pclaset(
'Upper', ib, n-j+ja-1, zero, zero, a, i, j+1,
198 $ desca )
199
200
201
202 CALL pclarfb(
'Right',
'Conjugate transpose',
'Forward',
203 $ 'Rowwise', m-i+ia, n-j+ja, ib, work( ipv ), 1,
204 $ jv, descv, work( ipt ), a, i, j, desca,
205 $ work( ipw ) )
206
207 descv( rsrc_ ) = mod( descv( rsrc_ ) + nprow - 1, nprow )
208
209 10 CONTINUE
210
211
212
213 ib = in - ia + 1
214
215
216
217 CALL pclarft(
'Forward',
'Rowwise', n, ib, a, ia, ja, desca, tau,
218 $ work( ipt ), work( ipw ) )
219
220
221
222 CALL pclacpy(
'Upper', ib, n, a, ia, ja, desca, work( ipv ), 1,
223 $ icoff+1, descv )
224 CALL pclaset(
'Lower', ib, n, zero, one, work, 1, icoff+1, descv )
225
226
227
228
229 CALL pclaset(
'Upper', ib, n-1, zero, zero, a, ia, ja+1, desca )
230
231
232
233 CALL pclarfb(
'Right',
'Conjugate transpose',
'Forward',
234 $ 'Rowwise', m, n, ib, work( ipv ), 1, icoff+1, descv,
235 $ work( ipt ), a, ia, ja, desca, work( ipw ) )
236
237 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', rowbtop )
238 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', colbtop )
239
240 RETURN
241
242
243
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)