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 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 REAL ONE, ZERO
254 parameter( one = 1.0e+0, zero = 0.0e+0 )
255
256
257 INTEGER I, IACOL, IAROW, ICTXT, II, IPY, IW, J, JJ,
258 $ JWY, K, MYCOL, MYROW, NPCOL, NPROW
259 REAL ALPHA, TAU
260 INTEGER DESCD( DLEN_ ), DESCE( DLEN_ ),
261 $ DESCTP( DLEN_ ), DESCTQ( DLEN_ ),
262 $ DESCW( DLEN_ ), DESCWY( DLEN_ )
263
264
267 $ psscal
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 psgemv( '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 psgemv( '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 pselset( a, i-1, j, desca, alpha )
321 END IF
322
323
324
325 CALL pslarfg( m-k+1, alpha, i, j, a, i+1, j, desca, 1,
326 $ tauq )
327 CALL pselset( d, 1, j, descd, alpha )
328 CALL pselset( a, i, j, desca, one )
329
330
331
332 CALL psgemv( '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 psgemv( '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 psgemv( '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 psgemv( '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 psgemv( '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 pselget(
'Rowwise',
' ', tau, tauq, 1, j, desctq )
349 CALL psscal( n-k, tau, work( ipy ), 1, jwy, descwy,
350 $ descwy( m_ ) )
351 CALL pscopy( n-k, work( ipy ), 1, jwy, descwy, descwy( m_ ),
352 $ y, iy+k-1, jy+k, descy, descy( m_ ) )
353
354
355
356 CALL psgemv( '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 psgemv( '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 pselset( a, i, j, desca, alpha )
363
364
365
366 CALL pslarfg( n-k, alpha, i, j+1, a, i,
367 $
min( j+2, n+ja-1 ), desca, desca( m_ ), taup )
368 CALL pselset( e, i, 1, desce, alpha )
369 CALL pselset( a, i, j+1, desca, one )
370
371
372
373 CALL psgemv( '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 psgemv( '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 psgemv( '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 psgemv( '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 psgemv( '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 pselget(
'Columnwise',
' ', tau, taup, i, 1, desctp )
390 CALL psscal( 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 psgemv( '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 psgemv( '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 pselset( a, i, j-1, desca, alpha )
417 END IF
418
419
420
421 CALL pslarfg( n-k+1, alpha, i, j, a, i, j+1, desca,
422 $ desca( m_ ), taup )
423 CALL pselset( d, i, 1, descd, alpha )
424 CALL pselset( a, i, j, desca, one )
425
426
427
428 CALL psgemv( '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 psgemv( '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 psgemv( '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 psgemv( '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 psgemv( '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 pselget(
'Columnwise',
' ', tau, taup, i, 1, desctp )
445 CALL psscal( m-k, tau, x, ix+k, jx+k-1, descx, 1 )
446
447
448
449 CALL psgemv( '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 psgemv( '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 pselset( a, i, j, desca, alpha )
456
457
458
459 CALL pslarfg( m-k, alpha, i+1, j, a,
min( i+2, m+ia-1 ),
460 $ j, desca, 1, tauq )
461 CALL pselset( e, 1, j, desce, alpha )
462 CALL pselset( a, i+1, j, desca, one )
463
464
465
466 CALL psgemv( '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 psgemv( '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 psgemv( '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 psgemv( '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 psgemv( '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 pselget(
'Rowwise',
' ', tau, tauq, 1, j, desctq )
483 CALL psscal( n-k, tau, work( ipy ), 1, jwy, descwy,
484 $ descwy( m_ ) )
485 CALL pscopy( 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 pselget(scope, top, alpha, a, ia, ja, desca)
subroutine pselset(a, ia, ja, desca, alpha)
subroutine pslarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)