pdtreecomb.f

Go to the documentation of this file.

      SUBROUTINE pdtreecomb( ICTXT, SCOPE, N, MINE, RDEST0, CDEST0,
     $                       SUBPTR )
*
*  -- ScaLAPACK tools routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      CHARACTER          SCOPE
      INTEGER            CDEST0, ICTXT, N, RDEST0
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   MINE( * )
*     ..
*     .. Subroutine Arguments ..
      EXTERNAL           subptr
*     ..
*
*  Purpose
*  =======
*
*  PDTREECOMB does a 1-tree parallel combine operation on scalars,
*  using the subroutine indicated by SUBPTR to perform the required
*  computation.
*
*  Arguments
*  =========
*
*  ICTXT   (global input) INTEGER
*          The BLACS context handle, indicating the global context of
*          the operation. The context itself is global.
*
*  SCOPE   (global input) CHARACTER
*          The scope of the operation:  'Rowwise', 'Columnwise', or
*          'All'.
*
*  N       (global input) INTEGER
*          The number of elements in MINE.  N = 1 for the norm-2
*          computation and 2 for the sum of square.
*
*  MINE    (local input/global output) DOUBLE PRECISION array of
*          dimension at least equal to N. The local data to use in the
*          combine.
*
*  RDEST0  (global input) INTEGER
*          The process row to receive the answer. If RDEST0 = -1,
*          every process in the scope gets the answer.
*
*  CDEST0  (global input) INTEGER
*          The process column to receive the answer. If CDEST0 = -1,
*          every process in the scope gets the answer.
*
*  SUBPTR  (local input) Pointer to the subroutine to call to perform
*          the required combine.
*
*  =====================================================================
*
*     .. Local Scalars ..
      LOGICAL            BCAST, RSCOPE, CSCOPE
      INTEGER            CMSSG, DEST, DIST, HISDIST, I, IAM, MYCOL,
     $                   myrow, mydist, mydist2, np, npcol, nprow,
     $                   rmssg, tcdest, trdest
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   HIS( 2 )
*     ..
*     .. External Subroutines ..
      EXTERNAL           blacs_gridinfo, dgebr2d, dgebs2d,
     $                   dgerv2d, dgesd2d
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           lsame
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          mod
*     ..
*     .. Executable Statements ..
*
        dest = 0
*
*     See if everyone wants the answer (need to broadcast the answer)
*
      bcast = ( ( rdest0.EQ.-1 ).OR.( cdest0.EQ.-1 ) )
      IF( bcast ) THEN
         trdest = 0
         tcdest = 0
      ELSE
         trdest = rdest0
         tcdest = cdest0
      END IF
*
*     Get grid parameters.
*
      CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
*
*     Figure scope-dependant variables, or report illegal scope
*
      rscope = lsame( scope, 'R' )
      cscope = lsame( scope, 'C' )
*
      IF( rscope ) THEN
         IF( bcast ) THEN
            trdest = myrow
         ELSE IF( myrow.NE.trdest ) THEN
            RETURN
         END IF
         np = npcol
         mydist = mod( npcol + mycol - tcdest, npcol )
      ELSE IF( cscope ) THEN
         IF( bcast ) THEN
            tcdest = mycol
         ELSE IF( mycol.NE.tcdest ) THEN
            RETURN
         END IF
         np = nprow
         mydist = mod( nprow + myrow - trdest, nprow )
      ELSE IF( lsame( scope, 'A' ) ) THEN
         np = nprow * npcol
         iam = myrow*npcol + mycol
         dest = trdest*npcol + tcdest
         mydist = mod( np + iam - dest, np )
      ELSE
         RETURN
      END IF
*
      IF( np.LT.2 )
     $   RETURN
*
      mydist2 = mydist
      rmssg = myrow
      cmssg = mycol
      i = 1
*
   10 CONTINUE
*
         IF( mod( mydist, 2 ).NE.0 ) THEN
*
*           If I am process that sends information
*
            dist = i * ( mydist - mod( mydist, 2 ) )
*
*           Figure coordinates of dest of message
*
            IF( rscope ) THEN
               cmssg = mod( tcdest + dist, np )
            ELSE IF( cscope ) THEN
               rmssg = mod( trdest + dist, np )
            ELSE
               cmssg = mod( dest + dist, np )
               rmssg = cmssg / npcol
               cmssg = mod( cmssg, npcol )
            END IF
*
            CALL dgesd2d( ictxt, n, 1, mine, n, rmssg, cmssg )
*
            GO TO 20
*
         ELSE
*
*           If I am a process receiving information, figure coordinates
*           of source of message
*
            dist = mydist2 + i
            IF( rscope ) THEN
               cmssg = mod( tcdest + dist, np )
               hisdist = mod( np + cmssg - tcdest, np )
            ELSE IF( cscope ) THEN
               rmssg = mod( trdest + dist, np )
               hisdist = mod( np + rmssg - trdest, np )
            ELSE
               cmssg = mod( dest + dist, np )
               rmssg = cmssg / npcol
               cmssg = mod( cmssg, npcol )
               hisdist = mod( np + rmssg*npcol+cmssg - dest, np )
            END IF
*
            IF( mydist2.LT.hisdist ) THEN
*
*              If I have anyone sending to me
*
               CALL dgerv2d( ictxt, n, 1, his, n, rmssg, cmssg )
               CALL subptr( mine, his )
*
            END IF
            mydist = mydist / 2
*
         END IF
         i = i * 2
*
      IF( i.LT.np )
     $   GO TO 10
*
   20 CONTINUE
*
      IF( bcast ) THEN
         IF( mydist2.EQ.0 ) THEN
            CALL dgebs2d( ictxt, scope, ' ', n, 1, mine, n )
         ELSE
            CALL dgebr2d( ictxt, scope, ' ', n, 1, mine, n,
     $                    trdest, tcdest )
         END IF
      END IF
*
      RETURN
*
*     End of PDTREECOMB
*
      END
*
      SUBROUTINE dcombamax( V1, V2 )
*
*  -- ScaLAPACK tools routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Array Arguments ..
      DOUBLE PRECISION   V1( 2 ), V2( 2 )
*     ..
*
*  Purpose
*  =======
*
*  DCOMBAMAX finds the element having max. absolute value as well
*  as its corresponding globl index.
*
*  Arguments
*  =========
*
*  V1        (local input/local output) DOUBLE PRECISION array of
*            dimension 2.  The first maximum absolute value element and
*            its global index. V1(1) = AMAX, V1(2) = INDX.
*
*  V2        (local input) DOUBLE PRECISION array of dimension 2.
*            The second maximum absolute value element and its global
*            index. V2(1) = AMAX, V2(2) = INDX.
*
*  =====================================================================
*
*     .. Intrinsic Functions ..
      INTRINSIC          abs
*     ..
*     .. Executable Statements ..
*
      IF( abs( v1( 1 ) ).LT.abs( v2( 1 ) ) ) THEN
         v1( 1 ) = v2( 1 )
         v1( 2 ) = v2( 2 )
      END IF
*
      RETURN
*
*     End of DCOMBAMAX
*
      END
*
      SUBROUTINE dcombssq( V1, V2 )
*
*  -- ScaLAPACK tools routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Array Arguments ..
      DOUBLE PRECISION   V1( 2 ), V2( 2 )
*     ..
*
*  Purpose
*  =======
*
*  DCOMBSSQ does a scaled sum of squares on two scalars.
*
*  Arguments
*  =========
*
*  V1        (local input/local output) DOUBLE PRECISION array of
*            dimension 2.  The first scaled sum. V1(1) = SCALE,
*            V1(2) = SUMSQ.
*
*  V2        (local input) DOUBLE PRECISION array of dimension 2.
*            The second scaled sum. V2(1) = SCALE, V2(2) = SUMSQ.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO
      parameter( zero = 0.0d+0 )
*     ..
*     .. Executable Statements ..
*
      IF( v1( 1 ).GE.v2( 1 ) ) THEN
         IF( v1( 1 ).NE.zero )
     $      v1( 2 ) = v1( 2 ) + ( v2( 1 ) / v1( 1 ) )**2 * v2( 2 )
      ELSE
         v1( 2 ) = v2( 2 ) + ( v1( 1 ) / v2( 1 ) )**2 * v1( 2 )
         v1( 1 ) = v2( 1 )
      END IF
*
      RETURN
*
*     End of DCOMBSSQ
*
      END
*
      SUBROUTINE dcombnrm2( X, Y )
*
*  -- ScaLAPACK tools routine (version 1.7) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     May 1, 1997
*
*     .. Scalar Arguments ..
      DOUBLE PRECISION   X, Y
*     ..
*
*  Purpose
*  =======
*
*  DCOMBNRM2 combines local norm 2 results, taking care not to cause
*  unnecessary overflow.
*
*  Arguments
*  =========
*
*  X       (local input) DOUBLE PRECISION
*  Y       (local input) DOUBLE PRECISION
*          X and Y specify the values x and y. X and Y are supposed to
*          be >= 0.
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      parameter( one = 1.0d+0, zero = 0.0d+0 )
*     ..
*     .. Local Scalars ..
      DOUBLE PRECISION   W, Z
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          max, min, sqrt
*     ..
*     .. Executable Statements ..
*
      w = max( x, y )
      z = min( x, y )
*
      IF( z.EQ.zero ) THEN
         x = w
      ELSE
         x = w*sqrt( one+( z / w )**2 )
      END IF
*
      RETURN
*
*     End of DCOMBNRM2
*
      END