2
3
4
5
6
7
8
9 INTEGER IA, JA, M, N
10
11
12 INTEGER DESCA( * )
13 DOUBLE PRECISION 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
123
124 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
125 $ LLD_, MB_, M_, NB_, N_, RSRC_
126 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
127 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
128 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
129 DOUBLE PRECISION ONE, ZERO
130 parameter( one = 1.0d+0, zero = 0.0d+0 )
131
132
133 CHARACTER COLBTOP, ROWBTOP
134 INTEGER IACOL, IAROW, ICTXT, IIA, IPT, IPV, IPW, IROFF,
135 $ IV, J, JB, JJA, JN, K, MP, MYCOL, MYROW, NPCOL,
136 $ NPROW
137
138
139 INTEGER DESCV( DLEN_ )
140
141
144 $ pb_topset
145
146
147 INTEGER ICEIL, NUMROC
149
150
152
153
154
155
156
157 ictxt = desca( ctxt_ )
158 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
159
161 jn =
min(
iceil( ja+n-k, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
162
163 iroff = mod( ia-1, desca( mb_ ) )
164 CALL infog2l( ia, ja+n-k, desca, nprow, npcol, myrow, mycol,
165 $ iia, jja, iarow, iacol )
166 mp =
numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
167 ipv = 1
168 ipt = ipv + mp * desca( nb_ )
169 ipw = ipt + desca( nb_ ) * desca( nb_ )
170 CALL pb_topget( ictxt, 'Broadcast', 'Rowwise', rowbtop )
171 CALL pb_topget( ictxt, 'Broadcast', 'Columnwise', colbtop )
172 CALL pb_topset( ictxt, 'Broadcast', 'Rowwise', 'I-ring' )
173 CALL pb_topset( ictxt, 'Broadcast', 'Columnwise', ' ' )
174
175 CALL descset( descv, m+iroff, desca( nb_ ), desca( mb_ ),
176 $ desca( nb_ ), iarow, iacol, ictxt,
max( 1, mp ) )
177
178
179
180 iv = 1 + m - k + iroff
181 jb = jn - ja - n + k + 1
182
183
184
185 CALL pdlarft(
'Backward',
'Columnwise', m-n+jn-ja+1, jb, a, ia,
186 $ ja+n-k, desca, tau, work( ipt ), work( ipw ) )
187
188
189
190 CALL pdlacpy(
'All', m-n+jn-ja+1, jb, a, ia, ja+n-k, desca,
191 $ work( ipv ), iroff+1, 1, descv )
192 CALL pdlaset(
'Lower', jb, jb, zero, one, work( ipv ), iv,
193 $ 1, descv )
194
195
196
197
198 CALL pdlaset(
'All', m-k, jb, zero, zero, a, ia, ja+n-k,
199 $ desca )
200 CALL pdlaset(
'Upper', jb, jb-1, zero, zero, a, ia+m-k,
201 $ ja+n-k+1, desca )
202
203
204
205 CALL pdlarfb(
'Left',
'No transpose',
'Backward',
'Columnwise',
206 $ m-n+jn-ja+1, jn-ja+1, jb, work( ipv ), iroff+1, 1,
207 $ descv, work( ipt ), a, ia, ja, desca, work( ipw ) )
208
209 descv( csrc_ ) = mod( descv( csrc_ ) + 1, npcol )
210
211
212
213 DO 10 j = jn+1, ja+n-1, desca( nb_ )
214 jb =
min( ja+n-j, desca( nb_ ) )
215 iv = 1 + m - n + j - ja + iroff
216
217
218
219 CALL pdlarft(
'Backward',
'Columnwise', m-n+j+jb-ja, jb, a, ia,
220 $ j, desca, tau, work( ipt ), work( ipw ) )
221
222
223
224 CALL pdlacpy(
'All', m-n+j+jb-ja, jb, a, ia, j, desca,
225 $ work( ipv ), iroff+1, 1, descv )
226 CALL pdlaset(
'Lower', jb, jb, zero, one, work( ipv ), iv,
227 $ 1, descv )
228
229
230
231
232 CALL pdlaset(
'All', m-n+j-ja, jb, zero, zero, a, ia, j,
233 $ desca )
234 CALL pdlaset(
'Upper', jb, jb-1, zero, zero, a, ia+m-n+j-ja,
235 $ j+1, desca )
236
237
238
239 CALL pdlarfb(
'Left',
'No transpose',
'Backward',
'Columnwise',
240 $ m-n+j+jb-ja, j+jb-ja, jb, work( ipv ), iroff+1,
241 $ 1, descv, work( ipt ), a, ia, ja, desca,
242 $ work( ipw ) )
243
244 descv( csrc_ ) = mod( descv( csrc_ ) + 1, npcol )
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 pdlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
subroutine pdlacpy(uplo, m, n, a, ia, ja, desca, b, ib, jb, descb)
subroutine pdlarfb(side, trans, direct, storev, m, n, k, v, iv, jv, descv, t, c, ic, jc, descc, work)
subroutine pdlarft(direct, storev, n, k, v, iv, jv, descv, tau, t, work)