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 REAL D( * ), E( * )
15 COMPLEX 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 ONE, ZERO
256 parameter( one = ( 1.0e+0, 0.0e+0 ),
257 $ zero = ( 0.0e+0, 0.0e+0 ) )
258
259
260 INTEGER I, IACOL, IAROW, ICTXT, II, IPY, IW, J, JJ,
261 $ JWY, K, MYCOL, MYROW, NPCOL, NPROW
262 COMPLEX ALPHA, TAU
263 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ),
264 $ DESCTP( DLEN_ ), DESCTQ( DLEN_ ),
265 $ DESCW( DLEN_ ), DESCWY( DLEN_ )
266
267
271
272
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 pcgemv( '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 pcgemv( '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 pcelset( a, i-1, j, desca, alpha )
324 END IF
325
326
327
328 CALL pclarfg( m-k+1, alpha, i, j, a, i+1, j, desca, 1,
329 $ tauq )
330 CALL pselset( d, 1, j, descd, real( alpha ) )
331 CALL pcelset( a, i, j, desca, one )
332
333
334
335 CALL pcgemv( '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 pcgemv( '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 pcgemv( '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 pcgemv( '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 pcgemv( '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 pcelget(
'Rowwise',
' ', tau, tauq, 1, j, desctq )
352 CALL pcscal( n-k, tau, work( ipy ), 1, jwy, descwy,
353 $ descwy( m_ ) )
354 CALL pclacgv( n-k, work( ipy ), 1, jwy, descwy,
355 $ descwy( m_ ) )
356 CALL pccopy( n-k, work( ipy ), 1, jwy, descwy, descwy( m_ ),
357 $ y, iy+k-1, jy+k, descy, descy( m_ ) )
358
359
360
361 CALL pclacgv( n-k, a, i, j+1, desca, desca( m_ ) )
362 CALL pclacgv( k, a, i, ja, desca, desca( m_ ) )
363 CALL pcgemv( '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 pclacgv( k, a, i, ja, desca, desca( m_ ) )
367 CALL pclacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
368 CALL pcgemv( '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 pclacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
373
374
375
376 CALL pclarfg( n-k, alpha, i, j+1, a, i,
377 $
min( j+2, n+ja-1 ), desca, desca( m_ ), taup )
378 CALL pselset( e, i, 1, desce, real( alpha ) )
379 CALL pcelset( a, i, j+1, desca, one )
380
381
382
383 CALL pcgemv( '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 pcgemv( '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 pcgemv( '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 pcgemv( '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 pcgemv( '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 pcelget(
'Columnwise',
' ', tau, taup, i, 1, desctp )
400 CALL pcscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
401 CALL pclacgv( 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 pclacgv( n-k+1, a, i, j, desca, desca( m_ ) )
421 IF( k.GT.1 ) THEN
422 CALL pclacgv( k-1, a, i, ja, desca, desca( m_ ) )
423 CALL pcgemv( '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 pclacgv( k-1, a, i, ja, desca, desca( m_ ) )
428 CALL pclacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
429 CALL pcgemv( '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 pclacgv( k-1, x, ix+k-1, jx, descx, descx( m_ ) )
434 CALL pcelset( a, i, j-1, desca,
cmplx( real( alpha ) ) )
435 END IF
436
437
438
439 CALL pclarfg( n-k+1, alpha, i, j, a, i, j+1, desca,
440 $ desca( m_ ), taup )
441 CALL pselset( d, i, 1, descd, real( alpha ) )
442 CALL pcelset( a, i, j, desca, one )
443
444
445
446 CALL pcgemv( '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 pcgemv( '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 pcgemv( '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 pcgemv( '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 pcgemv( '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 pcelget(
'Columnwise',
' ', tau, taup, i, 1, desctp )
463 CALL pcscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
464 CALL pclacgv( n-k+1, a, i, j, desca, desca( m_ ) )
465
466
467
468 CALL pcgemv( '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 pcgemv( '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 pcelset( a, i, j, desca, alpha )
475
476
477
478 CALL pclarfg( m-k, alpha, i+1, j, a,
min( i+2, m+ia-1 ),
479 $ j, desca, 1, tauq )
480 CALL pselset( e, 1, j, desce, real( alpha ) )
481 CALL pcelset( a, i+1, j, desca, one )
482
483
484
485 CALL pcgemv( '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 pcgemv( '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 pcgemv( '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 pcgemv( '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 pcgemv( '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 pcelget(
'Rowwise',
' ', tau, tauq, 1, j, desctq )
502 CALL pcscal( n-k, tau, work( ipy ), 1, jwy, descwy,
503 $ descwy( m_ ) )
504 CALL pclacgv( n-k, work( ipy ), 1, jwy, descwy,
505 $ descwy( m_ ) )
506 CALL pccopy( 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 pcelget(scope, top, alpha, a, ia, ja, desca)
subroutine pcelset(a, ia, ja, desca, alpha)
subroutine pclacgv(n, x, ix, jx, descx, incx)
subroutine pclarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)
subroutine pselset(a, ia, ja, desca, alpha)