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 D( * ), E( * )
15 COMPLEX*16 A( * ), TAUP( * ), TAUQ( * ), X( * ), Y( * ),
16 $ WORK( * )
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
249
250 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
251 $ LLD_, MB_, M_, NB_, N_, RSRC_
252 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
253 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
254 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
255 COMPLEX*16 ONE, ZERO
256 parameter( one = ( 1.0d+0, 0.0d+0 ),
257 $ zero = ( 0.0d+0, 0.0d+0 ) )
258
259
260 INTEGER I, IACOL, IAROW, ICTXT, II, IPY, IW, J, JJ,
261 $ JWY, K, MYCOL, MYROW, NPCOL, NPROW
262 COMPLEX*16 ALPHA, TAU
263 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ),
264 $ DESCTP( DLEN_ ), DESCTQ( DLEN_ ),
265 $ DESCW( DLEN_ ), DESCWY( DLEN_ )
266
267
271
272
273 INTRINSIC dcmplx,
min, mod
274
275
276
277
278
279 IF( m.LE.0 .OR. n.LE.0 )
280 $ RETURN
281
282 ictxt = desca( ctxt_ )
283 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
284 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, ii, jj,
285 $ iarow, iacol )
286 ipy = desca( mb_ ) + 1
287 iw = mod( ia-1, desca( nb_ ) ) + 1
288 alpha = zero
289
290 CALL descset( descwy, 1, n+mod( ia-1, descy( nb_ ) ), 1,
291 $ desca( nb_ ), iarow, iacol, ictxt, 1 )
292 CALL descset( descw, desca( mb_ ), 1, desca( mb_ ), 1, iarow,
293 $ iacol, ictxt, desca( mb_ ) )
294 CALL descset( desctq, 1, ja+
min(m,n)-1, 1, desca( nb_ ), iarow,
295 $ desca( csrc_ ), desca( ctxt_ ), 1 )
296 CALL descset( desctp, ia+
min(m,n)-1, 1, desca( mb_ ), 1,
297 $ desca( rsrc_ ), iacol, desca( ctxt_ ),
298 $ desca( lld_ ) )
299
300 IF( m.GE.n ) THEN
301
302
303
304 CALL descset( descd, 1, ja+
min(m,n)-1, 1, desca( nb_ ), myrow,
305 $ desca( csrc_ ), desca( ctxt_ ), 1 )
306 CALL descset( desce, ia+
min(m,n)-1, 1, desca( mb_ ), 1,
307 $ desca( rsrc_ ), mycol, desca( ctxt_ ),
308 $ desca( lld_ ) )
309 DO 10 k = 1, nb
310 i = ia + k - 1
311 j = ja + k - 1
312 jwy = iw + k
313
314
315
316 IF( k.GT.1 ) THEN
317 CALL pzgemv( 'No transpose', m-k+1, k-1, -one, a, i, ja,
318 $ desca, y, iy, jy+k-1, descy, 1, one, a, i,
319 $ j, desca, 1 )
320 CALL pzgemv( 'No transpose', m-k+1, k-1, -one, x, ix+k-1,
321 $ jx, descx, a, ia, j, desca, 1, one, a, i, j,
322 $ desca, 1 )
323 CALL pzelset( a, i-1, j, desca, alpha )
324 END IF
325
326
327
328 CALL pzlarfg( m-k+1, alpha, i, j, a, i+1, j, desca, 1,
329 $ tauq )
330 CALL pdelset( d, 1, j, descd, dble( alpha ) )
331 CALL pzelset( a, i, j, desca, one )
332
333
334
335 CALL pzgemv( 'Conjugate transpose', m-k+1, n-k, one, a, i,
336 $ j+1, desca, a, i, j, desca, 1, zero,
337 $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
338 CALL pzgemv( 'Conjugate transpose', m-k+1, k-1, one, a, i,
339 $ ja, desca, a, i, j, desca, 1, zero, work, iw,
340 $ 1, descw, 1 )
341 CALL pzgemv( 'Conjugate transpose', k-1, n-k, -one, y, iy,
342 $ jy+k, descy, work, iw, 1, descw, 1, one,
343 $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
344 CALL pzgemv( 'Conjugate transpose', m-k+1, k-1, one, x,
345 $ ix+k-1, jx, descx, a, i, j, desca, 1, zero,
346 $ work, iw, 1, descw, 1 )
347 CALL pzgemv( 'Conjugate transpose', k-1, n-k, -one, a, ia,
348 $ j+1, desca, work, iw, 1, descw, 1, one,
349 $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
350
351 CALL pzelget(
'Rowwise',
' ', tau, tauq, 1, j, desctq )
352 CALL pzscal( n-k, tau, work( ipy ), 1, jwy, descwy,
353 $ descwy( m_ ) )
354 CALL pzlacgv( n-k, work( ipy ), 1, jwy, descwy,
355 $ descwy( m_ ) )
356 CALL pzcopy( n-k, work( ipy ), 1, jwy, descwy, descwy( m_ ),
357 $ y, iy+k-1, jy+k, descy, descy( m_ ) )
358
359
360
361 CALL pzlacgv( n-k, a, i, j+1, desca, desca( m_ ) )
362 CALL pzlacgv( k, a, i, ja, desca, desca( m_ ) )
363 CALL pzgemv( 'Conjugate transpose', k, n-k, -one, y, iy,
364 $ jy+k, descy, a, i, ja, desca, desca( m_ ), one,
365 $ a, i, j+1, desca, desca( m_ ) )
366 CALL pzlacgv( k, a, i, ja, desca, desca( m_ ) )
367 CALL pzlacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
368 CALL pzgemv( 'Conjugate transpose', k-1, n-k, -one, a, ia,
369 $ j+1, desca, x, ix+k-1, jx, descx, descx( m_ ),
370 $ one, a, i, j+1, desca, desca( m_ ) )
371 CALL pzlacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
372 CALL pzelset( a, i, j, desca, dcmplx( dble( alpha ) ) )
373
374
375
376 CALL pzlarfg( n-k, alpha, i, j+1, a, i,
377 $
min( j+2, n+ja-1 ), desca, desca( m_ ), taup )
378 CALL pdelset( e, i, 1, desce, dble( alpha ) )
379 CALL pzelset( a, i, j+1, desca, one )
380
381
382
383 CALL pzgemv( 'No transpose', m-k, n-k, one, a, i+1, j+1,
384 $ desca, a, i, j+1, desca, desca( m_ ), zero, x,
385 $ ix+k, jx+k-1, descx, 1 )
386 CALL pzgemv( 'No transpose', k, n-k, one, y, iy, jy+k,
387 $ descy, a, i, j+1, desca, desca( m_ ), zero,
388 $ work, iw, 1, descw, 1 )
389 CALL pzgemv( 'No transpose', m-k, k, -one, a, i+1, ja,
390 $ desca, work, iw, 1, descw, 1, one, x, ix+k,
391 $ jx+k-1, descx, 1 )
392 CALL pzgemv( 'No transpose', k-1, n-k, one, a, ia, j+1,
393 $ desca, a, i, j+1, desca, desca( m_ ), zero,
394 $ work, iw, 1, descw, 1 )
395 CALL pzgemv( 'No transpose', m-k, k-1, -one, x, ix+k, jx,
396 $ descx, work, iw, 1, descw, 1, one, x, ix+k,
397 $ jx+k-1, descx, 1 )
398
399 CALL pzelget(
'Columnwise',
' ', tau, taup, i, 1, desctp )
400 CALL pzscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
401 CALL pzlacgv( n-k, a, i, j+1, desca, desca( m_ ) )
402 10 CONTINUE
403
404 ELSE
405
406
407
408 CALL descset( descd, ia+
min(m,n)-1, 1, desca( mb_ ), 1,
409 $ desca( rsrc_ ), mycol, desca( ctxt_ ),
410 $ desca( lld_ ) )
411 CALL descset( desce, 1, ja+
min(m,n)-1, 1, desca( nb_ ), myrow,
412 $ desca( csrc_ ), desca( ctxt_ ), 1 )
413 DO 20 k = 1, nb
414 i = ia + k - 1
415 j = ja + k - 1
416 jwy = iw + k
417
418
419
420 CALL pzlacgv( n-k+1, a, i, j, desca, desca( m_ ) )
421 IF( k.GT.1 ) THEN
422 CALL pzlacgv( k-1, a, i, ja, desca, desca( m_ ) )
423 CALL pzgemv( 'Conjugate transpose', k-1, n-k+1, -one, y,
424 $ iy, jy+k-1, descy, a, i, ja, desca,
425 $ desca( m_ ), one, a, i, j, desca,
426 $ desca( m_ ) )
427 CALL pzlacgv( k-1, a, i, ja, desca, desca( m_ ) )
428 CALL pzlacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
429 CALL pzgemv( 'Conjugate transpose', k-1, n-k+1, -one, a,
430 $ ia, j, desca, x, ix+k-1, jx, descx,
431 $ descx( m_ ), one, a, i, j, desca,
432 $ desca( m_ ) )
433 CALL pzlacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
434 CALL pzelset( a, i, j-1, desca, dcmplx( dble( alpha ) ) )
435 END IF
436
437
438
439 CALL pzlarfg( n-k+1, alpha, i, j, a, i, j+1, desca,
440 $ desca( m_ ), taup )
441 CALL pdelset( d, i, 1, descd, dble( alpha ) )
442 CALL pzelset( a, i, j, desca, one )
443
444
445
446 CALL pzgemv( 'No transpose', m-k, n-k+1, one, a, i+1, j,
447 $ desca, a, i, j, desca, desca( m_ ), zero, x,
448 $ ix+k, jx+k-1, descx, 1 )
449 CALL pzgemv( 'No transpose', k-1, n-k+1, one, y, iy, jy+k-1,
450 $ descy, a, i, j, desca, desca( m_ ), zero,
451 $ work, iw, 1, descw, 1 )
452 CALL pzgemv( 'No transpose', m-k, k-1, -one, a, i+1, ja,
453 $ desca, work, iw, 1, descw, 1, one, x, ix+k,
454 $ jx+k-1, descx, 1 )
455 CALL pzgemv( 'No transpose', k-1, n-k+1, one, a, ia, j,
456 $ desca, a, i, j, desca, desca( m_ ), zero,
457 $ work, iw, 1, descw, 1 )
458 CALL pzgemv( 'No transpose', m-k, k-1, -one, x, ix+k, jx,
459 $ descx, work, iw, 1, descw, 1, one, x, ix+k,
460 $ jx+k-1, descx, 1 )
461
462 CALL pzelget(
'Columnwise',
' ', tau, taup, i, 1, desctp )
463 CALL pzscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
464 CALL pzlacgv( n-k+1, a, i, j, desca, desca( m_ ) )
465
466
467
468 CALL pzgemv( 'No transpose', m-k, k-1, -one, a, i+1, ja,
469 $ desca, y, iy, jy+k-1, descy, 1, one, a, i+1, j,
470 $ desca, 1 )
471 CALL pzgemv( 'No transpose', m-k, k, -one, x, ix+k, jx,
472 $ descx, a, ia, j, desca, 1, one, a, i+1, j,
473 $ desca, 1 )
474 CALL pzelset( a, i, j, desca, alpha )
475
476
477
478 CALL pzlarfg( m-k, alpha, i+1, j, a,
min( i+2, m+ia-1 ),
479 $ j, desca, 1, tauq )
480 CALL pdelset( e, 1, j, desce, dble( alpha ) )
481 CALL pzelset( a, i+1, j, desca, one )
482
483
484
485 CALL pzgemv( 'Conjugate transpose', m-k, n-k, one, a, i+1,
486 $ j+1, desca, a, i+1, j, desca, 1, zero,
487 $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
488 CALL pzgemv( 'Conjugate transpose', m-k, k-1, one, a, i+1,
489 $ ja, desca, a, i+1, j, desca, 1, zero, work, iw,
490 $ 1, descw, 1 )
491 CALL pzgemv( 'Conjugate transpose', k-1, n-k, -one, y, iy,
492 $ jy+k, descy, work, iw, 1, descw, 1, one,
493 $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
494 CALL pzgemv( 'Conjugate transpose', m-k, k, one, x, ix+k,
495 $ jx, descx, a, i+1, j, desca, 1, zero, work, iw,
496 $ 1, descw, 1 )
497 CALL pzgemv( 'Conjugate transpose', k, n-k, -one, a, ia,
498 $ j+1, desca, work, iw, 1, descw, 1, one,
499 $ work( ipy ), 1, jwy, descwy, descwy( m_ ) )
500
501 CALL pzelget(
'Rowwise',
' ', tau, tauq, 1, j, desctq )
502 CALL pzscal( n-k, tau, work( ipy ), 1, jwy, descwy,
503 $ descwy( m_ ) )
504 CALL pzlacgv( n-k, work( ipy ), 1, jwy, descwy,
505 $ descwy( m_ ) )
506 CALL pzcopy( n-k, work( ipy ), 1, jwy, descwy, descwy( m_ ),
507 $ y, iy+k-1, jy+k, descy, descy( m_ ) )
508 20 CONTINUE
509 END IF
510
511 RETURN
512
513
514
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 pdelset(a, ia, ja, desca, alpha)
subroutine pzelget(scope, top, alpha, a, ia, ja, desca)
subroutine pzelset(a, ia, ja, desca, alpha)
subroutine pzlacgv(n, x, ix, jx, descx, incx)
subroutine pzlarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)