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