* ************************************************************************* 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