3
4
5
6
7
8
9
10 INTEGER IA, IX, IY, JA, JX, JY, M, N, NB
11
12
13 INTEGER DESCA( * ), DESCX( * ), DESCY( * )
14 DOUBLE PRECISION A( * ), D( * ), E( * ), TAUP( * ),
15 $ TAUQ( * ), X( * ), Y( * ), WORK( * )
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
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
249 $ LLD_, MB_, M_, NB_, N_, RSRC_
250 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
251 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
252 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
253 DOUBLE PRECISION ONE, ZERO
254 parameter( one = 1.0d+0, zero = 0.0d+0 )
255
256
257 INTEGER I, IACOL, IAROW, ICTXT, II, IPY, IW, J, JJ,
258 $ JWY, K, MYCOL, MYROW, NPCOL, NPROW
259 DOUBLE PRECISION ALPHA, TAU
260 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ),
261 $ DESCTP( DLEN_ ), DESCTQ( DLEN_ ),
262 $ DESCW( DLEN_ ), DESCWY( DLEN_ )
263
264
267 $ pdscal
268
269
271
272
273
274
275
276 IF( m.LE.0 .OR. n.LE.0 )
277 $ RETURN
278
279 ictxt = desca( ctxt_ )
280 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
281 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, ii, jj,
282 $ iarow, iacol )
283 ipy = desca( mb_ ) + 1
284 iw = mod( ia-1, desca( nb_ ) ) + 1
285 alpha = zero
286
287 CALL descset( descwy, 1, n+mod( ia-1, descy( nb_ ) ), 1,
288 $ desca( nb_ ), iarow, iacol, ictxt, 1 )
289 CALL descset( descw, desca( mb_ ), 1, desca( mb_ ), 1, iarow,
290 $ iacol, ictxt, desca( mb_ ) )
291 CALL descset( desctq, 1, ja+
min(m,n)-1, 1, desca( nb_ ), iarow,
292 $ desca( csrc_ ), desca( ctxt_ ), 1 )
293 CALL descset( desctp, ia+
min(m,n)-1, 1, desca( mb_ ), 1,
294 $ desca( rsrc_ ), iacol, desca( ctxt_ ),
295 $ desca( lld_ ) )
296
297 IF( m.GE.n ) THEN
298
299
300
301 CALL descset( descd, 1, ja+
min(m,n)-1, 1, desca( nb_ ), myrow,
302 $ desca( csrc_ ), desca( ctxt_ ), 1 )
303 CALL descset( desce, ia+
min(m,n)-1, 1, desca( mb_ ), 1,
304 $ desca( rsrc_ ), mycol, desca( ctxt_ ),
305 $ desca( lld_ ) )
306 DO 10 k = 1, nb
307 i = ia + k - 1
308 j = ja + k - 1
309 jwy = iw + k
310
311
312
313 IF( k.GT.1 ) THEN
314 CALL pdgemv( 'No transpose', m-k+1, k-1, -one, a, i, ja,
315 $ desca, y, iy, jy+k-1, descy, 1, one, a, i,
316 $ j, desca, 1 )
317 CALL pdgemv( 'No transpose', m-k+1, k-1, -one, x, ix+k-1,
318 $ jx, descx, a, ia, j, desca, 1, one, a, i, j,
319 $ desca, 1 )
320 CALL pdelset( a, i-1, j, desca, alpha )
321 END IF
322
323
324
325 CALL pdlarfg( m-k+1, alpha, i, j, a, i+1, j, desca, 1,
326 $ tauq )
327 CALL pdelset( d, 1, j, descd, alpha )
328 CALL pdelset( a, i, j, desca, one )
329
330
331
332 CALL pdgemv( 'Transpose', m-k+1, n-k, one, a, i, j+1, desca,
333 $ a, i, j, desca, 1, zero, work( ipy ), 1, jwy,
334 $ descwy, descwy( m_ ) )
335 CALL pdgemv( 'Transpose', m-k+1, k-1, one, a, i, ja, desca,
336 $ a, i, j, desca, 1, zero, work, iw, 1, descw,
337 $ 1 )
338 CALL pdgemv( 'Transpose', k-1, n-k, -one, y, iy, jy+k,
339 $ descy, work, iw, 1, descw, 1, one, work( ipy ),
340 $ 1, jwy, descwy, descwy( m_ ) )
341 CALL pdgemv( 'Transpose', m-k+1, k-1, one, x, ix+k-1, jx,
342 $ descx, a, i, j, desca, 1, zero, work, iw, 1,
343 $ descw, 1 )
344 CALL pdgemv( 'Transpose', k-1, n-k, -one, a, ia, j+1, desca,
345 $ work, iw, 1, descw, 1, one, work( ipy ), 1,
346 $ jwy, descwy, descwy( m_ ) )
347
348 CALL pdelget(
'Rowwise',
' ', tau, tauq, 1, j, desctq )
349 CALL pdscal( n-k, tau, work( ipy ), 1, jwy, descwy,
350 $ descwy( m_ ) )
351 CALL pdcopy( n-k, work( ipy ), 1, jwy, descwy, descwy( m_ ),
352 $ y, iy+k-1, jy+k, descy, descy( m_ ) )
353
354
355
356 CALL pdgemv( 'Transpose', k, n-k, -one, y, iy, jy+k, descy,
357 $ a, i, ja, desca, desca( m_ ), one, a, i, j+1,
358 $ desca, desca( m_ ) )
359 CALL pdgemv( 'Transpose', k-1, n-k, -one, a, ia, j+1, desca,
360 $ x, ix+k-1, jx, descx, descx( m_ ), one, a, i,
361 $ j+1, desca, desca( m_ ) )
362 CALL pdelset( a, i, j, desca, alpha )
363
364
365
366 CALL pdlarfg( n-k, alpha, i, j+1, a, i,
367 $
min( j+2, n+ja-1 ), desca, desca( m_ ), taup )
368 CALL pdelset( e, i, 1, desce, alpha )
369 CALL pdelset( a, i, j+1, desca, one )
370
371
372
373 CALL pdgemv( 'No transpose', m-k, n-k, one, a, i+1, j+1,
374 $ desca, a, i, j+1, desca, desca( m_ ), zero, x,
375 $ ix+k, jx+k-1, descx, 1 )
376 CALL pdgemv( 'No transpose', k, n-k, one, y, iy, jy+k,
377 $ descy, a, i, j+1, desca, desca( m_ ), zero,
378 $ work, iw, 1, descw, 1 )
379 CALL pdgemv( 'No transpose', m-k, k, -one, a, i+1, ja,
380 $ desca, work, iw, 1, descw, 1, one, x, ix+k,
381 $ jx+k-1, descx, 1 )
382 CALL pdgemv( 'No transpose', k-1, n-k, one, a, ia, j+1,
383 $ desca, a, i, j+1, desca, desca( m_ ), zero,
384 $ work, iw, 1, descw, 1 )
385 CALL pdgemv( 'No transpose', m-k, k-1, -one, x, ix+k, jx,
386 $ descx, work, iw, 1, descw, 1, one, x, ix+k,
387 $ jx+k-1, descx, 1 )
388
389 CALL pdelget(
'Columnwise',
' ', tau, taup, i, 1, desctp )
390 CALL pdscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
391 10 CONTINUE
392
393 ELSE
394
395
396
397 CALL descset( descd, ia+
min(m,n)-1, 1, desca( mb_ ), 1,
398 $ desca( rsrc_ ), mycol, desca( ctxt_ ),
399 $ desca( lld_ ) )
400 CALL descset( desce, 1, ja+
min(m,n)-1, 1, desca( nb_ ), myrow,
401 $ desca( csrc_ ), desca( ctxt_ ), 1 )
402 DO 20 k = 1, nb
403 i = ia + k - 1
404 j = ja + k - 1
405 jwy = iw + k
406
407
408
409 IF( k.GT.1 ) THEN
410 CALL pdgemv( 'Transpose', k-1, n-k+1, -one, y, iy,
411 $ jy+k-1, descy, a, i, ja, desca, desca( m_ ),
412 $ one, a, i, j, desca, desca( m_ ) )
413 CALL pdgemv( 'Transpose', k-1, n-k+1, -one, a, ia, j,
414 $ desca, x, ix+k-1, jx, descx, descx( m_ ),
415 $ one, a, i, j, desca, desca( m_ ) )
416 CALL pdelset( a, i, j-1, desca, alpha )
417 END IF
418
419
420
421 CALL pdlarfg( n-k+1, alpha, i, j, a, i, j+1, desca,
422 $ desca( m_ ), taup )
423 CALL pdelset( d, i, 1, descd, alpha )
424 CALL pdelset( a, i, j, desca, one )
425
426
427
428 CALL pdgemv( 'No transpose', m-k, n-k+1, one, a, i+1, j,
429 $ desca, a, i, j, desca, desca( m_ ), zero, x,
430 $ ix+k, jx+k-1, descx, 1 )
431 CALL pdgemv( 'No transpose', k-1, n-k+1, one, y, iy, jy+k-1,
432 $ descy, a, i, j, desca, desca( m_ ), zero,
433 $ work, iw, 1, descw, 1 )
434 CALL pdgemv( 'No transpose', m-k, k-1, -one, a, i+1, ja,
435 $ desca, work, iw, 1, descw, 1, one, x, ix+k,
436 $ jx+k-1, descx, 1 )
437 CALL pdgemv( 'No transpose', k-1, n-k+1, one, a, ia, j,
438 $ desca, a, i, j, desca, desca( m_ ), zero,
439 $ work, iw, 1, descw, 1 )
440 CALL pdgemv( 'No transpose', m-k, k-1, -one, x, ix+k, jx,
441 $ descx, work, iw, 1, descw, 1, one, x, ix+k,
442 $ jx+k-1, descx, 1 )
443
444 CALL pdelget(
'Columnwise',
' ', tau, taup, i, 1, desctp )
445 CALL pdscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
446
447
448
449 CALL pdgemv( 'No transpose', m-k, k-1, -one, a, i+1, ja,
450 $ desca, y, iy, jy+k-1, descy, 1, one, a, i+1, j,
451 $ desca, 1 )
452 CALL pdgemv( 'No transpose', m-k, k, -one, x, ix+k, jx,
453 $ descx, a, ia, j, desca, 1, one, a, i+1, j,
454 $ desca, 1 )
455 CALL pdelset( a, i, j, desca, alpha )
456
457
458
459 CALL pdlarfg( m-k, alpha, i+1, j, a,
min( i+2, m+ia-1 ),
460 $ j, desca, 1, tauq )
461 CALL pdelset( e, 1, j, desce, alpha )
462 CALL pdelset( a, i+1, j, desca, one )
463
464
465
466 CALL pdgemv( 'Transpose', m-k, n-k, one, a, i+1, j+1, desca,
467 $ a, i+1, j, desca, 1, zero, work( ipy ), 1,
468 $ jwy, descwy, descwy( m_ ) )
469 CALL pdgemv( 'Transpose', m-k, k-1, one, a, i+1, ja, desca,
470 $ a, i+1, j, desca, 1, zero, work, iw, 1, descw,
471 $ 1 )
472 CALL pdgemv( 'Transpose', k-1, n-k, -one, y, iy, jy+k,
473 $ descy, work, iw, 1, descw, 1, one, work( ipy ),
474 $ 1, jwy, descwy, descwy( m_ ) )
475 CALL pdgemv( 'Transpose', m-k, k, one, x, ix+k, jx, descx,
476 $ a, i+1, j, desca, 1, zero, work, iw, 1, descw,
477 $ 1 )
478 CALL pdgemv( 'Transpose', k, n-k, -one, a, ia, j+1, desca,
479 $ work, iw, 1, descw, 1, one, work( ipy ), 1,
480 $ jwy, descwy, descwy( m_ ) )
481
482 CALL pdelget(
'Rowwise',
' ', tau, tauq, 1, j, desctq )
483 CALL pdscal( n-k, tau, work( ipy ), 1, jwy, descwy,
484 $ descwy( m_ ) )
485 CALL pdcopy( n-k, work( ipy ), 1, jwy, descwy, descwy( m_ ),
486 $ y, iy+k-1, jy+k, descy, descy( m_ ) )
487 20 CONTINUE
488 END IF
489
490 RETURN
491
492
493
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
subroutine pdelget(scope, top, alpha, a, ia, ja, desca)
subroutine pdelset(a, ia, ja, desca, alpha)
subroutine pdlarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)