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

      subroutine psprdt( fcalc, lfc, a, hptr, v, accum,
     +                   elvar, elptr,  elst,  w1, w2 )

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

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

*     This routine computes the product of a symmetric partially 
*     separable matrix A times the vector V and adds it to the 
*     vector ACCUM.

*     Only the element functions of indices FCALC(i) for i=1 up to
*     i=LFC will be used.

*     Each of the element matrices is stored column-wise as lower
*     triangle, the upper triangle being defined implicitely by symmetry.

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

*     fcalc   ( int )

*       input  : array containing the indices of the element
*                functions that will be used in the product.
*       output : unmodified.

*     lfc     ( int )

*       input  : the number of element functions used.
*       output : unmodified.

*     a      ( dble )

*       input  : the matrix used for the matrix vector product.
*                Each of the element matrices are stored  
*                column-wise in A, as lower triangle.
*       output : unmodified.

*     hptr   ( int )

*       input  : array whose ith value is the position of the
*                first component of the ith element matrix
*                in A.
*       output : unmodified.

*     v      ( dble )

*       input  : the vector used for the matrix vector product.
*       output : unmodified.

*     accum   ( dble )

*       input  : meaningless.
*       output : vector in which is stored the result of the
*                product of the matrix A times the vector V.

*     elvar   ( int )

*       input  : array containing the indices of the variables 
*                in the first element, followed by those in the
*                second element, etc. 
*       output : unmodified.

*     elptr   ( int )

*       input  : array whose kth value is the position of the
*                first variable of element k, in the list ELVAR.
*       output : unmodified.

*     elst    ( int )

*       input  : array containing the status of the element
*                functions.
*       output : unmodified.

*     w1      ( dble )

*       input  : array used as workspace.
*       output : meaningless.

*     w2      ( dble )

*       input  : array used as workspace.
*       output : meaningless.

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

*     ddot, daxpy (blasd)
*     gathr0, scattr, range, elint.

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

*     D. Tuyttens 

*     Remark :
*     --------

*     The coherence of the information contained in the input is not checked
*     by the routine.

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

*     Routine parameters

      integer           elvar(*), elst(*), elptr(*), hptr(*),
     +                  fcalc(*), lfc
      double precision  a(*), v(*),  accum(*),  w1(*), w2(*)

*     Internal variables

      integer           iel, elbeg, eldim, intdim, aptr, beg,
     +                  na, i, l, n1, ifc
      logical           elint
      double precision  ddot

*     Common specifications

      integer           arcs, nodes, elem
      common / prbdim / arcs, nodes, elem

      double precision  zero, one, two, three, half, tenm1, tenm2, tenm4
      common / prbcst / zero, one, two, three, half, tenm1, tenm2, tenm4
*
*     Loop on the considered element functions.
*
      do 50 ifc = 1 , lfc
*
*       Get informations about element function.
*
        iel = fcalc(ifc)
        if( iel .gt. 0 ) then
          aptr  = hptr(iel)
          elbeg = elptr(iel)
          eldim = elptr(iel+1) - elbeg
          beg   = elvar(elbeg)
*
          if( eldim.eq.1 ) then
*
*           Treating the case of a one dimensional
*           element function.
*
            beg = elvar(elbeg)
            accum(beg) = accum(beg) + a(aptr)*v(beg)
*
          else
*
            call gathr0(eldim, w1, 1, v, elvar(elbeg) )
            if( elint(elst(iel)) ) then
*
*             Treating the internal dimension case.
*
              call range(iel, eldim, intdim, w1, w2, 0 )
              call range(iel, eldim, intdim, w1, w2, 1 )
              n1 = intdim + 1
              na = aptr + 1
              do 10 i = 1 , intdim-1
                l = intdim - i
                w1(i) = ddot(l, a(na), 1, w2(i+1), 1 )
                na = na + l + 1
 10           continue
              w1(intdim) = zero
              na = aptr
              do 20 i = 1 , intdim
                l = n1 - i
                call daxpy(l, w2(i), a(na), 1, w1(i), 1 )
                na = na + l
 20           continue
              call range(iel, eldim, intdim, w1, w2, 2 )
*
            else
*
*             Treating the elemental dimension case.
*
              na = aptr + 1
              do 30 i = 1 , eldim-1
                l = eldim - i
                w2(i) = ddot(l, a(na), 1, w1(i+1), 1 )
                na = na + l + 1
 30           continue
              w2(eldim) = zero
              na = aptr
              n1 = eldim + 1
              do 40 i = 1 , eldim
                l = n1 - i
                call daxpy(l, w1(i), a(na), 1, w2(i), 1 )
                na = na + l
 40           continue
*
            endif
*
*           Adds the result W2 of the product of element matrix 
*           IEL times V to vector ACCUM.
*
            call scattr(eldim, w2, 1, accum, elvar(elbeg) )
* 
          endif
        endif
 50   continue
*
      return
      end