3
4
5
6
7
8
9
10 CHARACTER DIRECT, STOREV
11 INTEGER IV, JV, K, N
12
13
14 INTEGER DESCV( * )
15 DOUBLE PRECISION TAU( * ), T( * ), V( * ), 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 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
174 $ LLD_, MB_, M_, NB_, N_, RSRC_
175 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
176 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
177 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
178 DOUBLE PRECISION ONE, ZERO
179 parameter( one = 1.0d+0, zero = 0.0d+0 )
180
181
182 LOGICAL FORWARD
183 INTEGER ICOFF, ICTXT, II, IIV, IROFF, IVCOL, IVROW,
184 $ ITMP0, ITMP1, IW, JJ, JJV, LDV, MICOL, MIROW,
185 $ MYCOL, MYROW, NP, NPCOL, NPROW, NQ
186 DOUBLE PRECISION VII
187
188
189 EXTERNAL blacs_gridinfo, dcopy, dgemv, dgsum2d,
191
192
193 LOGICAL LSAME
194 INTEGER INDXG2P, NUMROC
196
197
198 INTRINSIC mod
199
200
201
202
203
204 IF( n.LE.0 .OR. k.LE.0 )
205 $ RETURN
206
207 ictxt = descv( ctxt_ )
208 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
209
210 forward =
lsame( direct,
'F' )
211 CALL infog2l( iv, jv, descv, nprow, npcol, myrow, mycol,
212 $ iiv, jjv, ivrow, ivcol )
213
214 IF(
lsame( storev,
'C' ) .AND. mycol.EQ.ivcol )
THEN
215
216 iw = 1
217 ldv = descv( lld_ )
218 iroff = mod( iv-1, descv( mb_ ) )
219
220 IF( forward ) THEN
221
222
223
224 np =
numroc( n+iroff, descv( mb_ ), myrow, ivrow, nprow )
225 IF( myrow.EQ.ivrow ) THEN
226 np = np - iroff
227 ii = iiv + 1
228 ELSE
229 ii = iiv
230 END IF
231 IF( iroff+1.EQ.descv( mb_ ) ) THEN
232 mirow = mod( ivrow+1, nprow )
233 ELSE
234 mirow = ivrow
235 END IF
236 itmp0 = 0
237
238 DO 10 jj = jjv+1, jjv+k-1
239
240 IF( myrow.EQ.mirow ) THEN
241 vii = v( ii+(jj-1)*ldv )
242 v( ii+(jj-1)*ldv ) = one
243 END IF
244
245
246
247
248 itmp0 = itmp0 + 1
249 IF( np-ii+iiv.GT.0 ) THEN
250 CALL dgemv( 'Transpose', np-ii+iiv, itmp0,
251 $ -tau( jj ), v( ii+(jjv-1)*ldv ), ldv,
252 $ v( ii+(jj-1)*ldv ), 1, zero,
253 $ work( iw ), 1 )
254 ELSE
255 CALL dlaset( 'All', itmp0, 1, zero, zero, work( iw ),
256 $ itmp0 )
257 END IF
258
259 iw = iw + itmp0
260 IF( myrow.EQ.mirow ) THEN
261 v( ii+(jj-1)*ldv ) = vii
262 ii = ii + 1
263 END IF
264
265 IF( mod( iv+itmp0, descv( mb_ ) ).EQ.0 )
266 $ mirow = mod( mirow+1, nprow )
267
268 10 CONTINUE
269
270 CALL dgsum2d( ictxt, 'Columnwise', ' ', iw-1, 1, work, iw-1,
271 $ ivrow, mycol )
272
273 IF( myrow.EQ.ivrow ) THEN
274
275 iw = 1
276 itmp0 = 0
277 itmp1 = 1
278
279 t( itmp1 ) = tau( jjv )
280
281 DO 20 jj = jjv+1, jjv+k-1
282
283
284
285 itmp0 = itmp0 + 1
286 itmp1 = itmp1 + descv( nb_ )
287 CALL dcopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
288 iw = iw + itmp0
289
290 CALL dtrmv( 'Upper', 'No transpose', 'Non-unit',
291 $ itmp0, t, descv( nb_ ), t( itmp1 ), 1 )
292 t(itmp1+itmp0) = tau( jj )
293
294 20 CONTINUE
295
296 END IF
297
298 ELSE
299
300
301
302 np =
numroc( n+iroff-1, descv( mb_ ), myrow, ivrow, nprow )
303 IF( myrow.EQ.ivrow )
304 $ np = np - iroff
305 mirow =
indxg2p( iv+n-2, descv( mb_ ), myrow,
306 $ descv( rsrc_ ), nprow )
307 ii = iiv + np - 1
308 itmp0 = 0
309
310 DO 30 jj = jjv+k-2, jjv, -1
311
312 IF( myrow.EQ.mirow ) THEN
313 vii = v( ii+(jj-1)*ldv )
314 v( ii+(jj-1)*ldv ) = one
315 END IF
316
317
318
319
320 itmp0 = itmp0 + 1
321 IF( ii-iiv+1.GT.0 ) THEN
322 CALL dgemv( 'Transpose', ii-iiv+1, itmp0, -tau( jj ),
323 $ v( iiv+jj*ldv ), ldv,
324 $ v( iiv+(jj-1)*ldv ), 1, zero,
325 $ work( iw ), 1 )
326 ELSE
327 CALL dlaset( 'All', itmp0, 1, zero, zero, work( iw ),
328 $ itmp0 )
329 END IF
330
331 iw = iw + itmp0
332 IF( myrow.EQ.mirow ) THEN
333 v( ii+(jj-1)*ldv ) = vii
334 ii = ii - 1
335 END IF
336
337 IF( mod( iv+n-itmp0-2, descv(mb_) ).EQ.0 )
338 $ mirow = mod( mirow+nprow-1, nprow )
339
340 30 CONTINUE
341
342 CALL dgsum2d( ictxt, 'Columnwise', ' ', iw-1, 1, work, iw-1,
343 $ ivrow, mycol )
344
345 IF( myrow.EQ.ivrow ) THEN
346
347 iw = 1
348 itmp0 = 0
349 itmp1 = k + 1 + (k-1) * descv( nb_ )
350
351 t( itmp1-1 ) = tau( jjv+k-1 )
352
353 DO 40 jj = jjv+k-2, jjv, -1
354
355
356
357 itmp0 = itmp0 + 1
358 itmp1 = itmp1 - descv( nb_ ) - 1
359 CALL dcopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
360 iw = iw + itmp0
361
362 CALL dtrmv( 'Lower', 'No transpose', 'Non-unit',
363 $ itmp0, t( itmp1+descv( nb_ ) ),
364 $ descv( nb_ ), t( itmp1 ), 1 )
365 t( itmp1-1 ) = tau( jj )
366
367 40 CONTINUE
368
369 END IF
370
371 END IF
372
373 ELSE IF(
lsame( storev,
'R' ) .AND. myrow.EQ.ivrow )
THEN
374
375 iw = 1
376 ldv = descv( lld_ )
377 icoff = mod( jv-1, descv( nb_ ) )
378
379 IF( forward ) THEN
380
381
382
383 nq =
numroc( n+icoff, descv( nb_ ), mycol, ivcol, npcol )
384 IF( mycol.EQ.ivcol ) THEN
385 nq = nq - icoff
386 jj = jjv + 1
387 ELSE
388 jj = jjv
389 END IF
390 IF( icoff+1.EQ.descv( nb_ ) ) THEN
391 micol = mod( ivcol+1, npcol )
392 ELSE
393 micol = ivcol
394 END IF
395 itmp0 = 0
396
397 DO 50 ii = iiv+1, iiv+k-1
398
399 IF( mycol.EQ.micol ) THEN
400 vii = v( ii+(jj-1)*ldv )
401 v( ii+(jj-1)*ldv ) = one
402 END IF
403
404
405
406
407 itmp0 = itmp0 + 1
408 IF( nq-jj+jjv.GT.0 ) THEN
409 CALL dgemv( 'No transpose', itmp0, nq-jj+jjv,
410 $ -tau(ii), v( iiv+(jj-1)*ldv ), ldv,
411 $ v( ii+(jj-1)*ldv ), ldv, zero,
412 $ work( iw ), 1 )
413 ELSE
414 CALL dlaset( 'All', itmp0, 1, zero, zero,
415 $ work( iw ), itmp0 )
416 END IF
417
418 iw = iw + itmp0
419 IF( mycol.EQ.micol ) THEN
420 v( ii+(jj-1)*ldv ) = vii
421 jj = jj + 1
422 END IF
423
424 IF( mod( jv+itmp0, descv( nb_ ) ).EQ.0 )
425 $ micol = mod( micol+1, npcol )
426
427 50 CONTINUE
428
429 CALL dgsum2d( ictxt, 'Rowwise', ' ', iw-1, 1, work, iw-1,
430 $ myrow, ivcol )
431
432 IF( mycol.EQ.ivcol ) THEN
433
434 iw = 1
435 itmp0 = 0
436 itmp1 = 1
437
438 t( itmp1 ) = tau( iiv )
439
440 DO 60 ii = iiv+1, iiv+k-1
441
442
443
444 itmp0 = itmp0 + 1
445 itmp1 = itmp1 + descv( mb_ )
446 CALL dcopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
447 iw = iw + itmp0
448
449 CALL dtrmv( 'Upper', 'No transpose', 'Non-unit',
450 $ itmp0, t, descv( mb_ ), t( itmp1 ), 1 )
451 t( itmp1+itmp0 ) = tau( ii )
452
453 60 CONTINUE
454
455 END IF
456
457 ELSE
458
459
460
461 nq =
numroc( n+icoff-1, descv( nb_ ), mycol, ivcol, npcol )
462 IF( mycol.EQ.ivcol )
463 $ nq = nq - icoff
464 micol =
indxg2p( jv+n-2, descv( nb_ ), mycol,
465 $ descv( csrc_ ), npcol )
466 jj = jjv + nq - 1
467 itmp0 = 0
468
469 DO 70 ii = iiv+k-2, iiv, -1
470
471 IF( mycol.EQ.micol ) THEN
472 vii = v( ii+(jj-1)*ldv )
473 v( ii+(jj-1)*ldv ) = one
474 END IF
475
476
477
478
479 itmp0 = itmp0 + 1
480 IF( jj-jjv+1.GT.0 ) THEN
481 CALL dgemv( 'No transpose', itmp0, jj-jjv+1,
482 $ -tau( ii ), v( ii+1+(jjv-1)*ldv ), ldv,
483 $ v( ii+(jjv-1)*ldv ), ldv, zero,
484 $ work( iw ), 1 )
485 ELSE
486 CALL dlaset( 'All', itmp0, 1, zero, zero,
487 $ work( iw ), itmp0 )
488 END IF
489
490 iw = iw + itmp0
491 IF( mycol.EQ.micol ) THEN
492 v( ii+(jj-1)*ldv ) = vii
493 jj = jj - 1
494 END IF
495
496 IF( mod( jv+n-itmp0-2, descv( nb_ ) ).EQ.0 )
497 $ micol = mod( micol+npcol-1, npcol )
498
499 70 CONTINUE
500
501 CALL dgsum2d( ictxt, 'Rowwise', ' ', iw-1, 1, work, iw-1,
502 $ myrow, ivcol )
503
504 IF( mycol.EQ.ivcol ) THEN
505
506 iw = 1
507 itmp0 = 0
508 itmp1 = k + 1 + (k-1) * descv( mb_ )
509
510 t( itmp1-1 ) = tau( iiv+k-1 )
511
512 DO 80 ii = iiv+k-2, iiv, -1
513
514
515
516 itmp0 = itmp0 + 1
517 itmp1 = itmp1 - descv( mb_ ) - 1
518 CALL dcopy( itmp0, work( iw ), 1, t( itmp1 ), 1 )
519 iw = iw + itmp0
520
521 CALL dtrmv( 'Lower', 'No transpose', 'Non-unit',
522 $ itmp0, t( itmp1+descv( mb_ ) ),
523 $ descv( mb_ ), t( itmp1 ), 1 )
524 t( itmp1-1 ) = tau( ii )
525
526 80 CONTINUE
527
528 END IF
529
530 END IF
531
532 END IF
533
534 RETURN
535
536
537
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
integer function numroc(n, nb, iproc, isrcproc, nprocs)