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
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 REAL ONE, ZERO
130 parameter( one = 1.0e+0, zero = 0.0e+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 pslarft(
'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 pslacpy(
'All', m-n+jn-ja+1, jb, a, ia, ja+n-k, desca,
191 $ work( ipv ), iroff+1, 1, descv )
192 CALL pslaset(
'Lower', jb, jb, zero, one, work( ipv ), iv,
193 $ 1, descv )
194
195
196
197
198 CALL pslaset(
'All', m-k, jb, zero, zero, a, ia, ja+n-k,
199 $ desca )
200 CALL pslaset(
'Upper', jb, jb-1, zero, zero, a, ia+m-k,
201 $ ja+n-k+1, desca )
202
203
204
205 CALL pslarfb(
'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 pslarft(
'Backward',
'Columnwise', m-n+j+jb-ja, jb, a, ia,
220 $ j, desca, tau, work( ipt ), work( ipw ) )
221
222
223
224 CALL pslacpy(
'All', m-n+j+jb-ja, jb, a, ia, j, desca,
225 $ work( ipv ), iroff+1, 1, descv )
226 CALL pslaset(
'Lower', jb, jb, zero, one, work( ipv ), iv,
227 $ 1, descv )
228
229
230
231
232 CALL pslaset(
'All', m-n+j-ja, jb, zero, zero, a, ia, j,
233 $ desca )
234 CALL pslaset(
'Upper', jb, jb-1, zero, zero, a, ia+m-n+j-ja,
235 $ j+1, desca )
236
237
238
239 CALL pslarfb(
'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 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)