1 SUBROUTINE pdgeqpf( M, N, A, IA, JA, DESCA, IPIV, TAU, WORK,
10 INTEGER IA, JA, INFO, LWORK, M, N
13 INTEGER DESCA( * ), IPIV( * )
14 DOUBLE PRECISION A( * ), TAU( * ), WORK( * )
185 INTEGER BLOCK_CYCLIC_2D, CSRC_, CTXT_, DLEN_, DTYPE_,
186 $ lld_, mb_, m_, nb_, n_, rsrc_
187 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
188 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
189 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
190 DOUBLE PRECISION ONE, ZERO
191 parameter( one = 1.0d+0, zero = 0.0d+0 )
195 INTEGER I, IACOL, IAROW, ICOFF, ICTXT, ICURROW,
196 $ icurcol, ii, iia, ioffa, ipn, ipcol, ipw,
197 $ iroff, itemp, j, jb, jj, jja, jjpvt, jn, kb,
198 $ k, kk, kstart, kstep, lda, ll, lwmin, mn, mp,
199 $ mycol, myrow, npcol, nprow, nq, nq0, pvt
200 DOUBLE PRECISION AJJ, ALPHA, TEMP, TEMP2, TOL3Z
203 INTEGER DESCN( DLEN_ ), IDUM1( 1 ), IDUM2( 1 )
207 $ dgebr2d, dgebs2d, dgerv2d,
208 $ dgesd2d, dlarfg, dswap, igerv2d,
214 INTEGER ICEIL, INDXG2P, NUMROC
215 EXTERNAL iceil, indxg2p, numroc
216 DOUBLE PRECISION DLAMCH
220 INTRINSIC abs, dble, idint,
max,
min, mod, sqrt
226 ictxt = desca( ctxt_ )
227 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
232 IF( nprow.EQ.-1 )
THEN
235 CALL chk1mat( m, 1, n, 2, ia, ja, desca, 6, info )
237 iroff = mod( ia-1, desca( mb_ ) )
238 icoff = mod( ja-1, desca( nb_ ) )
239 iarow = indxg2p( ia, desca( mb_ ), myrow, desca( rsrc_ ),
241 iacol = indxg2p( ja, desca( nb_ ), mycol, desca( csrc_ ),
243 mp = numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
244 nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
245 nq0 = numroc( ja+n-1, desca( nb_ ), mycol, desca( csrc_ ),
247 lwmin =
max( 3, mp + nq ) + nq0 + nq
249 work( 1 ) = dble( lwmin )
250 lquery = ( lwork.EQ.-1 )
251 IF( lwork.LT.lwmin .AND. .NOT.lquery )
254 IF( lwork.EQ.-1 )
THEN
260 CALL pchk1mat( m, 1, n, 2, ia, ja, desca, 6, 1, idum1, idum2,
265 CALL pxerbla( ictxt,
'PDGEQPF', -info )
267 ELSE IF( lquery )
THEN
273 IF( m.EQ.0 .OR. n.EQ.0 )
276 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
283 tol3z = sqrt( dlamch(
'Epsilon') )
288 jn =
min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
289 kstep = npcol * desca( nb_ )
291 IF( mycol.EQ.iacol )
THEN
296 DO 10 ll = jja, jja+jb-1
297 ipiv( ll ) = ja + ll - jja
299 kstart = jn + kstep - desca( nb_ )
303 DO 30 kk = jja+jb, jja+nq-1, desca( nb_ )
304 kb =
min( jja+nq-kk, desca( nb_ ) )
305 DO 20 ll = kk, kk+kb-1
306 ipiv( ll ) = kstart+ll-kk+1
308 kstart = kstart + kstep
311 kstart = jn + ( mod( mycol-iacol+npcol, npcol )-1 )*
313 DO 50 kk = jja, jja+nq-1, desca( nb_ )
314 kb =
min( jja+nq-kk, desca( nb_ ) )
315 DO 40 ll = kk, kk+kb-1
316 ipiv( ll ) = kstart+ll-kk+1
318 kstart = kstart + kstep
324 CALL descset( descn, 1, desca( n_ ), 1, desca( nb_ ), myrow,
325 $ desca( csrc_ ), ictxt, 1 )
330 IF( mycol.EQ.iacol )
THEN
332 CALL pdnrm2( m, work( jj+kk ), a, ia, ja+kk, desca, 1 )
333 work( nq+jj+kk ) = work( jj+kk )
337 icurcol = mod( iacol+1, npcol )
341 DO 80 j = jn+1, ja+n-1, desca( nb_ )
342 jb =
min( ja+n-j, desca( nb_ ) )
344 IF( mycol.EQ.icurcol )
THEN
346 CALL pdnrm2( m, work( jj+kk ), a, ia, j+kk, desca, 1 )
347 work( nq+jj+kk ) = work( jj+kk )
351 icurcol = mod( icurcol+1, npcol )
356 DO 120 j = ja, ja+mn-1
359 CALL infog1l( j, desca( nb_ ), npcol, mycol, desca( csrc_ ),
363 CALL pdamax( k, temp, pvt, work( ipn ), 1, j, descn,
369 CALL infog1l( pvt, desca( nb_ ), npcol, mycol,
370 $ desca( csrc_ ), jjpvt, ipcol )
371 IF( icurcol.EQ.ipcol )
THEN
372 IF( mycol.EQ.icurcol )
THEN
373 CALL dswap( mp, a( iia+(jj-1)*lda ), 1,
374 $ a( iia+(jjpvt-1)*lda ), 1 )
375 itemp = ipiv( jjpvt )
376 ipiv( jjpvt ) = ipiv( jj )
378 work( ipn+jjpvt-1 ) = work( ipn+jj-1 )
379 work( ipn+nq+jjpvt-1 ) = work( ipn+nq+jj-1 )
382 IF( mycol.EQ.icurcol )
THEN
384 CALL dgesd2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
386 work( ipw ) = dble( ipiv( jj ) )
387 work( ipw+1 ) = work( ipn + jj - 1 )
388 work( ipw+2 ) = work( ipn + nq + jj - 1 )
389 CALL dgesd2d( ictxt, 3, 1, work( ipw ), 3, myrow,
392 CALL dgerv2d( ictxt, mp, 1, a( iia+(jj-1)*lda ), lda,
394 CALL igerv2d( ictxt, 1, 1, ipiv( jj ), 1, myrow,
397 ELSE IF( mycol.EQ.ipcol )
THEN
399 CALL dgesd2d( ictxt, mp, 1, a( iia+(jjpvt-1)*lda ),
400 $ lda, myrow, icurcol )
401 CALL igesd2d( ictxt, 1, 1, ipiv( jjpvt ), 1, myrow,
404 CALL dgerv2d( ictxt, mp, 1, a( iia+(jjpvt-1)*lda ),
405 $ lda, myrow, icurcol )
406 CALL dgerv2d( ictxt, 3, 1, work( ipw ), 3, myrow,
408 ipiv( jjpvt ) = idint( work( ipw ) )
409 work( ipn+jjpvt-1 ) = work( ipw+1 )
410 work( ipn+nq+jjpvt-1 ) = work( ipw+2 )
420 CALL infog1l( i, desca( mb_ ), nprow, myrow, desca( rsrc_ ),
422 IF( desca( m_ ).EQ.1 )
THEN
423 IF( myrow.EQ.icurrow )
THEN
424 IF( mycol.EQ.icurcol )
THEN
425 ioffa = ii+(jj-1)*desca( lld_ )
427 CALL dlarfg( 1, ajj, a( ioffa ), 1, tau( jj ) )
429 alpha = one - tau( jj )
430 CALL dgebs2d( ictxt,
'Rowwise',
' ', 1, 1, alpha,
432 CALL dscal( nq-jj, alpha, a( ioffa+desca( lld_ ) ),
435 CALL dgebs2d( ictxt,
'Columnwise',
' ', 1, 1,
440 CALL dgebr2d( ictxt,
'Rowwise',
' ', 1, 1, alpha,
441 $ 1, icurrow, icurcol )
442 CALL dscal( nq-jj+1, alpha, a( i ), desca( lld_ ) )
445 ELSE IF( mycol.EQ.icurcol )
THEN
446 CALL dgebr2d( ictxt,
'Columnwise',
' ', 1, 1, tau( jj ),
447 $ 1, icurrow, icurcol )
452 CALL pdlarfg( m-j+ja, ajj, i, j, a,
min( i+1, ia+m-1 ), j,
454 IF( j.LT.ja+n-1 )
THEN
458 CALL pdelset( a, i, j, desca, one )
459 CALL pdlarf(
'Left', m-j+ja, ja+n-1-j, a, i, j, desca,
460 $ 1, tau, a, i, j+1, desca, work( ipw ) )
462 CALL pdelset( a, i, j, desca, ajj )
468 IF( mycol.EQ.icurcol )
470 IF( mod( j, desca( nb_ ) ).EQ.0 )
471 $ icurcol = mod( icurcol+1, npcol )
472 IF( (jja+nq-jj).GT.0 )
THEN
473 IF( myrow.EQ.icurrow )
THEN
474 CALL dgebs2d( ictxt,
'Columnwise',
' ', 1, jja+nq-jj,
475 $ a( ii+(
min( jja+nq-1, jj )-1 )*lda ),
477 CALL dcopy( jja+nq-jj, a( ii+(
min( jja+nq-1, jj )
478 $ -1)*lda ), lda, work( ipw+
min( jja+nq-1,
481 CALL dgebr2d( ictxt,
'Columnwise',
' ', jja+nq-jj, 1,
482 $ work( ipw+
min( jja+nq-1, jj )-1 ),
483 $
max( 1, nq ), icurrow, mycol )
487 jn =
min( iceil( j+1, desca( nb_ ) ) * desca( nb_ ),
489 IF( mycol.EQ.icurcol )
THEN
490 DO 90 ll = jj-1, jj + jn - j - 2
491 IF( work( ipn+ll ).NE.zero )
THEN
492 temp = abs( work( ipw+ll ) ) / work( ipn+ll )
493 temp =
max( zero, ( one+temp )*( one-temp ) )
494 temp2 = temp*( work( ipn+ll ) / work( ipn+nq+ll ) )**2
495 IF( temp2.LE.tol3z )
THEN
496 IF( ia+m-1.GT.i )
THEN
497 CALL pdnrm2( ia+m-i-1, work( ipn+ll ), a, i+1,
498 $ j+ll-jj+2, desca, 1 )
499 work( ipn+nq+ll ) = work( ipn+ll )
501 work( ipn+ll ) = zero
502 work( ipn+nq+ll ) = zero
505 work( ipn+ll ) = work( ipn+ll ) * sqrt( temp )
511 icurcol = mod( icurcol+1, npcol )
513 DO 110 k = jn+1, ja+n-1, desca( nb_ )
514 kb =
min( ja+n-k, desca( nb_ ) )
516 IF( mycol.EQ.icurcol )
THEN
517 DO 100 ll = jj-1, jj+kb-2
518 IF( work( ipn+ll ).NE.zero )
THEN
519 temp = abs( work( ipw+ll ) ) / work( ipn+ll )
520 temp =
max( zero, ( one+temp )*( one-temp ) )
522 $ ( work( ipn+ll ) / work( ipn+nq+ll ) )**2
523 IF( temp2.LE.tol3z )
THEN
524 IF( ia+m-1.GT.i )
THEN
525 CALL pdnrm2( ia+m-i-1, work( ipn+ll ), a,
526 $ i+1, k+ll-jj+1, desca, 1 )
527 work( ipn+nq+ll ) = work( ipn+ll )
529 work( ipn+ll ) = zero
530 work( ipn+nq+ll ) = zero
533 work( ipn+ll ) = work( ipn+ll ) * sqrt( temp )
539 icurcol = mod( icurcol+1, npcol )
545 work( 1 ) = dble( lwmin )