3
4
5
6
7
8
9
10 INTEGER I, L, LWORK, M
11 COMPLEX H33, H43H34, H44
12
13
14 INTEGER DESCA( * )
15 COMPLEX A( * ), BUF( * )
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 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
161 $ LLD_, MB_, M_, NB_, N_, RSRC_
162 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
163 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
164 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
165
166
167 INTEGER CONTXT, DOWN, HBL, IBUF1, IBUF2, IBUF3, IBUF4,
168 $ IBUF5, ICOL1, II, IRCV1, IRCV2, IRCV3, IRCV4,
169 $ IRCV5, IROW1, ISRC, ISTR1, ISTR2, ISTR3, ISTR4,
170 $ ISTR5, JJ, JSRC, LDA, LEFT, MODKM1, MYCOL,
171 $ MYROW, NPCOL, NPROW, NUM, RIGHT, UP
172 REAL S, TST1, ULP
173 COMPLEX CDUM, H00, H10, H11, H12, H21, H22, H33S, H44S,
174 $ V1, V2, V3
175
176
177 INTEGER ILCM
178 REAL PSLAMCH
180
181
183 $ cgerv2d, cgesd2d
184
185
186 INTRINSIC abs, real, aimag, mod
187
188
189 REAL CABS1
190
191
192 cabs1( cdum ) = abs( real( cdum ) ) + abs( aimag( cdum ) )
193
194
195
196 hbl = desca( mb_ )
197 contxt = desca( ctxt_ )
198 lda = desca( lld_ )
199 ulp =
pslamch( contxt,
'PRECISION' )
200 CALL blacs_gridinfo( contxt, nprow, npcol, myrow, mycol )
201 left = mod( mycol+npcol-1, npcol )
202 right = mod( mycol+1, npcol )
203 up = mod( myrow+nprow-1, nprow )
204 down = mod( myrow+1, nprow )
205 num = nprow*npcol
206
207
208
209
210
211
212
213 istr1 = 0
214 istr2 = ( ( i-l-1 ) / hbl )
215 IF( istr2*hbl.LT.( i-l-1 ) )
216 $ istr2 = istr2 + 1
217 ii = istr2 /
ilcm( nprow, npcol )
218 IF( ii*
ilcm( nprow, npcol ).LT.istr2 )
THEN
219 istr2 = ii + 1
220 ELSE
221 istr2 = ii
222 END IF
223 IF( lwork.LT.7*istr2 ) THEN
224 CALL pxerbla( contxt,
'PCLACONSB', 10 )
225 RETURN
226 END IF
227 istr3 = 3*istr2
228 istr4 = istr3 + istr2
229 istr5 = istr3 + istr3
230 CALL infog2l( i-2, i-2, desca, nprow, npcol, myrow, mycol, irow1,
231 $ icol1, ii, jj )
232 modkm1 = mod( i-3+hbl, hbl )
233
234
235
236
237
238 ibuf1 = 0
239 ibuf2 = 0
240 ibuf3 = 0
241 ibuf4 = 0
242 ibuf5 = 0
243 ircv1 = 0
244 ircv2 = 0
245 ircv3 = 0
246 ircv4 = 0
247 ircv5 = 0
248 DO 10 m = i - 2, l, -1
249 IF( ( modkm1.EQ.0 ) .AND. ( down.EQ.ii ) .AND.
250 $ ( right.EQ.jj ) .AND. ( m.GT.l ) ) THEN
251
252
253
254 IF( ( down.NE.myrow ) .OR. ( right.NE.mycol ) ) THEN
255 CALL infog2l( m-1, m-1, desca, nprow, npcol, myrow,
256 $ mycol, irow1, icol1, isrc, jsrc )
257 ibuf1 = ibuf1 + 1
258 buf( istr1+ibuf1 ) = a( ( icol1-1 )*lda+irow1 )
259 END IF
260 END IF
261 IF( ( modkm1.EQ.0 ) .AND. ( myrow.EQ.ii ) .AND.
262 $ ( right.EQ.jj ) .AND. ( m.GT.l ) ) THEN
263
264
265
266 IF( npcol.GT.1 ) THEN
267 CALL infog2l( m, m-1, desca, nprow, npcol, myrow, mycol,
268 $ irow1, icol1, isrc, jsrc )
269 ibuf5 = ibuf5 + 1
270 buf( istr5+ibuf5 ) = a( ( icol1-1 )*lda+irow1 )
271 END IF
272 END IF
273 IF( ( modkm1.EQ.hbl-1 ) .AND. ( up.EQ.ii ) .AND.
274 $ ( mycol.EQ.jj ) ) THEN
275
276
277
278 IF( nprow.GT.1 ) THEN
279 CALL infog2l( m+1, m, desca, nprow, npcol, myrow, mycol,
280 $ irow1, icol1, isrc, jsrc )
281 ibuf2 = ibuf2 + 1
282 buf( istr2+ibuf2 ) = a( ( icol1-1 )*lda+irow1 )
283 END IF
284 END IF
285 IF( ( modkm1.EQ.hbl-1 ) .AND. ( myrow.EQ.ii ) .AND.
286 $ ( left.EQ.jj ) ) THEN
287
288
289
290 IF( npcol.GT.1 ) THEN
291 CALL infog2l( m, m+1, desca, nprow, npcol, myrow, mycol,
292 $ irow1, icol1, isrc, jsrc )
293 ibuf3 = ibuf3 + 1
294 buf( istr3+ibuf3 ) = a( ( icol1-1 )*lda+irow1 )
295 END IF
296 END IF
297 IF( ( modkm1.EQ.hbl-1 ) .AND. ( up.EQ.ii ) .AND.
298 $ ( left.EQ.jj ) ) THEN
299
300
301
302
303 IF( ( up.NE.myrow ) .OR. ( left.NE.mycol ) ) THEN
304 CALL infog2l( m+1, m+1, desca, nprow, npcol, myrow,
305 $ mycol, irow1, icol1, isrc, jsrc )
306 ibuf4 = ibuf4 + 2
307 buf( istr4+ibuf4-1 ) = a( ( icol1-1 )*lda+irow1 )
308 buf( istr4+ibuf4 ) = a( ( icol1-1 )*lda+irow1+1 )
309 END IF
310 END IF
311 IF( ( modkm1.EQ.hbl-2 ) .AND. ( up.EQ.ii ) .AND.
312 $ ( mycol.EQ.jj ) ) THEN
313
314
315
316 IF( nprow.GT.1 ) THEN
317 CALL infog2l( m+2, m+1, desca, nprow, npcol, myrow,
318 $ mycol, irow1, icol1, isrc, jsrc )
319 ibuf2 = ibuf2 + 1
320 buf( istr2+ibuf2 ) = a( ( icol1-1 )*lda+irow1 )
321 END IF
322 END IF
323
324
325
326 IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
327 IF( ( modkm1.EQ.0 ) .AND. ( m.GT.l ) .AND.
328 $ ( ( nprow.GT.1 ) .OR. ( npcol.GT.1 ) ) ) THEN
329
330
331
332 ircv1 = ircv1 + 1
333 END IF
334 IF( ( modkm1.EQ.0 ) .AND. ( npcol.GT.1 ) .AND. ( m.GT.l ) )
335 $ THEN
336
337
338
339 ircv5 = ircv5 + 1
340 END IF
341 IF( ( modkm1.EQ.hbl-1 ) .AND. ( nprow.GT.1 ) ) THEN
342
343
344
345 ircv2 = ircv2 + 1
346 END IF
347 IF( ( modkm1.EQ.hbl-1 ) .AND. ( npcol.GT.1 ) ) THEN
348
349
350
351 ircv3 = ircv3 + 1
352 END IF
353 IF( ( modkm1.EQ.hbl-1 ) .AND.
354 $ ( ( nprow.GT.1 ) .OR. ( npcol.GT.1 ) ) ) THEN
355
356
357
358 ircv4 = ircv4 + 2
359 END IF
360 IF( ( modkm1.EQ.hbl-2 ) .AND. ( nprow.GT.1 ) ) THEN
361
362
363
364 ircv2 = ircv2 + 1
365 END IF
366 END IF
367
368
369
370 IF( modkm1.EQ.0 ) THEN
371 ii = ii - 1
372 jj = jj - 1
373 IF( ii.LT.0 )
374 $ ii = nprow - 1
375 IF( jj.LT.0 )
376 $ jj = npcol - 1
377 END IF
378 modkm1 = modkm1 - 1
379 IF( modkm1.LT.0 )
380 $ modkm1 = hbl - 1
381 10 CONTINUE
382
383
384
385
386 IF( ibuf1.GT.0 ) THEN
387 CALL cgesd2d( contxt, ibuf1, 1, buf( istr1+1 ), ibuf1, down,
388 $ right )
389 END IF
390 IF( ibuf2.GT.0 ) THEN
391 CALL cgesd2d( contxt, ibuf2, 1, buf( istr2+1 ), ibuf2, up,
392 $ mycol )
393 END IF
394 IF( ibuf3.GT.0 ) THEN
395 CALL cgesd2d( contxt, ibuf3, 1, buf( istr3+1 ), ibuf3, myrow,
396 $ left )
397 END IF
398 IF( ibuf4.GT.0 ) THEN
399 CALL cgesd2d( contxt, ibuf4, 1, buf( istr4+1 ), ibuf4, up,
400 $ left )
401 END IF
402 IF( ibuf5.GT.0 ) THEN
403 CALL cgesd2d( contxt, ibuf5, 1, buf( istr5+1 ), ibuf5, myrow,
404 $ right )
405 END IF
406
407
408
409 IF( ircv1.GT.0 ) THEN
410 CALL cgerv2d( contxt, ircv1, 1, buf( istr1+1 ), ircv1, up,
411 $ left )
412 END IF
413 IF( ircv2.GT.0 ) THEN
414 CALL cgerv2d( contxt, ircv2, 1, buf( istr2+1 ), ircv2, down,
415 $ mycol )
416 END IF
417 IF( ircv3.GT.0 ) THEN
418 CALL cgerv2d( contxt, ircv3, 1, buf( istr3+1 ), ircv3, myrow,
419 $ right )
420 END IF
421 IF( ircv4.GT.0 ) THEN
422 CALL cgerv2d( contxt, ircv4, 1, buf( istr4+1 ), ircv4, down,
423 $ right )
424 END IF
425 IF( ircv5.GT.0 ) THEN
426 CALL cgerv2d( contxt, ircv5, 1, buf( istr5+1 ), ircv5, myrow,
427 $ left )
428 END IF
429
430
431
432 ibuf1 = 0
433 ibuf2 = 0
434 ibuf3 = 0
435 ibuf4 = 0
436 ibuf5 = 0
437 CALL infog2l( i-2, i-2, desca, nprow, npcol, myrow, mycol, irow1,
438 $ icol1, ii, jj )
439 modkm1 = mod( i-3+hbl, hbl )
440 IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) .AND.
441 $ ( modkm1.NE.hbl-1 ) ) THEN
442 CALL infog2l( i-2, i-1, desca, nprow, npcol, myrow, mycol,
443 $ irow1, icol1, isrc, jsrc )
444 END IF
445
446
447
448 DO 20 m = i - 2, l, -1
449
450
451
452
453
454 IF( ( myrow.EQ.ii ) .AND. ( mycol.EQ.jj ) ) THEN
455 IF( modkm1.EQ.0 ) THEN
456 h22 = a( ( icol1-1 )*lda+irow1+1 )
457 h11 = a( ( icol1-2 )*lda+irow1 )
458 v3 = a( ( icol1-1 )*lda+irow1+2 )
459 h21 = a( ( icol1-2 )*lda+irow1+1 )
460 h12 = a( ( icol1-1 )*lda+irow1 )
461 IF( m.GT.l ) THEN
462 IF( num.GT.1 ) THEN
463 ibuf1 = ibuf1 + 1
464 h00 = buf( istr1+ibuf1 )
465 ELSE
466 h00 = a( ( icol1-3 )*lda+irow1-1 )
467 END IF
468 IF( npcol.GT.1 ) THEN
469 ibuf5 = ibuf5 + 1
470 h10 = buf( istr5+ibuf5 )
471 ELSE
472 h10 = a( ( icol1-3 )*lda+irow1 )
473 END IF
474 END IF
475 END IF
476 IF( modkm1.EQ.hbl-1 ) THEN
477 CALL infog2l( m, m, desca, nprow, npcol, myrow, mycol,
478 $ irow1, icol1, isrc, jsrc )
479 h11 = a( ( icol1-1 )*lda+irow1 )
480 IF( num.GT.1 ) THEN
481 ibuf4 = ibuf4 + 2
482 h22 = buf( istr4+ibuf4-1 )
483 v3 = buf( istr4+ibuf4 )
484 ELSE
485 h22 = a( icol1*lda+irow1+1 )
486 v3 = a( ( icol1+1 )*lda+irow1+1 )
487 END IF
488 IF( nprow.GT.1 ) THEN
489 ibuf2 = ibuf2 + 1
490 h21 = buf( istr2+ibuf2 )
491 ELSE
492 h21 = a( ( icol1-1 )*lda+irow1+1 )
493 END IF
494 IF( npcol.GT.1 ) THEN
495 ibuf3 = ibuf3 + 1
496 h12 = buf( istr3+ibuf3 )
497 ELSE
498 h12 = a( icol1*lda+irow1 )
499 END IF
500 IF( m.GT.l ) THEN
501 h00 = a( ( icol1-2 )*lda+irow1-1 )
502 h10 = a( ( icol1-2 )*lda+irow1 )
503 END IF
504
505
506
507 icol1 = icol1 + 1
508 END IF
509 IF( modkm1.EQ.hbl-2 ) THEN
510 h22 = a( ( icol1-1 )*lda+irow1+1 )
511 h11 = a( ( icol1-2 )*lda+irow1 )
512 IF( nprow.GT.1 ) THEN
513 ibuf2 = ibuf2 + 1
514 v3 = buf( istr2+ibuf2 )
515 ELSE
516 v3 = a( ( icol1-1 )*lda+irow1+2 )
517 END IF
518 h21 = a( ( icol1-2 )*lda+irow1+1 )
519 h12 = a( ( icol1-1 )*lda+irow1 )
520 IF( m.GT.l ) THEN
521 h00 = a( ( icol1-3 )*lda+irow1-1 )
522 h10 = a( ( icol1-3 )*lda+irow1 )
523 END IF
524 END IF
525 IF( ( modkm1.LT.hbl-2 ) .AND. ( modkm1.GT.0 ) ) THEN
526 h22 = a( ( icol1-1 )*lda+irow1+1 )
527 h11 = a( ( icol1-2 )*lda+irow1 )
528 v3 = a( ( icol1-1 )*lda+irow1+2 )
529 h21 = a( ( icol1-2 )*lda+irow1+1 )
530 h12 = a( ( icol1-1 )*lda+irow1 )
531 IF( m.GT.l ) THEN
532 h00 = a( ( icol1-3 )*lda+irow1-1 )
533 h10 = a( ( icol1-3 )*lda+irow1 )
534 END IF
535 END IF
536 h44s = h44 - h11
537 h33s = h33 - h11
538 v1 = ( h33s*h44s-h43h34 ) / h21 + h12
539 v2 = h22 - h11 - h33s - h44s
540 s = cabs1( v1 ) + cabs1( v2 ) + cabs1( v3 )
541 v1 = v1 / s
542 v2 = v2 / s
543 v3 = v3 / s
544 IF( m.EQ.l )
545 $ GO TO 30
546 tst1 = cabs1( v1 )*( cabs1( h00 )+cabs1( h11 )+
547 $ cabs1( h22 ) )
548 IF( cabs1( h10 )*( cabs1( v2 )+cabs1( v3 ) ).LE.ulp*tst1 )
549 $ GO TO 30
550
551
552
553 irow1 = irow1 - 1
554 icol1 = icol1 - 1
555 END IF
556 IF( m.EQ.l ) THEN
557
558
559
560 GO TO 30
561 END IF
562
563
564
565 IF( modkm1.EQ.0 ) THEN
566 ii = ii - 1
567 jj = jj - 1
568 IF( ii.LT.0 )
569 $ ii = nprow - 1
570 IF( jj.LT.0 )
571 $ jj = npcol - 1
572 END IF
573 modkm1 = modkm1 - 1
574 IF( modkm1.LT.0 )
575 $ modkm1 = hbl - 1
576 20 CONTINUE
577 30 CONTINUE
578
579 CALL igamx2d( contxt, 'ALL', ' ', 1, 1, m, 1, l, l, -1, -1, -1 )
580
581 RETURN
582
583
584
integer function ilcm(m, n)
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
real function pslamch(ictxt, cmach)
subroutine pxerbla(ictxt, srname, info)