!     Last change:  X     6 Jan 2005   10:35 am



MODULE pade2
CONTAINS

FUNCTION gauss(a_in,b_in) RESULT(x)

   IMPLICIT NONE

   REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: a_in
   REAL(KIND=8), DIMENSION(:), INTENT(IN) :: b_in

   REAL(KIND=8):: add, dog, cat, choice,z
   REAL(KIND=8), DIMENSION(1:SIZE(b_in),1:SIZE(b_in)) :: a
   REAL(KIND=8), DIMENSION(1:SIZE(b_in)) :: x,s,b
   INTEGER, DIMENSION(1:SIZE(b_in)) ::p
   INTEGER::i,j,k, sel,n

   n=SIZE(b_in)
   a=a_in
   b=b_in

   DO i=1,n
      p(i)=i
   END DO

   s=MAXVAL(ABS(a),DIM=2)


   DO k=1,n-1
      sel=k
      DO i=k,n
          dog=ABS(a(p(i),k)/s(p(i)))
          cat=ABS(a(p(sel),k)/s(p(sel)))
          IF (dog>cat) THEN
             sel=i
          ELSE
          END IF
       END DO
       choice=p(k)
       p(k)=p(sel)
       p(sel)=choice
       DO i=k+1,n
          z=a(p(i),k)/a(p(k),k)
          a(p(i),k)=z
          DO j=k+1,n
             a(p(i),j)=a(p(i),j)-z*a(p(k),j)
          END DO
       END DO
   END DO

   DO k=1,n-1
      DO i=k+1,n
         b(p(i))=b(p(i))-a(p(i),k)*b(p(k))
      END DO
   END DO

   DO i=n,1,-1
       add=b(p(i))
       DO j=i+1,n
           add=add-a(p(i),j)*x(j)
       END DO
       x(i)=add/a(p(i),i)
   END DO


   RETURN
END FUNCTION gauss

!____________________________________________________________________
FUNCTION padstep(y1,h) RESULT(x)

   IMPLICIT NONE

   REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: y1
   REAL(KIND=8), INTENT(IN) :: h

   INTEGER::   m, eqs, i,k,j,ii
   REAL(KIND=8), DIMENSION(1:SIZE(y1,1)) :: mult,add, numsum, denomsum,x
   REAL(KIND=8), DIMENSION(1:SIZE(y1,1),0:((SIZE(y1,2)-1))/2 ):: b,L,U
   REAL(KIND=8), DIMENSION(1:SIZE(y1,1),0:SIZE(y1,2)):: y
   REAL(KIND=8), DIMENSION(1:SIZE(y1,1),0:((SIZE(y1,2)-1))/2 ,0:((SIZE(y1,2)-1))/2 ) :: A



   m=SIZE(y1,2)
   m=(m-1)/2

   eqs=SIZE(y1,1)

   DO i=1,eqs
      DO j=0,2*m
        y(i,j)=y1(i,j+1)
      END DO
   END DO


 ! NOW START PADE

   !FILL MATRIX AND RHS FOR PADE
   DO ii=1,eqs
     DO i=1,m
       DO j=1,m
           A(ii,i,j)=y(ii,m+i-j)
       END DO
       b(ii,i)=-y(ii,m+i)
     END DO
   END DO

   !FULL FORWARD ELIMINATION
   DO ii=1,eqs
      DO k=1,m-1
         DO i=k+1,m
            mult(ii)=a(ii,i,k)/a(ii,k,k)
            a(ii,i,k)=0.0
            DO j=k+1,m
                a(ii,i,j)=a(ii,i,j)-mult(ii)*a(ii,k,j)
            END DO
            b(ii,i)=b(ii,i)-mult(ii)*b(ii,k)
         END DO
      END DO
   END DO

   ! BACK SUBS TO SET DENOM COEF

   DO ii=1,eqs
       L(ii,m)=b(ii,m)/a(ii,m,m)
       DO i=m-1,1,-1
           add(ii)=b(ii,i)
           DO j=m,i+1,-1
               add(ii)=add(ii)-a(ii,i,j)*L(ii,j)
           END DO
           L(ii,i)=add(ii)/a(ii,i,i)
       END DO
    END DO



    !SET NUMER COEF USING CAUCHY PRODUCT

    U=0.0
    DO ii=1,eqs
          L(ii,0)=1.0
    END DO

    DO ii=1,eqs
       DO i=0,m
          DO k=0,i
             U(ii,i)=U(ii,i)+y(ii,k)*L(ii,i-k)
          END DO
       END DO
    END DO


    ! Evaluating upper and lower polynomials in horner form
    numsum=0.0
    denomsum=0.0

    DO ii=1,eqs
       DO i=m,1,-1
          numsum(ii)=(numsum(ii)+U(ii,i))*h
          denomsum(ii)=(denomsum(ii)+L(ii,i))*h
       END DO
       numsum(ii)=numsum(ii)+U(ii,0)
       denomsum(ii)=denomsum(ii)+L(ii,0)
    END DO

    DO ii=1,eqs
         x(ii)=numsum(ii)/denomsum(ii)
    END DO
    !x(3)=y(1,0)**(7.0D0/3.0D0)
    !x(4)=1.0D0/y(1,0)




   RETURN
END FUNCTION padstep




!********************************************************************
!____________________________________________________________________
FUNCTION padterm(y,h,ii) RESULT(x)

   IMPLICIT NONE

   REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: y
   REAL(KIND=8), INTENT(IN) :: h
   INTEGER, INTENT(IN)::ii


   INTEGER::   m,k,j,i
   REAL(KIND=8)::  numsum, denomsum,x
   REAL(KIND=8), DIMENSION(0:((SIZE(y,2)-1))/2 ):: L,U
   REAL(KIND=8), DIMENSION(0:SIZE(y,2)):: yi
   REAL(KIND=8), DIMENSION(1:((SIZE(y,2)-1))/2 ,1:((SIZE(y,2)-1))/2 ) :: A
   REAL(KIND=8), DIMENSION(1:((SIZE(y,2)-1))/2 ):: b, z



   m=SIZE(y,2)
   m=(m-1)/2

   DO j=0,2*m
        yi(j)=y(ii,j+1)
        !WRITE(*,*)'yi(',j,')=',yi(j)
   END DO

    ! NOW START PADE

   !FILL MATRIX AND RHS FOR PADE

   DO i=1,m
       DO j=1,m
           A(i,j)=yi(m+i-j)
       END DO
       b(i)=-yi(m+i)
   END DO


   z=gauss(A,b)

   DO i=1,m
      L(i)=z(i)
   END DO
   !SET NUMER COEF USING CAUCHY PRODUCT

    U=0.0

    L(0)=1.0

    DO i=0,m
       DO k=0,i
             U(i)=U(i)+yi(k)*L(i-k)
       END DO
    END DO

    ! Evaluating upper and lower polynomials in horner form
    numsum=0.0
    denomsum=0.0


    DO i=m,1,-1
          numsum=(numsum+U(i))*h
          denomsum=(denomsum+L(i))*h
    END DO
    numsum=numsum+U(0)
    denomsum=denomsum+L(0)



    x=numsum/denomsum

    !x(3)=y(1,0)**(7.0D0/3.0D0)
    !x(4)=1.0D0/y(1,0)

    RETURN
END FUNCTION padterm




FUNCTION cp(a,ii,b,jj,i) RESULT(x)

   IMPLICIT NONE

   REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: a
   REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: b
   INTEGER, INTENT(IN)::ii, jj, i
   REAL(KIND=8)::x
   INTEGER:: k



   x=0.0

   DO k=0,i
     x=x+a(ii,k+1)*b(jj,i-k+1)
   END DO





    RETURN
END FUNCTION cp


!********************************************************************
!____________________________________________________________________
FUNCTION picterm(y,h,ii) RESULT(x)

   IMPLICIT NONE

   REAL(KIND=8), DIMENSION(:,:), INTENT(IN) :: y
   REAL(KIND=8), INTENT(IN) :: h
   INTEGER, INTENT(IN):: ii


   INTEGER::   m,i
   REAL(KIND=8)::  x


   m=SIZE(y,2)-1


   x=0.0
   DO i=m,1,-1
          x=(x+y(ii,i+1))*h
    END DO
    x=x+y(ii,1)


    RETURN
END FUNCTION picterm



END MODULE pade2