3
4
5
6
7
8
9
10 CHARACTER HOWMNY, SIDE
11 INTEGER INFO, M, MM, N
12
13
14 LOGICAL SELECT( * )
15 INTEGER DESCT( * ), DESCVL( * ), DESCVR( * )
16 DOUBLE PRECISION RWORK( * )
17 COMPLEX*16 T( * ), VL( * ), VR( * ), WORK( * )
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 DOUBLE PRECISION ZERO, ONE
205 parameter( zero = 0.0d+0, one = 1.0d+0 )
206 COMPLEX*16 CZERO, CONE
207 parameter( czero = ( 0.0d+0, 0.0d+0 ),
208 $ cone = ( 1.0d+0, 0.0d+0 ) )
209 INTEGER BLOCK_CYCLIC_2D, DLEN_, DTYPE_, CTXT_, M_, N_,
210 $ MB_, NB_, RSRC_, CSRC_, LLD_
211 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
212 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
213 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
214
215
216 LOGICAL ALLV, BOTHV, LEFTV, OVER, RIGHTV, SOMEV
217 INTEGER CONTXT, CSRC, I, ICOL, II, IROW, IS, ITMP1,
218 $ ITMP2, J, K, KI, LDT, LDVL, LDVR, LDW, MB,
219 $ MYCOL, MYROW, NB, NPCOL, NPROW, RSRC
220 REAL SELF
221 DOUBLE PRECISION OVFL, REMAXD, SCALE, SMLNUM, ULP, UNFL
222 COMPLEX*16 CDUM, REMAXC, SHIFT
223
224
225 INTEGER DESCW( DLEN_ )
226 DOUBLE PRECISION SMIN( 1 )
227
228
229 LOGICAL LSAME
230 DOUBLE PRECISION PDLAMCH
232
233
234 EXTERNAL blacs_gridinfo,
descinit, dgsum2d, igamn2d,
237 $ zgsum2d
238
239
240 INTRINSIC abs, dble, dcmplx, dconjg, dimag,
max
241
242
243 DOUBLE PRECISION CABS1
244
245
246 cabs1( cdum ) = abs( dble( cdum ) ) + abs( dimag( cdum ) )
247
248
249
250
251 IF( block_cyclic_2d*csrc_*ctxt_*dlen_*dtype_*lld_*mb_*m_*nb_*n_*
252 $ rsrc_.LT.0 )RETURN
253
254 contxt = desct( ctxt_ )
255 rsrc = desct( rsrc_ )
256 csrc = desct( csrc_ )
257 mb = desct( mb_ )
258 nb = desct( nb_ )
259 ldt = desct( lld_ )
260 ldw = ldt
261 ldvr = descvr( lld_ )
262 ldvl = descvl( lld_ )
263
264 CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
265 self = myrow*npcol + mycol
266
267
268
269 bothv =
lsame( side,
'B' )
270 rightv =
lsame( side,
'R' ) .OR. bothv
271 leftv =
lsame( side,
'L' ) .OR. bothv
272
273 allv =
lsame( howmny,
'A' )
274 over =
lsame( howmny,
'B' ) .OR.
lsame( howmny,
'O' )
275 somev =
lsame( howmny,
'S' )
276
277
278
279
280 IF( somev ) THEN
281 m = 0
282 DO 10 j = 1, n
283 IF( SELECT( j ) )
284 $ m = m + 1
285 10 CONTINUE
286 ELSE
287 m = n
288 END IF
289
290 info = 0
291 IF( .NOT.rightv .AND. .NOT.leftv ) THEN
292 info = -1
293 ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
294 info = -2
295 ELSE IF( n.LT.0 ) THEN
296 info = -4
297 ELSE IF( mm.LT.m ) THEN
298 info = -11
299 END IF
300 CALL igamn2d( contxt, 'ALL', ' ', 1, 1, info, 1, itmp1, itmp2, -1,
301 $ -1, -1 )
302 IF( info.LT.0 ) THEN
303 CALL pxerbla( contxt,
'PZTREVC', -info )
304 RETURN
305 END IF
306
307
308
309 IF( n.EQ.0 )
310 $ RETURN
311
312
313
314 unfl =
pdlamch( contxt,
'Safe minimum' )
315 ovfl = one / unfl
316 CALL pdlabad( contxt, unfl, ovfl )
317 ulp =
pdlamch( contxt,
'Precision' )
318 smlnum = unfl*( n / ulp )
319
320
321
322 DO 20 i = 1, n
323 CALL infog2l( i, i, desct, nprow, npcol, myrow, mycol, irow,
324 $ icol, itmp1, itmp2 )
325 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
326 work( ldw+irow ) = t( ( icol-1 )*ldt+irow )
327 END IF
328 20 CONTINUE
329
330
331
332
333
334 rwork( 1 ) = zero
335 DO 30 j = 2, n
336 CALL pdzasum( j-1, rwork( j ), t, 1, j, desct, 1 )
337 30 CONTINUE
338
339
340 CALL dgsum2d( contxt, 'Row', ' ', n, 1, rwork, n, -1, -1 )
341
342 IF( rightv ) THEN
343
344
345
346
347
348 CALL descinit( descw, n, 1, nb, 1, rsrc, csrc, contxt, ldw,
349 $ info )
350
351 is = m
352 DO 70 ki = n, 1, -1
353
354 IF( somev ) THEN
355 IF( .NOT.SELECT( ki ) )
356 $ GO TO 70
357 END IF
358
359 smin( 1 ) = zero
360 shift = czero
361 CALL infog2l( ki, ki, desct, nprow, npcol, myrow, mycol,
362 $ irow, icol, itmp1, itmp2 )
363 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
364 shift = t( ( icol-1 )*ldt+irow )
365 smin( 1 ) =
max( ulp*( cabs1( shift ) ), smlnum )
366 END IF
367 CALL dgsum2d( contxt, 'ALL', ' ', 1, 1, smin, 1, -1, -1 )
368 CALL zgsum2d( contxt, 'ALL', ' ', 1, 1, shift, 1, -1, -1 )
369
370 CALL infog2l( 1, 1, descw, nprow, npcol, myrow, mycol, irow,
371 $ icol, itmp1, itmp2 )
372 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
373 work( 1 ) = cone
374 END IF
375
376
377
378
379 IF( ki.GT.1 ) THEN
380 CALL pzcopy( ki-1, t, 1, ki, desct, 1, work, 1, 1, descw,
381 $ 1 )
382 END IF
383 DO 40 k = 1, ki - 1
384 CALL infog2l( k, 1, descw, nprow, npcol, myrow, mycol,
385 $ irow, icol, itmp1, itmp2 )
386 IF( myrow.EQ.itmp1 .AND. mycol.EQ.itmp2 ) THEN
387 work( irow ) = -work( irow )
388 END IF
389 40 CONTINUE
390
391
392
393
394 DO 50 k = 1, ki - 1
395 CALL infog2l( k, k, desct, nprow, npcol, myrow, mycol,
396 $ irow, icol, itmp1, itmp2 )
397 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
398 t( ( icol-1 )*ldt+irow ) = t( ( icol-1 )*ldt+irow ) -
399 $ shift
400 IF( cabs1( t( ( icol-1 )*ldt+irow ) ).LT.smin( 1 ) )
401 $ THEN
402 t( ( icol-1 )*ldt+irow ) = dcmplx( smin( 1 ) )
403 END IF
404 END IF
405 50 CONTINUE
406
407 IF( ki.GT.1 ) THEN
408 CALL pzlattrs(
'Upper',
'No transpose',
'Non-unit',
'Y',
409 $ ki-1, t, 1, 1, desct, work, 1, 1, descw,
410 $ scale, rwork, info )
411 CALL infog2l( ki, 1, descw, nprow, npcol, myrow, mycol,
412 $ irow, icol, itmp1, itmp2 )
413 IF( myrow.EQ.itmp1 .AND. mycol.EQ.itmp2 ) THEN
414 work( irow ) = dcmplx( scale )
415 END IF
416 END IF
417
418
419
420 IF( .NOT.over ) THEN
421 CALL pzcopy( ki, work, 1, 1, descw, 1, vr, 1, is, descvr,
422 $ 1 )
423
424 CALL pzamax( ki, remaxc, ii, vr, 1, is, descvr, 1 )
425 remaxd = one /
max( cabs1( remaxc ), unfl )
426 CALL pzdscal( ki, remaxd, vr, 1, is, descvr, 1 )
427
428 CALL pzlaset(
' ', n-ki, 1, czero, czero, vr, ki+1, is,
429 $ descvr )
430 ELSE
431 IF( ki.GT.1 )
432 $ CALL pzgemv( 'N', n, ki-1, cone, vr, 1, 1, descvr,
433 $ work, 1, 1, descw, 1, dcmplx( scale ),
434 $ vr, 1, ki, descvr, 1 )
435
436 CALL pzamax( n, remaxc, ii, vr, 1, ki, descvr, 1 )
437 remaxd = one /
max( cabs1( remaxc ), unfl )
438 CALL pzdscal( n, remaxd, vr, 1, ki, descvr, 1 )
439 END IF
440
441
442
443 DO 60 k = 1, ki - 1
444 CALL infog2l( k, k, desct, nprow, npcol, myrow, mycol,
445 $ irow, icol, itmp1, itmp2 )
446 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
447 t( ( icol-1 )*ldt+irow ) = work( ldw+irow )
448 END IF
449 60 CONTINUE
450
451 is = is - 1
452 70 CONTINUE
453 END IF
454
455 IF( leftv ) THEN
456
457
458
459
460
461 CALL descinit( descw, n, 1, mb, 1, rsrc, csrc, contxt, ldw,
462 $ info )
463
464 is = 1
465 DO 110 ki = 1, n
466
467 IF( somev ) THEN
468 IF( .NOT.SELECT( ki ) )
469 $ GO TO 110
470 END IF
471
472 smin( 1 ) = zero
473 shift = czero
474 CALL infog2l( ki, ki, desct, nprow, npcol, myrow, mycol,
475 $ irow, icol, itmp1, itmp2 )
476 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
477 shift = t( ( icol-1 )*ldt+irow )
478 smin( 1 ) =
max( ulp*( cabs1( shift ) ), smlnum )
479 END IF
480 CALL dgsum2d( contxt, 'ALL', ' ', 1, 1, smin, 1, -1, -1 )
481 CALL zgsum2d( contxt, 'ALL', ' ', 1, 1, shift, 1, -1, -1 )
482
483 CALL infog2l( n, 1, descw, nprow, npcol, myrow, mycol, irow,
484 $ icol, itmp1, itmp2 )
485 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
486 work( irow ) = cone
487 END IF
488
489
490
491 IF( ki.LT.n ) THEN
492 CALL pzcopy( n-ki, t, ki, ki+1, desct, n, work, ki+1, 1,
493 $ descw, 1 )
494 END IF
495 DO 80 k = ki + 1, n
496 CALL infog2l( k, 1, descw, nprow, npcol, myrow, mycol,
497 $ irow, icol, itmp1, itmp2 )
498 IF( myrow.EQ.itmp1 .AND. mycol.EQ.itmp2 ) THEN
499 work( irow ) = -dconjg( work( irow ) )
500 END IF
501 80 CONTINUE
502
503
504
505
506 DO 90 k = ki + 1, n
507 CALL infog2l( k, k, desct, nprow, npcol, myrow, mycol,
508 $ irow, icol, itmp1, itmp2 )
509 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
510 t( ( icol-1 )*ldt+irow ) = t( ( icol-1 )*ldt+irow ) -
511 $ shift
512 IF( cabs1( t( ( icol-1 )*ldt+irow ) ).LT.smin( 1 ) )
513 $ t( ( icol-1 )*ldt+irow ) = dcmplx( smin( 1 ) )
514 END IF
515 90 CONTINUE
516
517 IF( ki.LT.n ) THEN
518 CALL pzlattrs(
'Upper',
'Conjugate transpose',
'Nonunit',
519 $ 'Y', n-ki, t, ki+1, ki+1, desct, work,
520 $ ki+1, 1, descw, scale, rwork, info )
521 CALL infog2l( ki, 1, descw, nprow, npcol, myrow, mycol,
522 $ irow, icol, itmp1, itmp2 )
523 IF( myrow.EQ.itmp1 .AND. mycol.EQ.itmp2 ) THEN
524 work( irow ) = dcmplx( scale )
525 END IF
526 END IF
527
528
529
530 IF( .NOT.over ) THEN
531 CALL pzcopy( n-ki+1, work, ki, 1, descw, 1, vl, ki, is,
532 $ descvl, 1 )
533
534 CALL pzamax( n-ki+1, remaxc, ii, vl, ki, is, descvl, 1 )
535 remaxd = one /
max( cabs1( remaxc ), unfl )
536 CALL pzdscal( n-ki+1, remaxd, vl, ki, is, descvl, 1 )
537
538 CALL pzlaset(
' ', ki-1, 1, czero, czero, vl, 1, is,
539 $ descvl )
540 ELSE
541 IF( ki.LT.n )
542 $ CALL pzgemv( 'N', n, n-ki, cone, vl, 1, ki+1, descvl,
543 $ work, ki+1, 1, descw, 1, dcmplx( scale ),
544 $ vl, 1, ki, descvl, 1 )
545
546 CALL pzamax( n, remaxc, ii, vl, 1, ki, descvl, 1 )
547 remaxd = one /
max( cabs1( remaxc ), unfl )
548 CALL pzdscal( n, remaxd, vl, 1, ki, descvl, 1 )
549 END IF
550
551
552
553 DO 100 k = ki + 1, n
554 CALL infog2l( k, k, desct, nprow, npcol, myrow, mycol,
555 $ irow, icol, itmp1, itmp2 )
556 IF( ( myrow.EQ.itmp1 ) .AND. ( mycol.EQ.itmp2 ) ) THEN
557 t( ( icol-1 )*ldt+irow ) = work( ldw+irow )
558 END IF
559 100 CONTINUE
560
561 is = is + 1
562 110 CONTINUE
563 END IF
564
565 RETURN
566
567
568
subroutine descinit(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld, info)
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
double precision function pdlamch(ictxt, cmach)
subroutine pdlabad(ictxt, small, large)
subroutine pxerbla(ictxt, srname, info)
subroutine pzlaset(uplo, m, n, alpha, beta, a, ia, ja, desca)
subroutine pzlattrs(uplo, trans, diag, normin, n, a, ia, ja, desca, x, ix, jx, descx, scale, cnorm, info)