3
4
5
6
7
8
9
10 INTEGER IA, JA, INFO, LRWORK, LWORK, M, N
11
12
13 INTEGER DESCA( * ), IPIV( * )
14 REAL RWORK( * )
15 COMPLEX A( * ), TAU( * ), 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 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
203 $ LLD_, MB_, M_, NB_, N_, RSRC_
204 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
205 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
206 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
207 REAL ONE, ZERO
208 parameter( one = 1.0e+0, zero = 0.0e+0 )
209
210
211 LOGICAL LQUERY
212 INTEGER I, IACOL, IAROW, ICOFF, ICTXT, ICURROW,
213 $ ICURCOL, II, IIA, IOFFA, IPCOL, IROFF, ITEMP,
214 $ J, JB, JJ, JJA, JJPVT, JN, KB, K, KK, KSTART,
215 $ KSTEP, LDA, LL, LRWMIN, LWMIN, MN, MP, MYCOL,
216 $ MYROW, NPCOL, NPROW, NQ, NQ0, PVT
217 REAL TEMP, TEMP2, TOL3Z
218 COMPLEX AJJ, ALPHA
219
220
221 INTEGER DESCN( DLEN_ ), IDUM1( 2 ), IDUM2( 2 )
222
223
224 EXTERNAL blacs_gridinfo, ccopy, cgebr2d, cgebs2d,
225 $ cgerv2d, cgesd2d,
chk1mat, clarfg,
229
230
231 INTEGER ICEIL, INDXG2P, NUMROC
233 REAL SLAMCH
235
236
237 INTRINSIC abs,
cmplx, conjg, ifix,
max,
min, mod, sqrt
238
239
240
241
242
243 ictxt = desca( ctxt_ )
244 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
245
246
247
248 info = 0
249 IF( nprow.EQ.-1 ) THEN
250 info = -(600+ctxt_)
251 ELSE
252 CALL chk1mat( m, 1, n, 2, ia, ja, desca, 6, info )
253 IF( info.EQ.0 ) THEN
254 iroff = mod( ia-1, desca( mb_ ) )
255 icoff = mod( ja-1, desca( nb_ ) )
256 iarow =
indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
257 $ nprow )
258 iacol =
indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
259 $ npcol )
260 mp =
numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
261 nq =
numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
262 nq0 =
numroc( ja+n-1, desca( nb_ ), mycol, desca( csrc_ ),
263 $ npcol )
264 lwmin =
max( 3, mp + nq )
265 lrwmin = nq0 + nq
266
267 work( 1 ) =
cmplx( real( lwmin ) )
268 rwork( 1 ) = real( lrwmin )
269 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
270 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
271 info = -10
272 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
273 info = -12
274 END IF
275 END IF
276 IF( lwork.EQ.-1 ) THEN
277 idum1( 1 ) = -1
278 ELSE
279 idum1( 1 ) = 1
280 END IF
281 idum2( 1 ) = 10
282 IF( lrwork.EQ.-1 ) THEN
283 idum1( 2 ) = -1
284 ELSE
285 idum1( 2 ) = 1
286 END IF
287 idum2( 2 ) = 12
288 CALL pchk1mat( m, 1, n, 2, ia, ja, desca, 6, 2, idum1, idum2,
289 $ info )
290 END IF
291
292 IF( info.NE.0 ) THEN
293 CALL pxerbla( ictxt,
'PCGEQPF', -info )
294 RETURN
295 ELSE IF( lquery ) THEN
296 RETURN
297 END IF
298
299
300
301 IF( m.EQ.0 .OR. n.EQ.0 )
302 $ RETURN
303
304 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
305 $ iarow, iacol )
306 IF( myrow.EQ.iarow )
307 $ mp = mp - iroff
308 IF( mycol.EQ.iacol )
309 $ nq = nq - icoff
311 tol3z = sqrt(
slamch(
'Epsilon') )
312
313
314
315 lda = desca( lld_ )
316 jn =
min(
iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
317 kstep = npcol * desca( nb_ )
318
319 IF( mycol.EQ.iacol ) THEN
320
321
322
323 jb = jn - ja + 1
324 DO 10 ll = jja, jja+jb-1
325 ipiv( ll ) = ja + ll - jja
326 10 CONTINUE
327 kstart = jn + kstep - desca( nb_ )
328
329
330
331 DO 30 kk = jja+jb, jja+nq-1, desca( nb_ )
332 kb =
min( jja+nq-kk, desca( nb_ ) )
333 DO 20 ll = kk, kk+kb-1
334 ipiv( ll ) = kstart+ll-kk+1
335 20 CONTINUE
336 kstart = kstart + kstep
337 30 CONTINUE
338 ELSE
339 kstart = jn + ( mod( mycol-iacol+npcol, npcol )-1 )*
340 $ desca( nb_ )
341 DO 50 kk = jja, jja+nq-1, desca( nb_ )
342 kb =
min( jja+nq-kk, desca( nb_ ) )
343 DO 40 ll = kk, kk+kb-1
344 ipiv( ll ) = kstart+ll-kk+1
345 40 CONTINUE
346 kstart = kstart + kstep
347 50 CONTINUE
348 END IF
349
350
351
352 CALL descset( descn, 1, desca( n_ ), 1, desca( nb_ ), myrow,
353 $ desca( csrc_ ), ictxt, 1 )
354
355 jj = jja
356 IF( mycol.EQ.iacol ) THEN
357 DO 60 kk = 0, jb-1
358 CALL pscnrm2( m, rwork( jj+kk ), a, ia, ja+kk, desca, 1 )
359 rwork( nq+jj+kk ) = rwork( jj+kk )
360 60 CONTINUE
361 jj = jj + jb
362 END IF
363 icurcol = mod( iacol+1, npcol )
364
365
366
367 DO 80 j = jn+1, ja+n-1, desca( nb_ )
368 jb =
min( ja+n-j, desca( nb_ ) )
369
370 IF( mycol.EQ.icurcol ) THEN
371 DO 70 kk = 0, jb-1
372 CALL pscnrm2( m, rwork( jj+kk ), a, ia, j+kk, desca, 1 )
373 rwork( nq+jj+kk ) = rwork( jj+kk )
374 70 CONTINUE
375 jj = jj + jb
376 END IF
377 icurcol = mod( icurcol+1, npcol )
378 80 CONTINUE
379
380
381
382 DO 120 j = ja, ja+mn-1
383 i = ia + j - ja
384
385 CALL infog1l( j, desca( nb_ ), npcol, mycol, desca( csrc_ ),
386 $ jj, icurcol )
387 k = ja + n - j
388 IF( k.GT.1 ) THEN
389 CALL psamax( k, temp, pvt, rwork, 1, j, descn,
390 $ descn( m_ ) )
391 ELSE
392 pvt = j
393 END IF
394 IF( j.NE.pvt ) THEN
395 CALL infog1l( pvt, desca( nb_ ), npcol, mycol,
396 $ desca( csrc_ ), jjpvt, ipcol )
397 IF( icurcol.EQ.ipcol ) THEN
398 IF( mycol.EQ.icurcol ) THEN
399 CALL cswap( mp, a( iia+(jj-1)*lda ), 1,
400 $ a( iia+(jjpvt-1)*lda ), 1 )
401 itemp = ipiv( jjpvt )
402 ipiv( jjpvt ) = ipiv( jj )
403 ipiv( jj ) = itemp
404 rwork( jjpvt ) = rwork( jj )
405 rwork( nq+jjpvt ) = rwork( nq+jj )
406 END IF
407 ELSE
408 IF( mycol.EQ.icurcol ) THEN
409
410 CALL cgesd2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
411 $ myrow, ipcol )
412 work( 1 ) =
cmplx( real( ipiv( jj ) ) )
413 work( 2 ) =
cmplx( rwork( jj ) )
414 work( 3 ) =
cmplx( rwork( jj + nq ) )
415 CALL cgesd2d( ictxt, 3, 1, work, 3, myrow, ipcol )
416
417 CALL cgerv2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
418 $ myrow, ipcol )
419 CALL igerv2d( ictxt, 1, 1, ipiv( jj ), 1, myrow,
420 $ ipcol )
421
422 ELSE IF( mycol.EQ.ipcol ) THEN
423
424 CALL cgesd2d( ictxt, mp, 1, a( iia+(jjpvt-1)*lda ),
425 $ lda, myrow, icurcol )
426 CALL igesd2d( ictxt, 1, 1, ipiv( jjpvt ), 1, myrow,
427 $ icurcol )
428
429 CALL cgerv2d( ictxt, mp, 1, a( iia+(jjpvt-1)*lda ),
430 $ lda, myrow, icurcol )
431 CALL cgerv2d( ictxt, 3, 1, work, 3, myrow, icurcol )
432 ipiv( jjpvt ) = ifix( real( work( 1 ) ) )
433 rwork( jjpvt ) = real( work( 2 ) )
434 rwork( jjpvt+nq ) = real( work( 3 ) )
435
436 END IF
437
438 END IF
439
440 END IF
441
442
443
444 CALL infog1l( i, desca( mb_ ), nprow, myrow, desca( rsrc_ ),
445 $ ii, icurrow )
446 IF( desca( m_ ).EQ.1 ) THEN
447 IF( myrow.EQ.icurrow ) THEN
448 IF( mycol.EQ.icurcol ) THEN
449 ioffa = ii+(jj-1)*desca( lld_ )
450 ajj = a( ioffa )
451 CALL clarfg( 1, ajj, a( ioffa ), 1, tau( jj ) )
452 IF( n.GT.1 ) THEN
453 alpha =
cmplx( one ) - conjg( tau( jj ) )
454 CALL cgebs2d( ictxt, 'Rowwise', ' ', 1, 1, alpha,
455 $ 1 )
456 CALL cscal( nq-jj, alpha, a( ioffa+desca( lld_ ) ),
457 $ desca( lld_ ) )
458 END IF
459 CALL cgebs2d( ictxt, 'Columnwise', ' ', 1, 1,
460 $ tau( jj ), 1 )
461 a( ioffa ) = ajj
462 ELSE
463 IF( n.GT.1 ) THEN
464 CALL cgebr2d( ictxt, 'Rowwise', ' ', 1, 1, alpha,
465 $ 1, icurrow, icurcol )
466 CALL cscal( nq-jj+1, alpha, a( i ), desca( lld_ ) )
467 END IF
468 END IF
469 ELSE IF( mycol.EQ.icurcol ) THEN
470 CALL cgebr2d( ictxt, 'Columnwise', ' ', 1, 1, tau( jj ),
471 $ 1, icurrow, icurcol )
472 END IF
473
474 ELSE
475
476 CALL pclarfg( m-j+ja, ajj, i, j, a,
min( i+1, ia+m-1 ), j,
477 $ desca, 1, tau )
478 IF( j.LT.ja+n-1 ) THEN
479
480
481
483 CALL pclarfc(
'Left', m-j+ja, ja+n-1-j, a, i, j, desca,
484 $ 1, tau, a, i, j+1, desca, work )
485 END IF
486 CALL pcelset( a, i, j, desca, ajj )
487
488 END IF
489
490
491
492 IF( mycol.EQ.icurcol )
493 $ jj = jj + 1
494 IF( mod( j, desca( nb_ ) ).EQ.0 )
495 $ icurcol = mod( icurcol+1, npcol )
496 IF( (jja+nq-jj).GT.0 ) THEN
497 IF( myrow.EQ.icurrow ) THEN
498 CALL cgebs2d( ictxt, 'Columnwise', ' ', 1, jja+nq-jj,
499 $ a( ii+(
min( jja+nq-1, jj )-1 )*lda ),
500 $ lda )
501 CALL ccopy( jja+nq-jj, a( ii+(
min( jja+nq-1, jj )
502 $ -1)*lda ), lda, work(
min( jja+nq-1, jj ) ),
503 $ 1 )
504 ELSE
505 CALL cgebr2d( ictxt, 'Columnwise', ' ', jja+nq-jj, 1,
506 $ work(
min( jja+nq-1, jj ) ),
max( 1, nq ),
507 $ icurrow, mycol )
508 END IF
509 END IF
510
511 jn =
min(
iceil( j+1, desca( nb_ ) ) * desca( nb_ ),
512 $ ja + n - 1 )
513 IF( mycol.EQ.icurcol ) THEN
514 DO 90 ll = jj, jj + jn - j - 1
515 IF( rwork( ll ).NE.zero ) THEN
516 temp = abs( work( ll ) ) / rwork( ll )
517 temp =
max( zero, ( one+temp )*( one-temp ) )
518 temp2 = temp * ( rwork( ll ) / rwork( nq+ll ) )**2
519 IF( temp2.LE.tol3z ) THEN
520 IF( ia+m-1.GT.i ) THEN
521 CALL pscnrm2( ia+m-i-1, rwork( ll ), a,
522 $ i+1, j+ll-jj+1, desca, 1 )
523 rwork( nq+ll ) = rwork( ll )
524 ELSE
525 rwork( ll ) = zero
526 rwork( nq+ll ) = zero
527 END IF
528 ELSE
529 rwork( ll ) = rwork( ll ) * sqrt( temp )
530 END IF
531 END IF
532 90 CONTINUE
533 jj = jj + jn - j
534 END IF
535 icurcol = mod( icurcol+1, npcol )
536
537 DO 110 k = jn+1, ja+n-1, desca( nb_ )
538 kb =
min( ja+n-k, desca( nb_ ) )
539
540 IF( mycol.EQ.icurcol ) THEN
541 DO 100 ll = jj, jj+kb-1
542 IF( rwork(ll).NE.zero ) THEN
543 temp = abs( work( ll ) ) / rwork( ll )
544 temp =
max( zero, ( one+temp )*( one-temp ) )
545 temp2 = temp * ( rwork( ll ) / rwork( nq+ll ) )**2
546 IF( temp2.LE.tol3z ) THEN
547 IF( ia+m-1.GT.i ) THEN
548 CALL pscnrm2( ia+m-i-1, rwork( ll ), a,
549 $ i+1, k+ll-jj, desca, 1 )
550 rwork( nq+ll ) = rwork( ll )
551 ELSE
552 rwork( ll ) = zero
553 rwork( nq+ll ) = zero
554 END IF
555 ELSE
556 rwork( ll ) = rwork( ll ) * sqrt( temp )
557 END IF
558 END IF
559 100 CONTINUE
560 jj = jj + kb
561 END IF
562 icurcol = mod( icurcol+1, npcol )
563
564 110 CONTINUE
565
566 120 CONTINUE
567
568 work( 1 ) =
cmplx( real( lwmin ) )
569 rwork( 1 ) = real( lrwmin )
570
571 RETURN
572
573
574
subroutine chk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, info)
subroutine descset(desc, m, n, mb, nb, irsrc, icsrc, ictxt, lld)
integer function iceil(inum, idenom)
integer function indxg2p(indxglob, nb, iproc, isrcproc, nprocs)
subroutine infog1l(gindx, nb, nprocs, myroc, isrcproc, lindx, rocsrc)
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
integer function numroc(n, nb, iproc, isrcproc, nprocs)
subroutine pcelset(a, ia, ja, desca, alpha)
subroutine pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
subroutine pclarfc(side, m, n, v, iv, jv, descv, incv, tau, c, ic, jc, descc, work)
subroutine pclarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)
subroutine pxerbla(ictxt, srname, info)