*     *************************************************************************

      double precision function ddotcd( n, ind, x, y )

*     *************************************************************************

*     Purpose :
*     ---------

*     This routine computes the dot product of vector X
*     with vector Y.
*     The components of these vectors that are considered
*     are the components of indices IND(i) for i=1 up to
*     i=n.

*     Parameters :
*     ------------

*     n       ( int )

*       input  : number of components of X and Y that will
*                be considered.
*                It is also the length of array IND.
*       output : unmodified.

*     ind     ( int )

*       input  : array of length N, containing the indices
*                of the components of X and Y that will be
*                considered.
*       output : unmodified.

*     x       ( dble )

*       input  : a vector containing the componenets of
*                X that will be considered in the dot product.
*       output : unmodified.

*     y       ( dble )

*       input  : a vector containing the componenets of
*                Y that will be considered in the dot product.
*       output : unmodified.

*     ddotcd  ( dble )

*       input  : meaningless.
*       output : the dot product of vector X with vector Y,
*                only the components whose indices are IND(i)
*                for i=1 up to i=n are considered.

*     Routines used :
*     ---------------

*     None.

*     Programming :
*     -------------

*     D. Tuyttens

*     ========================================================================

*     Function parameters

      integer           n, ind(*)
      double precision  x(*), y(*)

*     Internal variables

      integer           i, j

      double precision  zero, one, two, three, half, tenm1, tenm2, tenm4
      common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4

      ddotcd = zero
      do 10 i = 1 , n
        j = ind(i)
        ddotcd = ddotcd + ( x(j)*y(j) )
 10   continue
*
      return
      end