void dset(int n, double da, double *dx, int incx)
{
  int i;

  for (i=0; i<=n-1; i=i+incx) 
     dx[i]=da;
}