spider_spline.f

Go to the documentation of this file.

      module spider_spline

      contains

      SUBROUTINE intrptau(PXIN,PYIN,PYINNEW,PY2,KNIN,PXOUT,PYOUT,PYOUTP,
     +  pyoutpp,knout,kopt,ptaus,pwork,pamat,mdamat,pbclft,pbcrgt)
C     =================================================================
C
C     NOTE: THIS ROUTINE INCLUDES THE STANDARD CUBIC SPLINE IF PTAUS=0:
C           THEN PYINNEW IS NOT USED AND MDAMAT=3 IS SUFFICIENT
C           (=> PYINNEW(1) OR PYINNEW=PYIN IS OK)
C
C     Interpolate (pxin,pyin) on (pxout,pyout) using
C     Hirshman fitted cubic spline with ptaus value or
C     standard cubic spline if PTAUS=0
C
C     KOPT = 0: ONLY INTERPOLATE FUNCTION INTO PYOUT
C     KOPT = 1: INTERPOLATE FUNCTION INTO PYOUT AND 1ST DER. INTO PYOUTP
C     KOPT = 2: AS KOPT=1 PLUS 2ND DER. INTO PYOUTPP
C
C     BOUNDARY CONDITIONS FOR CUBIC SPLINE: PBCLFT AND PBCRGT
C     THEN IF
C     = 2  : B.C.: 2ND DERIVATIVE (NATURAL)
C     = 1  : B.C.: 1ST DERIVATIVE FROM CUBIC LAGRANGE INTERPOLATION
C     = 0  : B.C.: 1ST DERIVATIVE = 0.0
C     = 3  : B.C.: 1ST DERIVATIVE FROM LINEAR INTERPOLATION
C     OTHERWISE: B.C: 1ST DERIVATIVE = PBCLFT OR PBCRGT
C

      IMPLICIT REAL*8(a-h,o-z)

      dimension pxin(knin), pyin(knin), pxout(knout), pyout(knout),
     +  pyoutp(knout), pyoutpp(knout), py2(knin), pwork(knin),
     +  pamat(mdamat,knin), pyinnew(knin)
C-----------------------------------------------------------------------
C
      zeps = 1.0d-12
C
CL    1. PREPARE SPLINE (DEFAULT <=> PKBCOPT.EQ.0)
C
C     LEFT
C
      zyp1 = pbclft
      IF (abs(pbclft - 1.d0) .LE. zeps) zyp1 = -1.d+32
      IF (abs(pbclft - 2.d0) .LE. zeps) zyp1 =  1.d+32
      IF (abs(pbclft - 3.d0) .LE. zeps) zyp1 = -1.d+34
C
C     RIGHT
C
      zypn = pbcrgt
      IF (abs(pbcrgt - 1.d0) .LE. zeps) zypn = -1.d+32
      IF (abs(pbcrgt - 2.d0) .LE. zeps) zypn=  1.d+32
      IF (abs(pbcrgt - 3.d0) .LE. zeps) zypn = -1.d+34
C
      CALL cubsplfit(pxin,pyin,pyinnew,knin,zyp1,zypn,py2,ptaus
     +  ,pwork,pamat,mdamat)
C
CL    2. COMPUTE INTERPOLATED VALUE AT EACH PXOUT
C
      DO i=1,knout
        IF (ptaus .EQ. 0.d0) THEN
c%OS        CALL SPLINT(PXIN,PYIN,PY2,KNIN,PXOUT(I),ZY,ZYP,ZYPP)
          CALL splintend(pxin,pyin,py2,knin,pxout(i),zy,zyp,zypp,1)
        ELSE
c%OS        CALL SPLINT(PXIN,PYINNEW,PY2,KNIN,PXOUT(I),ZY,ZYP,ZYPP)
          CALL splintend(pxin,pyinnew,py2,knin,pxout(i),zy,zyp,zypp,1)
        ENDIF
        pyout(i) = zy
        IF (kopt .GE. 1) pyoutp(i) = zyp
        IF (kopt .EQ. 2) pyoutpp(i) = zypp
      END DO
C
      RETURN
      END SUBROUTINE
C-----------------------------------------------------------------------
      SUBROUTINE cubsplfit(X,Y,YNEW,N,YP1,YPN,Y2,TAUS,WORK,AMAT,
     +  mdamat)
C
C     PREPARE SECOND DERIVATIVE OF CUBIC SPLINE INTERPOLATION AND NEW
C     VALUES OF Y AT NODES YNEW FITTED SUCH THAT CHI**2 + TAUS*F''**2
C     IS MINIMIZED ACCORDING TO HIRSHMAN ET AL, PHYS. PLASMAS 1 (1994) 2280.
C     TAUS = TAU*SIGMA_K OF PAPER ASSUMING SIGMA_K CONSTANT.
C
C     SETTING TAUS=0., ONE FINDS THE USUAL CUBIC SPLINE INT. WITH CHI**2=0
C     TAUS LARGE => FIT CLOSER TO STRAIGHT LINE (SECOND DERIV.=0)
C
C     IF LAPACK ROUTINES NOT AVAILABLE, USE NONSYM.F AND REMOVE "c%nonsym"
C
C     IF TAUS=0, YNEW NOT USED => YNEW(1) OR YNEW=Y IS OK
C

      IMPLICIT REAL*8(a-h,o-z)

      parameter(mnonsym=501)
      dimension x(n), y(n), y2(n), work(n),
     +  amat(mdamat,n), ynew(n)
     +  ,anonsym(mnonsym*7), ipivot(mnonsym)
c%nonsym
c%OS         INCLUDE 'CUCCCC.inc'
c.......................................................................
c.......................................................................
C*COMDECK CUCCCC
C ----------------------------------------------------------------------
C --     STATEMENT FUNCTION FOR CUBIC INTERPOLATION                   --
C --                         23.04.88            AR        CRPP       --
C --                                                                  --
C -- CUBIC INTERPOLATION OF A FUNCTION F(X)                           --
C -- THE EIGHT ARGUMENTS A1,A2,A3,A4,B1,B2,B3,B4 ARE DEFINED BY:      --
C -- F(B1) = A1 , F(B2) = A2 , F(B3) = A3 , F(B4) = A4                --
C ----------------------------------------------------------------------
C
         fa3(a1,a2,a3,a4,b1,b2,b3,b4) =
     f(a1-a2) / ((b1-b2)*(b2-b4)*(b2-b3)) +
     f(a1-a3) / ((b4-b3)*(b3-b1)*(b3-b2)) +
     f(a1-a4) / ((b1-b4)*(b2-b4)*(b3-b4))
         fa2(a1,a2,a3,a4,b1,b2,b3,b4) =
     f(a1-a2) / ((b2-b1)*(b3-b2)) +
     f(a3-a1) / ((b3-b1)*(b3-b2)) -
     f(b1+b2+b3) * fa3(a1,a2,a3,a4,b1,b2,b3,b4)
         fa1(a1,a2,a3,a4,b1,b2,b3,b4) =
     f(a1-a2) / (b1-b2) -
     f(b1+b2) * fa2(a1,a2,a3,a4,b1,b2,b3,b4) -
     f(b1*b1+b1*b2+b2*b2) * fa3(a1,a2,a3,a4,b1,b2,b3,b4)
         fa0(a1,a2,a3,a4,b1,b2,b3,b4) =
     f        a1 -
     f        b1 * (fa1(a1,a2,a3,a4,b1,b2,b3,b4) +
     f              b1 * (fa2(a1,a2,a3,a4,b1,b2,b3,b4) +
     f                    b1 * fa3(a1,a2,a3,a4,b1,b2,b3,b4)))
C ----------------------------------------------------------------------
C -- FCCCC0 GIVES THE VALUE OF THE FUNCTION AT POINT PX:              --
C -- FCCCC0(......,PX) = F(PX)                                        --
C ----------------------------------------------------------------------
        fcccc0(a1,a2,a3,a4,b1,b2,b3,b4,px) =
     f              fa0(a1,a2,a3,a4,b1,b2,b3,b4) +
     f              px * (fa1(a1,a2,a3,a4,b1,b2,b3,b4) +
     f                    px * (fa2(a1,a2,a3,a4,b1,b2,b3,b4) +
     f                          px * fa3(a1,a2,a3,a4,b1,b2,b3,b4)))
C ----------------------------------------------------------------------
C -- FCCCC1 GIVES THE VALUE OF THE FIRST DERIVATIVE OF F(X) AT PX:    --
C -- FCCCC1(......,PX) = DF/DX (PX)                                   --
C ----------------------------------------------------------------------
        fcccc1(a1,a2,a3,a4,b1,b2,b3,b4,px) =
     f              fa1(a1,a2,a3,a4,b1,b2,b3,b4) +
     f              px * (2.d0 * fa2(a1,a2,a3,a4,b1,b2,b3,b4) +
     f                    3.d0 * px * fa3(a1,a2,a3,a4,b1,b2,b3,b4))
C ----------------------------------------------------------------------
C -- FCCCC2 GIVES THE VALUE OF THE SECOND DERIVATIVE OF F(X) AT PX:   --
C -- FCCCC2(......,PX) = D2F/DX2 (PX)                                 --
C ----------------------------------------------------------------------
         fcccc2(a1,a2,a3,a4,b1,b2,b3,b4,px) =
     f             2.d0 * fa2(a1,a2,a3,a4,b1,b2,b3,b4) +
     f             6.d0 * fa3(a1,a2,a3,a4,b1,b2,b3,b4) * px
C ----------------------------------------------------------------------
C -- FCCCC3 GIVES THE VALUE OF THE THIRD DERIVATIVE OF F(X) AT PX:     -
C -- FCCCC3(......,PX) = D3F/DX3 (PX)                                  -
C ----------------------------------------------------------------------
         fcccc3(a1,a2,a3,a4,b1,b2,b3,b4,px) =
     f                      6.d0* fa3(a1,a2,a3,a4,b1,b2,b3,b4)
c.......................................................................
c.......................................................................
C
C-----------------------------------------------------------------------
C
      DO i=1,n
        y2(i) = 0.d0
      END DO
      DO i=1,mdamat
        DO j=1,n
          amat(i,j) = 0.d0
        END DO
      END DO
C
C     PREPARE 1 / H_K
C
      DO k=1,n-1
        work(k) = 1.d0/ (x(k+1) - x(k))
      END DO
C
      ztaueff = taus
C
C     PREPARE BAND WIDTH
C
      iup = 2
      IF (ztaueff .EQ. 0.d0) iup = 1
      iband = 2*iup + 1
c%OS      IDIAG = IUP + 1
c%OS      IF (MDAMAT .LT. IUP+1) THEN
c%OS        PRINT *,' MDAMAT= ',MDAMAT,' < IUP+1= ',IUP
c%OS        STOP
c%OS      ENDIF
c%nonsym
      idiag = iband
      IF (mnonsym*7 .LT. (3*iup+1)*n) THEN
        print *,' MMNONSYM*7= ',mnonsym*7,' < (3*IUP+1)*N= ',(3*iup+1)*n
        stop
      ENDIF
c%nonsym
C
C     CONSTRUCT MATRIX AND R.H.S
C
c.......................................................................
C     AS MATRIX SYMMETRIC, COMPUTE ONLY UPPER PART
C
C     K=2,N-2 (BULK PART)
C
      DO k=2,n-2
C     A(K,K)
        amat(idiag,k) = (1.d0/work(k)+1.d0/work(k-1)) / 3.d0
     +    + 2.d0*ztaueff*(work(k)**2 + work(k)*work(k-1)+work(k-1)**2)
C     A(K,K+1)
        amat(idiag-1,k+1) = 1.d0/work(k)/6.d0
     +    - ztaueff*work(k)*(work(k+1)+2.d0*work(k)+work(k-1))
C     A(K,K+2)
        IF (iup .EQ. 2) amat(idiag-2,k+2) = ztaueff*work(k+1)*work(k)
C     B(K)
        y2(k) = (y(k+1)-y(k))*work(k) - (y(k)-y(k-1))*work(k-1)
C
      END DO
C
C     FOR K=1, N-1, N: SAME AS ABOVE WITH WORK1/2(-1;0;N;N+1) = 0
C
C     K=1
C
      k = 1
      amat(idiag,k) = 1.d0/work(k) / 3.d0
     +  + 2.d0*ztaueff*work(k)**2
      amat(idiag-1,k+1) = 1.d0/work(k)/6.d0
     +  - ztaueff*work(k)*(work(k+1)+2.d0*work(k))
      IF (iup .EQ. 2) amat(idiag-2,k+2) = ztaueff*work(k+1)*work(k)
      y2(k) = (y(k+1)-y(k))*work(k)
C
C     K=N-1
C
      k = n-1
      amat(idiag,k) = (1.d0/work(k)+1.d0/work(k-1)) / 3.d0
     +  + 2.d0*ztaueff*(work(k)**2 + work(k)*work(k-1)+work(k-1)**2)
      amat(idiag-1,k+1) = 1.d0/work(k)/6.d0
     +  - ztaueff*work(k)*(2.d0*work(k)+work(k-1))
      y2(k) = (y(k+1)-y(k))*work(k) - (y(k)-y(k-1))*work(k-1)
C
C     K = N
C
      k = n
      amat(idiag,k) = 1.d0/work(k-1) / 3.d0
     +    + 2.d0*ztaueff*work(k-1)**2
      y2(k) = - (y(k)-y(k-1))*work(k-1)
C
C.......................................................................
C     BOUNDARY CONDITIONS
C
      IF (yp1 .GT. .99d30) THEN
C     SECOND DERIVATIVE = 0 (NATURAL B.C.) AND SYMMETRIZE
        y2(1) = 0.d0
        amat(idiag,1) = 1.d0
        amat(idiag-1,2) = 0.d0
        IF (iup .EQ. 2) amat(idiag-2,3) = 0.d0
        zyp1 = 0.d0
      ELSE IF (-yp1 .GT. .99d+34) THEN
        zyp1 = (y(2)-y(1))/(x(2)-x(1))
      ELSE IF (-yp1 .GT. .99d+30) THEN
        zyp1 = fcccc1(y(1),y(2),y(3),y(4),x(1),x(2),x(3),x(4),x(1))
      ELSE
        zyp1 = yp1
      ENDIF
      IF (ypn .GT. .99d30) THEN
C     SECOND DERIVATIVE = 0 (NATURAL B.C.) AND SYMMETRIZE
        y2(n) = 0.d0
        amat(idiag,n) = 1.d0
        amat(idiag-1,n) = 0.d0
        IF (iup .EQ. 2) amat(idiag-2,n) = 0.d0
        zypn = 0.d0
      ELSE IF (-ypn .GT. .99d+34) THEN
        zypn = (y(n)-y(n-1))/(x(n)-x(n-1))
      ELSE IF (-ypn .GT. .99d+30) THEN
        zypn = fcccc1(y(n-3),y(n-2),y(n-1),y(n),
     +    x(n-3),x(n-2),x(n-1),x(n),x(n))
      ELSE
        zypn = ypn
      ENDIF
C
      y2(1) = y2(1) - zyp1
      y2(n) = y2(n) + zypn
C
C     SOLVE SYSTEM
C
      idima = mdamat
      idimrhs = n
      irhs = 1
c%OSc   CRAY
c%OSc%OS      CALL SPBTRF('U',N,IUP,AMAT,IDIMA,INFO)
c%OSC   SUN, SGI
c%OS      CALL DPBTRF('U',N,IUP,AMAT,IDIMA,INFO)
c%OS      IF (INFO .EQ. 0) THEN
c%OSc   CRAY
c%OSc%OS        CALL SPBTRS('U',N,IUP,IRHS,AMAT,IDIMA,Y2,IDIMRHS,INFO2)
c%OSC   SUN, SGI
c%OS        CALL DPBTRS('U',N,IUP,IRHS,AMAT,IDIMA,Y2,IDIMRHS,INFO2)
c%OS      ELSE
c%OS        PRINT *,' ERROR IN SPBTRF: INFO = ',INFO
c%OS        STOP 'INFO'
c%OS      ENDIF
C
c%nonsym
      DO i=1,iband*n
        anonsym(i) = 0.0
      ENDDO
      DO i=1,n
C     UPPER PART
        DO j=max(1,i),min(i+iup,n)
          anonsym((i-1)*iband+j-i+iup+1) = amat(idiag+i-j,j)
        ENDDO
C     LOWER PART
        DO j=max(1,i-iup),min(i-1,n-1)
          anonsym((i-1)*iband+j-i+iup+1) = amat(idiag+j-i,i)
        ENDDO
      ENDDO
      CALL nonsym(anonsym,amat,y2,n,iup,iup,1.0d-06,info2)
c%nonsym
C
C     ADAPT Y AT NODES
C
      IF (ztaueff .NE. 0.d0) THEN
C
        DO k=2,n-1
          ynew(k) = y(k) - ztaueff* ((y2(k+1)-y2(k))*work(k)
     +      - (y2(k)-y2(k-1))*work(k-1))
        ENDDO
        ynew(1) = y(1) - ztaueff * (y2(2)-y2(1))*work(1)
        ynew(n) = y(n) + ztaueff * (y2(n)-y2(n-1))*work(n-1)
C
      ENDIF

      IF (info2 .LT. 0) THEN
        print *,' ERROR IN SPBTRS: INFO2 = ',info2
        stop 'INFO2'
      ENDIF
C
      RETURN
      END SUBROUTINE
C-----------------------------------------------------------------------
      SUBROUTINE spline(X,Y,N,YP1,YPN,Y2)

      IMPLICIT REAL*8(a-h,o-z)

      parameter(nmax=100)
      dimension x(n),y(n),y2(n),u(nmax)
C
      IF (n .GT.nmax) THEN
        print *,' NMAX TOO SMALL IN SPLINE: N,NMAX= ',n,nmax
        stop
      ENDIF
      IF (yp1.GT..99d30) THEN
        y2(1)=0.d0
        u(1)=0.d0
      ELSE
C%OS
        IF (-yp1 .GT. .99d+30) THEN
          yp1 = (y(2)-y(1))/(x(2)-x(1))
        ENDIF
C%OS
        y2(1)=-0.5d0
        u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      ENDIF
      DO 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.d0
        y2(i)=(sig-1.d0)/p
        u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
     *      /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
11    CONTINUE
      IF (ypn.GT..99d30) THEN
        qn=0.d0
        un=0.d0
      ELSE
C%OS
        IF (-ypn .GT. .99d+30) THEN
          ypn = (y(n)-y(n-1))/(x(n)-x(n-1))
        ENDIF
C%OS
        qn=0.5d0
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      ENDIF
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      DO 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    CONTINUE
      RETURN
      END SUBROUTINE
C-----------------------------------------------------------------------
      SUBROUTINE splint(XA,YA,Y2A,N,X,Y,YP,YPP)

      IMPLICIT REAL*8(a-h,o-z)

      dimension xa(n),ya(n),y2a(n)
      klo=1
      khi=n
1     IF (khi-klo.GT.1) THEN
        k=(khi+klo)/2
        IF(xa(k).GT.x)THEN
          khi=k
        ELSE
          klo=k
        ENDIF
      goto 1
      ENDIF
      h=xa(khi)-xa(klo)
      IF (h.EQ.0.) stop 'BAD XA INPUT.'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+
     *      ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
      yp=(ya(khi)-ya(klo))/h -
     -  ( (3.*a*a-1.)*y2a(klo) - (3.*b*b-1.)*y2a(khi) )*h/6.
      ypp=a*y2a(klo)+b*y2a(khi)
      RETURN
      END SUBROUTINE
C-----------------------------------------------------------------------
      SUBROUTINE splintend(XA,YA,Y2A,N,X,Y,YP,YPP,KOPT)

      IMPLICIT REAL*8(a-h,o-z)

      dimension xa(n),ya(n),y2a(n)
C     KOPT = 0: NOTHING SPECIAL IF X OUT OF BOUND
C     KOPT = 1: USE QUADRATIC INTERPOLATION IF X OUT OF BOUND
C     KOPT = 2: USE CUBIC INTERPOLATION IF X OUT OF BOUND
C-----------------------------------------------------------------------
C*COMDECK CUCDCD
C ----------------------------------------------------------------------
C --     STATEMENT FUNCTION FOR CUBIC INTERPOLATION                   --
C --                         19.01.87            AR        CRPP       --
C --                                                                  --
C -- CUBIC INTERPOLATION OF A FUNCTION F(X)                           --
C -- THE SIX ARGUMENTS X1,F1,P1,X2,F2,P2 ARE DEFINED AS FOLLOWS:      --
C -- F(X1) = F1 , F(X2) = F2 , DF/DX(X1) = P1 , DF/DX(X2) = P2        --
C ----------------------------------------------------------------------
C
         fc3(x1,f1,p1,x2,f2,p2) =
     =      (2.d0* (f2 - f1) / (x1 - x2) + (p1 + p2)) / 
     /      ((x1 - x2) * (x1 - x2))
         fc2(x1,f1,p1,x2,f2,p2) =
     =      (3.d0* (x1 + x2) * (f1 - f2) / (x1 - x2) - 
     *       p1 * (x1 + 2.d0* x2) - p2 * (x2 + 2.d0* x1)) /
     /      ((x1 - x2) * (x1 - x2))
         fc1(x1,f1,p1,x2,f2,p2) =
     =      (6.d0* x1 * x2 * (f2 - f1) / (x1 - x2) + 
     *       x2 * p1 * (2 * x1 + x2) + x1 * p2 * (x1 + 2.d0* x2)) /
     /      ((x1 - x2) * (x1 - x2))
         fc0(x1,f1,p1,x2,f2,p2) =
     =      (f1 * x2**2 + f2 * x1**2 - x1 * x2 * (x2 * p1 + x1 * p2) +
     +       2.d0* x1 * x2 * (f1 * x2 - f2 * x1) / (x1 - x2)) /
     /      ((x1 - x2) * (x1 - x2))
C ----------------------------------------------------------------------
C -- FCDCD0 GIVES THE VALUE OF THE FUNCTION AT POINT PX               --
C -- FCDCD0(......,PX) = F(PX)                                        --
C ----------------------------------------------------------------------
         fcdcd0(x1,f1,p1,x2,f2,p2,px) =
     =              fc0(x1,f1,p1,x2,f2,p2) +
     +              px * (fc1(x1,f1,p1,x2,f2,p2) +
     +                    px * (fc2(x1,f1,p1,x2,f2,p2) +
     +                          px * fc3(x1,f1,p1,x2,f2,p2)))
C ----------------------------------------------------------------------
C -- FCDCD1 GIVES THE VALUE OF THE FIRST DERIVATIVE OF F(X) AT PX:    --
C -- FCDCD1(......,PX) = DF/DX (PX)                                   --
C ----------------------------------------------------------------------
         fcdcd1(x1,f1,p1,x2,f2,p2,px) =
     =              fc1(x1,f1,p1,x2,f2,p2) +
     +              px * (2.d0* fc2(x1,f1,p1,x2,f2,p2) +
     +                    3.d0* px * fc3(x1,f1,p1,x2,f2,p2))
C ----------------------------------------------------------------------
C -- FCDCD2 GIVES THE VALUE OF THE SECOND DERIVATIVE OF F(X) AT PX:   --
C -- FCDCD2(......,PX) = D2F/DX2 (PX)                                 --
C ----------------------------------------------------------------------
         fcdcd2(x1,f1,p1,x2,f2,p2,px) =
     =             2.d0* fc2(x1,f1,p1,x2,f2,p2) +
     +             6.d0* fc3(x1,f1,p1,x2,f2,p2) * px
C ----------------------------------------------------------------------
C -- FCDCD3 GIVES THE VALUE OF THE THIRD DERIVATIVE OF F(X) AT PX:    --
C -- FCDCD3(......,PX) = D3F/DX3 (PX)                                 --
C ----------------------------------------------------------------------
         fcdcd3(x1,f1,p1,x2,f2,p2,px) =
     =                      6.d0* fc3(x1,f1,p1,x2,f2,p2)
C-----------------------------------------------------------------------
C*COMDECK QUAQDQ
C ----------------------------------------------------------------------
C --     STATEMENT FUNCTION FOR QUADRATIC INTERPOLATION               --
C --                         19.01.87            AR        CRPP       --
C --                                                                  --
C -- QUADRATIC INTERPOLATION OF A FUNCTION F(X)                       --
C -- THE FIVE PARAMETERS X1,F1,P1,X2,F2    ARE DEFINED AS FOLLOWS:    --
C -- F(X1) = F1 , DF/DX(X1) = P1 , F(X2) = F2                         --
C ----------------------------------------------------------------------
C
         fd2(x1,f1,p1,x2,f2) = ((f2-f1)/(x2-x1) - p1) / (x2-x1)
         fd1(x1,f1,p1,x2,f2) = p1 - 2.*x1*fd2(x1,f1,p1,x2,f2)
         fd0(x1,f1,p1,x2,f2) = f1 - x1*(x1*fd2(x1,f1,p1,x2,f2) +
     +                                     fd1(x1,f1,p1,x2,f2))
C ----------------------------------------------------------------------
C -- FQDQ0 GIVES THE VALUE OF THE FUNCTION AT POINT PX                --
C -- FQDQ0(......,PX) = F(PX)                                         --
C ----------------------------------------------------------------------
         fqdq0(x1,f1,p1,x2,f2,px) = fd0(x1,f1,p1,x2,f2) +
     f                              px * (fd1(x1,f1,p1,x2,f2) +
     f                                    px * fd2(x1,f1,p1,x2,f2))
C ----------------------------------------------------------------------
C -- FQDQ1 GIVES THE VALUE OF THE FIRST DERIVATIVE OF F(X) AT PX      --
C -- FQDQ1(......,PX) = DF/DX (PX)                                    --
C ----------------------------------------------------------------------
         fqdq1(x1,f1,p1,x2,f2,px) = fd1(x1,f1,p1,x2,f2) +
     f                              2.* px * fd2(x1,f1,p1,x2,f2)
C ----------------------------------------------------------------------
C -- FQDQ2 GIVES THE VALUE OF THE SECOND DERIVATIVE OF F(X) AT PX     --
C -- FQDQ2(......,PX) = D2F/DX2 (PX)                                  --
C ----------------------------------------------------------------------
         fqdq2(x1,f1,p1,x2,f2) = 2.d0* fd2(x1,f1,p1,x2,f2)
C-----------------------------------------------------------------------

      klo=1
      khi=n
1     IF (khi-klo.GT.1) THEN
        k=(khi+klo)/2
        IF(xa(k).GT.x)THEN
          khi=k
        ELSE
          klo=k
        ENDIF
      goto 1
      ENDIF
      h=xa(khi)-xa(klo)
      IF (h.EQ.0.) stop 'BAD XA INPUT.'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+
     *      ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.
      yp=(ya(khi)-ya(klo))/h -
     -  ( (3.*a*a-1.)*y2a(klo) - (3.*b*b-1.)*y2a(khi) )*h/6.
      ypp=a*y2a(klo)+b*y2a(khi)

      IF (kopt .EQ. 0) RETURN
C
C     ZEPSILON: DISTANCE OUTSIDE INTERVAL RELATIVE TO H
C
      zepsilon = 1.0d-05
c
c     x outside interval: use quadratic with y and yp at edge and y-1
c
      IF (kopt .EQ. 1) THEN
C     LEFT END
        if (b.lt. -zepsilon) then
          print *,' warning points outside interval at left.',
     +      ' Use quadratic interpolation'
          a=1.d0
          b=0.d0
          ypaklo=(ya(khi)-ya(klo))/h -
     -    ( (3.d0*a*a-1.d0)*y2a(klo) - (3.d0*b*b-1.d0)*y2a(khi) )*h/6.d0
          a=0.d0
          b=1.d0
          ypakhi=(ya(khi)-ya(klo))/h -
     -    ( (3.d0*a*a-1.d0)*y2a(klo) - (3.d0*b*b-1.d0)*y2a(khi) )*h/6.d0
          y=fqdq0(xa(klo),ya(klo),ypaklo,xa(khi),ya(khi),x)
          yp=fqdq1(xa(klo),ya(klo),ypaklo,xa(khi),ya(khi),x)
          ypp=fqdq2(xa(klo),ya(klo),ypaklo,xa(khi),ya(khi))
        endif
C     RIGHT END
        if (a .lt. -zepsilon) then
          print *,' warning points outside interval at right.',
     +      ' Use quadratic interpolation'
          a=1.d0
          b=0.d0
          ypaklo=(ya(khi)-ya(klo))/h -
     -    ( (3.d0*a*a-1.d0)*y2a(klo) - (3.d0*b*b-1.d0)*y2a(khi) )*h/6.d0
          a=0.d0
          b=1.d0
          ypakhi=(ya(khi)-ya(klo))/h -
     -    ( (3.d0*a*a-1.d0)*y2a(klo) - (3.d0*b*b-1.d0)*y2a(khi) )*h/6.d0
          y=fqdq0(xa(khi),ya(khi),ypakhi,xa(klo),ya(klo),x)
          yp=fqdq1(xa(khi),ya(khi),ypakhi,xa(klo),ya(klo),x)
          ypp=fqdq2(xa(khi),ya(khi),ypakhi,xa(klo),ya(klo))
        endif

        RETURN
      ENDIF
c
c     x outside interval: use CUBIC with y and yp
c
      IF (kopt .EQ. 2) THEN
C     LEFT OR RIGHT END
        if (a.lt.-zepsilon .OR. b.LT.-zepsilon) then
          print *,' warning points outside interval.',
     +      ' Use cubic interpolation'
          a=1.d0
          b=0.d0
          ypaklo=(ya(khi)-ya(klo))/h -
     -    ( (3.d0*a*a-1.d0)*y2a(klo) - (3.d0*b*b-1.d0)*y2a(khi) )*h/6.d0
          a=0.d0
          b=1.d0
          ypakhi=(ya(khi)-ya(klo))/h -
     -    ( (3.d0*a*a-1.d0)*y2a(klo) - (3.d0*b*b-1.d0)*y2a(khi) )*h/6.d0
          y=fcdcd0(xa(klo),ya(klo),ypaklo,xa(khi),ya(khi),ypakhi,x)
          yp=fcdcd1(xa(klo),ya(klo),ypaklo,xa(khi),ya(khi),ypakhi,x)
          ypp=fcdcd2(xa(klo),ya(klo),ypaklo,xa(khi),ya(khi),ypakhi,x)
        endif
        RETURN
      ENDIF

      RETURN
      END SUBROUTINE
C
C
C
      SUBROUTINE nonsym(A,X,C,N,MLEFT,MRIGHT,EPS,NCOND)
C     *************************************************
C
C     5.3  SOLVES A REAL VALUED NONSYMMETRIC LINEAR SYSTEM  A . X  =  C
C
C     VERSION 1               APRIL 1988       KA        LAUSANNE
C     VERSION 2               MAI 1992         AJ        LAUSANNE
C
C-----------------------------------------------------------------------

      IMPLICIT REAL*8(a-h,o-z)

      dimension   a(n*(mleft+mright+1)),    x(n),    c(n)
CMPLX COMPLEX     A,       X,       C,        ZPIVOT,   ZSUM,     ZTOP
      REAL*8        a,       x,       c,        zpivot,   zsum,     ztop
      DATA        imess  / 0 /
C-----------------------------------------------------------------------
C
C     A IS A BANDMATRIX OF RANK N WITH MLEFT/MRIGHT OFF-DIAGONAL
C     ELEMENTS TO THE LEFT/RIGHT.
C     :A: AND :C: ARE DESTROYED BY :NONCYM: 
C     AT THE RETURN :C: CONTAINS THE SOLUTION VECTOR, :X: THE SQUARE
C     ROOT OF THE DIAGONAL ELEMENTS
C
C     THE ELEMENTS ARE HORIZONTALLY NUMBERED ROW AFTER ROW
C
C     FOR EASY NUMBERING, ZERO ELEMENTS HAVE BEEN INTRODUCED IN THE
C     UPPER LEFT-HAND CORNER AND IN THE LOWER RIGHT-HAND CORNER.
C     THESE ELEMENTS
C     M U S T  B E  S E T  T O  Z E R O
C     BEFORE CALLING :NONCYM: 
C
C     EXAMPLE FOR NUMBERING   MLEFT=2,  MRIGHT=3
C     -------
C
C     A(1)   A(2) I A(3)   A(4)   A(5)   A(6)
C     A(7) I A(8)   A(9)   A(10)  A(11)  A(12)
C     I                 .
C     I                        .
C
C     A PIVOT IS CONSIDERED TO BE BAD IF IT IS BY A FACTOR :EPS: SMALLER
C     THAN ITS OFF-DIAGONAL ELEMENTS. IF A BAD PIVOT IS ENCOUNTERED
C     A MESSAGE IS ISSUED (AT MOST TEN TIMES IN A RUN).
C     THE FLAG :NCOND: IS SET TO -1, OTHERWISE 0. 
C
C-----------------------------------------------------------------------
CL    0.        INITIALIZATION
C
      moffdi=mleft+mright
      mband=moffdi+1
      ratio=0.d0
      ncond=0
C     PRELIMINARY CHECK OF PIVOTS
      ipiv1=mleft+1
      ipiv2=ipiv1+(n-1)*mband
      DO 10 j=ipiv1,ipiv2,mband
CMPLX IF(CABS(A(J)) .EQ. 0.D0) GO TO 510
        IF( abs(a(j)) .EQ. 0.d0) go to 510
 10   CONTINUE
C
C-----------------------------------------------------------------------
CL    1.        PRECONDITIONNING: GET ONES ON THE DIAGONAL
C
      DO 130 jp=1,n
        jpabs=(mleft+1)+(jp-1)*mband
CMPLX X(JP)=SQRT(CABS(A(JPABS)))
        x(jp)=sqrt(abs(a(jpabs)))
C
C     DIVIDE THE EQUATIONS BY THE SQUARE ROOT OF THE PIVOT
        DO 110 j=-mleft,mright
          a(jpabs+j)=a(jpabs+j)/x(jp)
 110    CONTINUE
        c(jp)=c(jp)/x(jp)
C
C     CHANGE VARIABLES, DIVIDE COLUMN BY SQUARE ROOT OF THE PIVOT
        DO 120 ja=jpabs-mright*moffdi,jpabs+mleft*moffdi,moffdi
          IF (ja.GE.ipiv1 .AND. ja.LE.ipiv2) a(ja)=a(ja)/x(jp)
 120    CONTINUE
 130  CONTINUE
C
C-----------------------------------------------------------------------
CL    2.        CONSTRUCT UPPER TRIANGULAR MATRIX 
C
C     INITIALIZATION OF VERTICAL COUNTER
      idown=mleft
C     PIVOT INDEX FOR FIRST INCOMPLETE COLUMN
      incmpl=(n-mleft)*mband+mleft+1 
C     FIRST PIVOT AND LAST BUT ONE
      ipiv1=mleft+1
      ipiv2=ipiv1+(n-2)*mband
C
      DO 220 jpivot=ipiv1,ipiv2,mband
C     COLUMN LENGTH 
        IF(jpivot.GE.incmpl) idown=idown-1
        zpivot=a(jpivot)
C     FIRST AND LAST HORIZONTAL ELEMENT INDEX
        ihor1=jpivot+1
        ihor2=jpivot+mright
C     CHECK THE PIVOT
        zmax=0.d0
        DO 205 jhoriz=ihor1,ihor2
CMPLX ZMAX=AMAX1(ZMAX,CABS(A(JHORIZ)))
          zmax=amax1(zmax, abs(a(jhoriz)))
 205    CONTINUE
CMPLX IF(CABS(ZPIVOT).EQ.0.D0) GO TO 510
        IF( abs(zpivot).EQ.0.d0) go to 510
CMPLX RATIO=AMAX1(RATIO,ZMAX/CABS(ZPIVOT))
        ratio=amax1(ratio,zmax/ abs(zpivot))
C
        DO 210 jhoriz=ihor1,ihor2
          ztop=-a(jhoriz)/zpivot
          a(jhoriz)=ztop
C     INITIALIZATION OF ELEMENT INDEX
          ielem=jhoriz
C     INITIALIZATION OF RECTANGULAR RULE CORNER ELEMENT
          icorn=jpivot
C     LOOP DOWN THE COLUMN
          DO 211 jdown=1,idown 
            ielem=ielem+moffdi
            icorn=icorn+moffdi
            a(ielem)=a(ielem)+ztop*a(icorn)
 211      CONTINUE
 210    CONTINUE
C     TREAT THE CONSTANTS
C     INDEX OF THE FIRST ONE
        iconst=jpivot/mband+1
        ztop=-c(iconst)/zpivot
        c(iconst)=ztop
C     INITIALIZATION OF RECTANGULAR RULE CORNER ELEMENT
        icorn=jpivot
C     LOOP DOWN THE CONSTANTS 
        DO 221 jdown=1,idown 
          iconst=iconst+1
          icorn=icorn+moffdi
          c(iconst)=c(iconst)+ztop*a(icorn)
 221    CONTINUE
 220  CONTINUE
C
C-----------------------------------------------------------------------
CL    3.        BACKSUBSTITUTION
C
 300  CONTINUE
C     INITIALIZATION OF SOLUTION INDEX
      isolut=n
C     LAST PIVOT
      jpivot=ipiv2+mband
C     CHECK LAST PIVOT
CMPLX IF(CABS(A(JPIVOT)).EQ.0.D0) RATIO=1.D100
c%OS  IF( ABS(A(JPIVOT)).EQ.0.D0) RATIO=1.D100
      IF( abs(a(jpivot)).EQ.0.d0) go to 510
C     LAST UNKNOWN
      c(isolut)=c(isolut)/a(jpivot)
C     INITIALIZATION OF HORIZONTAL RANGE
      ihor1=jpivot+1
      ihor2=jpivot+mright
C
      inm1=n-1
      DO 320 j=1,inm1
        isolut=isolut-1
        ihor1=ihor1-mband
        ihor2=ihor2-mband
        zsum=-c(isolut)
        iknown=isolut
        DO 310 jhoriz=ihor1,ihor2
          iknown=iknown+1
          IF(iknown.GT.n) go to 310
          zsum=zsum+a(jhoriz)*c(iknown)
 310    CONTINUE
        c(isolut)=zsum
 320  CONTINUE
C
C-----------------------------------------------------------------------
CL    4.        PRECONDITIONNING: BACK TO ORIGINAL VARIABLES
C
      DO 410 j=1,n
        c(j)=c(j)/x(j)
 410  CONTINUE
C
C-----------------------------------------------------------------------
CL    5.        HAVE WE ENCOUNTERED A BAD PIVOT ? 
C
C%OS  IF(1./RATIO .GT. EPS) RETURN
      IF(ratio .LT. 1.d0/eps) RETURN
      ncond=-1
C     COUNT THE NUMBER OF MESSAGES ISSUED
      imess=imess+1
      IF(imess.LE.10) print 9500, ratio
      RETURN
C
C     ZERO PIVOT ENCOUNTERED
 510  CONTINUE
      print 9510 
      stop 'ZERO PIVOT IN NONCYM'
C
C-----------------------------------------------------------------------
CL    9.        FORMATS
C

 9500 FORMAT(/////20x,10(1h*),21h  message from noncym,2x,10(1h*)///
     +  20x,'BAD PIVOT ENCOUNTERED, RATIO: ',1pe12.3///
     +  20x,43(1h*)///)
 9510 FORMAT(/////20x,10(1h*),'   ZERO PIVOT IN NONCYM'///)
C
      END SUBROUTINE

      end module spider_spline