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)