source: PSPA/parmelaPSPA/trunk/spl1d1.f @ 315

      subroutine spl1d1(n,x,f,w,iop,ij,a,b,c)
c     where n= number of points in the interpolation
c           x= origin of table of independent variable
c           f= origin of table of dependent variable
c           w= an array of dimension n which contains the calculated
c              second derivatives upon return
c         iop= an array of dimension 2 which contains combinations of
c              the integers 1 thru 5 used to specify the boundary
c              conditions
c           ij= spacing in the f and w tables
c        a,b,c= arrays of dimension n used for temporary storage
c
c--------------------------------------------------------------------------
c
      include 'ucom.h'
c
      dimension iop(2),x(n),f(n*ij),w(n*ij),a(n),b(n),c(n)
c--------------------------------------------------------------------------
c*
      k=n-1
      a(2)=-(x(2)-x(1))/6.
      b(2)=(x(3)-x(1))/3.
      w(ij+1)=(f(2*ij+1)-f(ij+1))/(x(3)-x(2))-(f(ij+1)-f(1))
     1/(x(2)-x(1))
      if (n-3)3,4,3
    3 do 10 i=3,k
      m=(i-1)*ij+1
      j1=m+ij
      j2=m-ij
      con=(x(i+1)-x(i-1))/3.
      don=(x(i)-x(i-1))/6.
      b(i)=con-(don**2)/b(i-1)
      e=(f(j1)-f(m))/(x(i+1)-x(i))-(f(m)-f(j2))/(x(i)-x(i-1))
      w(m)=e-(don*w(j2))/b(i-1)
   10 a(i)=-(don*a(i-1))/b(i-1)
    4 k1=(n-2)*ij+1
      c(n-1)=-((x(n)-x(n-1))/6.)/b(n-1)
      w(k1)=w(k1)/b(n-1)
      a(n-1)=a(n-1)/b(n-1)
      k2=k-1
      if (n-3)7,8,7
    7 do 20 i=2,k2
      j=n-i
      con=(x(j+1)-x(j))/6.
      a(j)=(a(j)-con*a(j+1))/b(j)
      c(j)=-(con*c(j+1))/b(j)
      k3=(j-1)*ij+1
      m=k3+ij
   20 w(k3)=(w(k3)-con*w(m))/b(j)
    8 k4=(n-1)*ij+1
      if (iop(1)-5) 201,200,201
  201 c1=w(1)
      if (iop(2)-5) 203,202,203
  203 c2=w(k4)
      go to 205
  200 if (n-4)300,302,302
  302 a1=x(1)-x(2)
      a2=x(1)-x(3)
      a3=x(1)-x(4)
      a4=x(2)-x(3)
      a5=x(2)-x(4)
      a6=x(3)-x(4)
c      w(1)=f(1)*(1./a1+1./a2+1./a3)-a2*a3*f(ij+1)/(a1*a4*a5)+
c     1 a1*a3*f(2*ij+1)/(a2*a4*a6 )-a1*a2*f(3*ij+1)/(a3*a5*a6)
      w(1)=f(1)*(1./a1+1./a2+1./a3)-a2*a3*f(ij+1)/(a1*a4*a5)
      w(1)=w(1)+a1*a3*f(2*ij+1)/(a2*a4*a6 )-a1*a2*f(3*ij+1)/(a3*a5*a6)
      go to 201
  202 if (n-4)300,303,303
  303 b1=x(n)-x(n-3)
      b2=x(n)-x(n-2)
      b3=x(n)-x(n-1)
      b4=x(n-1)-x(n-3)
      b5=x(n-1)-x(n-2)
      b6=x(n-2)-x(n-3)
      l1=k4-ij
      l2=l1-ij
      l3=l2-ij
c      w(k4)=-b2*b3*f(l3)/(b6*b4*b1)+b1*b3*f(l2)/(b6*b5*b2)
c     1 -b1*b2*f(l1)/(b4*b5*b3)+f(k4)*(1./b1+1./b2+1./b3)
      w(k4)=-b2*b3*f(l3)/(b6*b4*b1)+b1*b3*f(l2)/(b6*b5*b2)
      w(k4)=w(k4)-b1*b2*f(l1)/(b4*b5*b3)+f(k4)*(1./b1+1./b2+1./b3)
      go to 203
  205 do 50 i=1,k
      m=(i-1)*ij+1
   60 mk=iop(1)
      go to (62,64,66,68,66),mk
   62 if (i-1)71,63,71
   63 a(1)=-1.
      c(1)=0.
      go to500
   71 bob=0.
      go to 500
   64 if (i-1)73,76,73
   76 a(1)=-1.
      c(1)=0.
      w(1)=0.
      go to 500
   73 if (i-2)81,81,82
   81 bob=-c1
      go to 500
   82 bob=0.
      go to 500
   66 if (i-1)83,84,83
   84 a(1)=-(x(2)-x(1))/3.
      c(1)=0.
      w(1)=-c1+(f(ij+1)-f(1))/(x(2)-x(1))
      go to 500
   83 if (i-2)85,85,86
   85 bob=(x(2)-x(1))/6.
      go to 500
   86 bob=0.
      go to 500
   68 if (i-1)87,88,87
   88 a(1)=-1.
      c(1)=1.
      w(1)=0.
      go to 500
   87 bob=0.
  500 ml=iop(2)
      go to (120,130,140,150,140),ml
  120 if (i-1)121,122,121
  122 a(n)=0.
      c(n)=-1.
      go to 70
  121 bill=0.
      go to 70
  130 if (i-1)131,132,131
  132 a(n)=0.
      c(n)=-1.
      w(k4)=0.
      go to 70
  131 if (i-k)134,133,134
  133 bill=-c2
      go to 70
  134 bill=0.
      go to 70
  140 if (i-1)141,142,141
  142 a(n)=0.
      c(n)=(x(n-1)-x(n))/3.
      w(k4)=c2-(f(k4)-f(k1))/(x(n)-x(n-1))
      go to 70
  141 if (i-k)143,144,143
  144 bill=(x(n)-x(n-1))/6.
      go to 70
  143 bill=0.
      go to 70
  150 if (i-1)151,152,151
  152 a(n)=0.
      c(n)=(x(n-1)+x(1)-x(n)-x(2))/3.
      w(k4)=(f(ij+1)-f(1))/(x(2)-x(1))-(f(k4)-f(k1))/(x(n)-x(n-1))
      go to 70
  151 if (i-2)153,154,153
  154 bill=(x(2)-x(1))/6.
      go to 70
  153 if (i-k)155,156,155
  156 bill=(x(n)-x(n-1))/6.
      go to 70
  155 bill=0.
   70 if (i-1)80,50,80
   80 w(1)=w(1)-bob*w(m)
      w(k4)=w(k4)-bill*w(m)
      a(1)=a(1)-bob*a(i)
      a(n)=a(n)-bill*a(i)
      c(1)=c(1)-bob*c(i)
      c(n)=c(n)-bill*c(i)
   50 continue
      go to 100
  100 con=a(1)*c(n)-c(1)*a(n)
      d1=-w(1)
      d2=-w(k4)
      w(1)=(d1*c(n)-c(1)*d2)/con
      w(k4)=(a(1)*d2-d1*a(n))/con
      do 110 i=2,k
      m=(i-1)*ij+1
  110 w(m)=w(m)+a(i)*w(1)+c(i)*w(k4)
      go to 305
  300 write(ndiag,*)' hsp1d1,  n.lt.4 results incorrect '
      call appendparm
      stop ' Abnormal stop spl1d1 '
  305 return
      end
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++