1 SUBROUTINE pcgeqpf( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
2 $ LWORK, RWORK, LRWORK, INFO )
10 INTEGER IA, JA, INFO, LRWORK, LWORK, M, N
13 INTEGER DESCA( * ), IPIV( * )
15 COMPLEX A( * ), TAU( * ), WORK( * )
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 )
208 parameter( one = 1.0e+0, zero = 0.0e+0 )
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
221 INTEGER DESCN( DLEN_ ), IDUM1( 2 ), IDUM2( 2 )
224 EXTERNAL blacs_gridinfo, ccopy, cgebr2d, cgebs2d,
225 $ cgerv2d, cgesd2d,
chk1mat, clarfg,
231 INTEGER ICEIL, INDXG2P, NUMROC
232 EXTERNAL iceil, indxg2p, numroc
237 INTRINSIC abs,
cmplx, conjg, ifix,
max,
min, mod, sqrt
243 ictxt = desca( ctxt_ )
244 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
249 IF( nprow.EQ.-1 )
THEN
252 CALL chk1mat( m, 1, n, 2, ia, ja, desca, 6, info )
254 iroff = mod( ia-1, desca( mb_ ) )
255 icoff = mod( ja-1, desca( nb_ ) )
256 iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
258 iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
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_ ),
264 lwmin =
max( 3, mp + nq )
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
272 ELSE IF( lrwork.LT.lrwmin .AND. .NOT.lquery )
THEN
276 IF( lwork.EQ.-1 )
THEN
282 IF( lrwork.EQ.-1 )
THEN
288 CALL pchk1mat( m, 1, n, 2, ia, ja, desca, 6, 2, idum1, idum2,
293 CALL pxerbla( ictxt,
'PCGEQPF', -info )
295 ELSE IF( lquery )
THEN
301 IF( m.EQ.0 .OR. n.EQ.0 )
304 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
311 tol3z = sqrt( slamch(
'Epsilon') )
316 jn =
min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
317 kstep = npcol * desca( nb_ )
319 IF( mycol.EQ.iacol )
THEN
324 DO 10 ll = jja, jja+jb-1
325 ipiv( ll ) = ja + ll - jja
327 kstart = jn + kstep - desca( nb_ )
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
336 kstart = kstart + kstep
339 kstart = jn + ( mod( mycol-iacol+npcol, npcol )-1 )*
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
346 kstart = kstart + kstep
352 CALL descset( descn, 1, desca( n_ ), 1, desca( nb_ ), myrow,
353 $ desca( csrc_ ), ictxt, 1 )
356 IF( mycol.EQ.iacol )
THEN
358 CALL pscnrm2( m, rwork( jj+kk ), a, ia, ja+kk, desca, 1 )
359 rwork( nq+jj+kk ) = rwork( jj+kk )
363 icurcol = mod( iacol+1, npcol )
367 DO 80 j = jn+1, ja+n-1, desca( nb_ )
368 jb =
min( ja+n-j, desca( nb_ ) )
370 IF( mycol.EQ.icurcol )
THEN
372 CALL pscnrm2( m, rwork( jj+kk ), a, ia, j+kk, desca, 1 )
373 rwork( nq+jj+kk ) = rwork( jj+kk )
377 icurcol = mod( icurcol+1, npcol )
382 DO 120 j = ja, ja+mn-1
385 CALL infog1l( j, desca( nb_ ), npcol, mycol, desca( csrc_ ),
389 CALL psamax( k, temp, pvt, rwork, 1, j, descn,
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 )
404 rwork( jjpvt ) = rwork( jj )
405 rwork( nq+jjpvt ) = rwork( nq+jj )
408 IF( mycol.EQ.icurcol )
THEN
410 CALL cgesd2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
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 )
417 CALL cgerv2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
419 CALL igerv2d( ictxt, 1, 1, ipiv( jj ), 1, myrow,
422 ELSE IF( mycol.EQ.ipcol )
THEN
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,
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 ) )
444 CALL infog1l( i, desca( mb_ ), nprow, myrow, desca( rsrc_ ),
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_ )
451 CALL clarfg( 1, ajj, a( ioffa ), 1, tau( jj ) )
453 alpha =
cmplx( one ) - conjg( tau( jj ) )
454 CALL cgebs2d( ictxt,
'Rowwise',
' ', 1, 1, alpha,
456 CALL cscal( nq-jj, alpha, a( ioffa+desca( lld_ ) ),
459 CALL cgebs2d( ictxt,
'Columnwise',
' ', 1, 1,
464 CALL cgebr2d( ictxt,
'Rowwise',
' ', 1, 1, alpha,
465 $ 1, icurrow, icurcol )
466 CALL cscal( nq-jj+1, alpha, a( i ), desca( lld_ ) )
469 ELSE IF( mycol.EQ.icurcol )
THEN
470 CALL cgebr2d( ictxt,
'Columnwise',
' ', 1, 1, tau( jj ),
471 $ 1, icurrow, icurcol )
476 CALL pclarfg( m-j+ja, ajj, i, j, a,
min( i+1, ia+m-1 ), j,
478 IF( j.LT.ja+n-1 )
THEN
483 CALL pclarfc(
'Left', m-j+ja, ja+n-1-j, a, i, j, desca,
484 $ 1, tau, a, i, j+1, desca, work )
486 CALL pcelset( a, i, j, desca, ajj )
492 IF( mycol.EQ.icurcol )
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 ),
501 CALL ccopy( jja+nq-jj, a( ii+(
min( jja+nq-1, jj )
502 $ -1)*lda ), lda, work(
min( jja+nq-1, jj ) ),
505 CALL cgebr2d( ictxt,
'Columnwise',
' ', jja+nq-jj, 1,
506 $ work(
min( jja+nq-1, jj ) ),
max( 1, nq ),
511 jn =
min( iceil( j+1, desca( nb_ ) ) * desca( nb_ ),
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 )
526 rwork( nq+ll ) = zero
529 rwork( ll ) = rwork( ll ) * sqrt( temp )
535 icurcol = mod( icurcol+1, npcol )
537 DO 110 k = jn+1, ja+n-1, desca( nb_ )
538 kb =
min( ja+n-k, desca( nb_ ) )
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 )
553 rwork( nq+ll ) = zero
556 rwork( ll ) = rwork( ll ) * sqrt( temp )
562 icurcol = mod( icurcol+1, npcol )
568 work( 1 ) =
cmplx( real( lwmin ) )
569 rwork( 1 ) = real( lrwmin )