2
3
4
5
6
7
8
9 INTEGER IA, JA, M, N
10
11
12 INTEGER DESCA( * )
13 COMPLEX*16 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*16 ONE, ZERO
128 parameter( one = ( 1.0d+0, 0.0d+0 ),
129 $ zero = ( 0.0d+0, 0.0d+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 pzlarft(
'Forward',
'Rowwise', n-j+ja, ib, a, i, j, desca,
185 $ tau, work( ipt ), work( ipw ) )
186
187
188
189 CALL pzlacpy(
'Upper', ib, n-j+ja, a, i, j, desca, work( ipv ),
190 $ 1, jv, descv )
191 CALL pzlaset(
'Lower', ib, n-j+ja, zero, one, work( ipv ), 1,
192 $ jv, descv )
193
194
195
196
197 CALL pzlaset(
'Upper', ib, n-j+ja-1, zero, zero, a, i, j+1,
198 $ desca )
199
200
201
202 CALL pzlarfb(
'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 pzlarft(
'Forward',
'Rowwise', n, ib, a, ia, ja, desca, tau,
218 $ work( ipt ), work( ipw ) )
219
220
221
222 CALL pzlacpy(
'Upper', ib, n, a, ia, ja, desca, work( ipv ), 1,
223 $ icoff+1, descv )
224 CALL pzlaset(
'Lower', ib, n, zero, one, work, 1, icoff+1, descv )
225
226
227
228
229 CALL pzlaset(
'Upper', ib, n-1, zero, zero, a, ia, ja+1, desca )
230
231
232
233 CALL pzlarfb(
'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 pzlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
subroutine pzlacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
subroutine pzlarfb(side, trans, direct, storev, m, n, k, v, iv, jv, descv, t, c, ic, jc, descc, work)
subroutine pzlarft(direct, storev, n, k, v, iv, jv, descv, tau, t, work)