3
4
5
6
7
8
9
10 INTEGER IA, JA, INFO, LRWORK, LWORK, M, N
11
12
13 INTEGER DESCA( * ), IPIV( * )
14 DOUBLE PRECISION RWORK( * )
15 COMPLEX*16 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 DOUBLE PRECISION ONE, ZERO
208 parameter( one = 1.0d+0, zero = 0.0d+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 DOUBLE PRECISION TEMP, TEMP2, TOL3Z
218 COMPLEX*16 AJJ, ALPHA
219
220
221 INTEGER DESCN( DLEN_ ), IDUM1( 2 ), IDUM2( 2 )
222
223
228 $ zgebs2d, zgerv2d, zgesd2d, zlarfg,
229 $ zswap
230
231
232 INTEGER ICEIL, INDXG2P, NUMROC
234 DOUBLE PRECISION DLAMCH
236
237
238 INTRINSIC abs, dcmplx, dconjg, idint,
max,
min, mod, sqrt
239
240
241
242
243
244 ictxt = desca( ctxt_ )
245 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
246
247
248
249 info = 0
250 IF( nprow.EQ.-1 ) THEN
251 info = -(600+ctxt_)
252 ELSE
253 CALL chk1mat( m, 1, n, 2, ia, ja, desca, 6, info )
254 IF( info.EQ.0 ) THEN
255 iroff = mod( ia-1, desca( mb_ ) )
256 icoff = mod( ja-1, desca( nb_ ) )
257 iarow =
indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
258 $ nprow )
259 iacol =
indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
260 $ npcol )
261 mp =
numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
262 nq =
numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
263 nq0 =
numroc( ja+n-1, desca( nb_ ), mycol, desca( csrc_ ),
264 $ npcol )
265 lwmin =
max( 3, mp + nq )
266 lrwmin = nq0 + nq
267
268 work( 1 ) = dcmplx( dble( lwmin ) )
269 rwork( 1 ) = dble( lrwmin )
270 lquery = ( lwork.EQ.-1 .OR. lrwork.EQ.-1 )
271 IF( lwork.LT.lwmin .AND. .NOT.lquery ) THEN
272 info = -10
273 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery ) THEN
274 info = -12
275 END IF
276 END IF
277 IF( lwork.EQ.-1 ) THEN
278 idum1( 1 ) = -1
279 ELSE
280 idum1( 1 ) = 1
281 END IF
282 idum2( 1 ) = 10
283 IF( lrwork.EQ.-1 ) THEN
284 idum1( 2 ) = -1
285 ELSE
286 idum1( 2 ) = 1
287 END IF
288 idum2( 2 ) = 12
289 CALL pchk1mat( m, 1, n, 2, ia, ja, desca, 6, 2, idum1, idum2,
290 $ info )
291 END IF
292
293 IF( info.NE.0 ) THEN
294 CALL pxerbla( ictxt,
'PZGEQPF', -info )
295 RETURN
296 ELSE IF( lquery ) THEN
297 RETURN
298 END IF
299
300
301
302 IF( m.EQ.0 .OR. n.EQ.0 )
303 $ RETURN
304
305 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
306 $ iarow, iacol )
307 IF( myrow.EQ.iarow )
308 $ mp = mp - iroff
309 IF( mycol.EQ.iacol )
310 $ nq = nq - icoff
312 tol3z = sqrt(
dlamch(
'Epsilon') )
313
314
315
316 lda = desca( lld_ )
317 jn =
min(
iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
318 kstep = npcol * desca( nb_ )
319
320 IF( mycol.EQ.iacol ) THEN
321
322
323
324 jb = jn - ja + 1
325 DO 10 ll = jja, jja+jb-1
326 ipiv( ll ) = ja + ll - jja
327 10 CONTINUE
328 kstart = jn + kstep - desca( nb_ )
329
330
331
332 DO 30 kk = jja+jb, jja+nq-1, desca( nb_ )
333 kb =
min( jja+nq-kk, desca( nb_ ) )
334 DO 20 ll = kk, kk+kb-1
335 ipiv( ll ) = kstart+ll-kk+1
336 20 CONTINUE
337 kstart = kstart + kstep
338 30 CONTINUE
339 ELSE
340 kstart = jn + ( mod( mycol-iacol+npcol, npcol )-1 )*
341 $ desca( nb_ )
342 DO 50 kk = jja, jja+nq-1, desca( nb_ )
343 kb =
min( jja+nq-kk, desca( nb_ ) )
344 DO 40 ll = kk, kk+kb-1
345 ipiv( ll ) = kstart+ll-kk+1
346 40 CONTINUE
347 kstart = kstart + kstep
348 50 CONTINUE
349 END IF
350
351
352
353 CALL descset( descn, 1, desca( n_ ), 1, desca( nb_ ), myrow,
354 $ desca( csrc_ ), ictxt, 1 )
355
356 jj = jja
357 IF( mycol.EQ.iacol ) THEN
358 DO 60 kk = 0, jb-1
359 CALL pdznrm2( m, rwork( jj+kk ), a, ia, ja+kk, desca, 1 )
360 rwork( nq+jj+kk ) = rwork( jj+kk )
361 60 CONTINUE
362 jj = jj + jb
363 END IF
364 icurcol = mod( iacol+1, npcol )
365
366
367
368 DO 80 j = jn+1, ja+n-1, desca( nb_ )
369 jb =
min( ja+n-j, desca( nb_ ) )
370
371 IF( mycol.EQ.icurcol ) THEN
372 DO 70 kk = 0, jb-1
373 CALL pdznrm2( m, rwork( jj+kk ), a, ia, j+kk, desca, 1 )
374 rwork( nq+jj+kk ) = rwork( jj+kk )
375 70 CONTINUE
376 jj = jj + jb
377 END IF
378 icurcol = mod( icurcol+1, npcol )
379 80 CONTINUE
380
381
382
383 DO 120 j = ja, ja+mn-1
384 i = ia + j - ja
385
386 CALL infog1l( j, desca( nb_ ), npcol, mycol, desca( csrc_ ),
387 $ jj, icurcol )
388 k = ja + n - j
389 IF( k.GT.1 ) THEN
390 CALL pdamax( k, temp, pvt, rwork, 1, j, descn,
391 $ descn( m_ ) )
392 ELSE
393 pvt = j
394 END IF
395 IF( j.NE.pvt ) THEN
396 CALL infog1l( pvt, desca( nb_ ), npcol, mycol,
397 $ desca( csrc_ ), jjpvt, ipcol )
398 IF( icurcol.EQ.ipcol ) THEN
399 IF( mycol.EQ.icurcol ) THEN
400 CALL zswap( mp, a( iia+(jj-1)*lda ), 1,
401 $ a( iia+(jjpvt-1)*lda ), 1 )
402 itemp = ipiv( jjpvt )
403 ipiv( jjpvt ) = ipiv( jj )
404 ipiv( jj ) = itemp
405 rwork( jjpvt ) = rwork( jj )
406 rwork( nq+jjpvt ) = rwork( nq+jj )
407 END IF
408 ELSE
409 IF( mycol.EQ.icurcol ) THEN
410
411 CALL zgesd2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
412 $ myrow, ipcol )
413 work( 1 ) = dcmplx( dble( ipiv( jj ) ) )
414 work( 2 ) = dcmplx( rwork( jj ) )
415 work( 3 ) = dcmplx( rwork( jj + nq ) )
416 CALL zgesd2d( ictxt, 3, 1, work, 3, myrow, ipcol )
417
418 CALL zgerv2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
419 $ myrow, ipcol )
420 CALL igerv2d( ictxt, 1, 1, ipiv( jj ), 1, myrow,
421 $ ipcol )
422
423 ELSE IF( mycol.EQ.ipcol ) THEN
424
425 CALL zgesd2d( ictxt, mp, 1, a( iia+(jjpvt-1)*lda ),
426 $ lda, myrow, icurcol )
427 CALL igesd2d( ictxt, 1, 1, ipiv( jjpvt ), 1, myrow,
428 $ icurcol )
429
430 CALL zgerv2d( ictxt, mp, 1, a( iia+(jjpvt-1)*lda ),
431 $ lda, myrow, icurcol )
432 CALL zgerv2d( ictxt, 3, 1, work, 3, myrow, icurcol )
433 ipiv( jjpvt ) = idint( dble( work( 1 ) ) )
434 rwork( jjpvt ) = dble( work( 2 ) )
435 rwork( jjpvt+nq ) = dble( work( 3 ) )
436
437 END IF
438
439 END IF
440
441 END IF
442
443
444
445 CALL infog1l( i, desca( mb_ ), nprow, myrow, desca( rsrc_ ),
446 $ ii, icurrow )
447 IF( desca( m_ ).EQ.1 ) THEN
448 IF( myrow.EQ.icurrow ) THEN
449 IF( mycol.EQ.icurcol ) THEN
450 ioffa = ii+(jj-1)*desca( lld_ )
451 ajj = a( ioffa )
452 CALL zlarfg( 1, ajj, a( ioffa ), 1, tau( jj ) )
453 IF( n.GT.1 ) THEN
454 alpha =
cmplx( one ) - dconjg( tau( jj ) )
455 CALL zgebs2d( ictxt, 'Rowwise', ' ', 1, 1, alpha,
456 $ 1 )
457 CALL zscal( nq-jj, alpha, a( ioffa+desca( lld_ ) ),
458 $ desca( lld_ ) )
459 END IF
460 CALL zgebs2d( ictxt, 'Columnwise', ' ', 1, 1,
461 $ tau( jj ), 1 )
462 a( ioffa ) = ajj
463 ELSE
464 IF( n.GT.1 ) THEN
465 CALL zgebr2d( ictxt, 'Rowwise', ' ', 1, 1, alpha,
466 $ 1, icurrow, icurcol )
467 CALL zscal( nq-jj+1, alpha, a( i ), desca( lld_ ) )
468 END IF
469 END IF
470 ELSE IF( mycol.EQ.icurcol ) THEN
471 CALL zgebr2d( ictxt, 'Columnwise', ' ', 1, 1, tau( jj ),
472 $ 1, icurrow, icurcol )
473 END IF
474
475 ELSE
476
477 CALL pzlarfg( m-j+ja, ajj, i, j, a,
min( i+1, ia+m-1 ), j,
478 $ desca, 1, tau )
479 IF( j.LT.ja+n-1 ) THEN
480
481
482
483 CALL pzelset( a, i, j, desca, dcmplx( one ) )
484 CALL pzlarfc(
'Left', m-j+ja, ja+n-1-j, a, i, j, desca,
485 $ 1, tau, a, i, j+1, desca, work )
486 END IF
487 CALL pzelset( a, i, j, desca, ajj )
488
489 END IF
490
491
492
493 IF( mycol.EQ.icurcol )
494 $ jj = jj + 1
495 IF( mod( j, desca( nb_ ) ).EQ.0 )
496 $ icurcol = mod( icurcol+1, npcol )
497 IF( (jja+nq-jj).GT.0 ) THEN
498 IF( myrow.EQ.icurrow ) THEN
499 CALL zgebs2d( ictxt, 'Columnwise', ' ', 1, jja+nq-jj,
500 $ a( ii+(
min( jja+nq-1, jj )-1 )*lda ),
501 $ lda )
502 CALL zcopy( jja+nq-jj, a( ii+(
min( jja+nq-1, jj )
503 $ -1)*lda ), lda, work(
min( jja+nq-1, jj ) ),
504 $ 1 )
505 ELSE
506 CALL zgebr2d( ictxt, 'Columnwise', ' ', jja+nq-jj, 1,
507 $ work(
min( jja+nq-1, jj ) ),
max( 1, nq ),
508 $ icurrow, mycol )
509 END IF
510 END IF
511
512 jn =
min(
iceil( j+1, desca( nb_ ) ) * desca( nb_ ),
513 $ ja + n - 1 )
514 IF( mycol.EQ.icurcol ) THEN
515 DO 90 ll = jj, jj + jn - j - 1
516 IF( rwork( ll ).NE.zero ) THEN
517 temp = abs( work( ll ) ) / rwork( ll )
518 temp =
max( zero, ( one+temp )*( one-temp ) )
519 temp2 = temp * ( rwork( ll ) / rwork( nq+ll ) )**2
520 IF( temp2.LE.tol3z ) THEN
521 IF( ia+m-1.GT.i ) THEN
522 CALL pdznrm2( ia+m-i-1, rwork( ll ), a,
523 $ i+1, j+ll-jj+1, desca, 1 )
524 rwork( nq+ll ) = rwork( ll )
525 ELSE
526 rwork( ll ) = zero
527 rwork( nq+ll ) = zero
528 END IF
529 ELSE
530 rwork( ll ) = rwork( ll ) * sqrt( temp )
531 END IF
532 END IF
533 90 CONTINUE
534 jj = jj + jn - j
535 END IF
536 icurcol = mod( icurcol+1, npcol )
537
538 DO 110 k = jn+1, ja+n-1, desca( nb_ )
539 kb =
min( ja+n-k, desca( nb_ ) )
540
541 IF( mycol.EQ.icurcol ) THEN
542 DO 100 ll = jj, jj+kb-1
543 IF( rwork(ll).NE.zero ) THEN
544 temp = abs( work( ll ) ) / rwork( ll )
545 temp =
max( zero, ( one+temp )*( one-temp ) )
546 temp2 = temp * ( rwork( ll ) / rwork( nq+ll ) )**2
547 IF( temp2.LE.tol3z ) THEN
548 IF( ia+m-1.GT.i ) THEN
549 CALL pdznrm2( ia+m-i-1, rwork( ll ), a,
550 $ i+1, k+ll-jj, desca, 1 )
551 rwork( nq+ll ) = rwork( ll )
552 ELSE
553 rwork( ll ) = zero
554 rwork( nq+ll ) = zero
555 END IF
556 ELSE
557 rwork( ll ) = rwork( ll ) * sqrt( temp )
558 END IF
559 END IF
560 100 CONTINUE
561 jj = jj + kb
562 END IF
563 icurcol = mod( icurcol+1, npcol )
564
565 110 CONTINUE
566
567 120 CONTINUE
568
569 work( 1 ) = dcmplx( dble( lwmin ) )
570 rwork( 1 ) = dble( lrwmin )
571
572 RETURN
573
574
575
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 pchk1mat(ma, mapos0, na, napos0, ia, ja, desca, descapos0, nextra, ex, expos, info)
subroutine pxerbla(ictxt, srname, info)
subroutine pzelset(a, ia, ja, desca, alpha)
subroutine pzlarfc(side, m, n, v, iv, jv, descv, incv, tau, c, ic, jc, descc, work)
subroutine pzlarfg(n, alpha, iax, jax, x, ix, jx, descx, incx, tau)