radial_basis_function.f90

4总常用的插值程序的Fortran90源代码

3800  format(7x,'x',9x,'numer',8x,'anal',9x,'err', &
  10x,'du/dx (numer and anal)')
  write ( *,'(1p6e13.5)')z,sol,anal,err,dsdx,adfdx

  go to 20

end
function ffun ( i, x )
!
!*******************************************************************************
!
!! FFUN returns the value of one of six different functions of X.
!
!
!  Discussion:
!
!    The user chooses the function by specifying an input argument.
!
  real ffun
  integer i
  real x
!
  if ( i == 1 ) then
    ffun = x**2
  else if ( i == 2 ) then
    ffun = x * ( x - 1.0 )**2
  else if ( i == 3 ) then
    ffun = sin ( 2.0 * 3.14159265358979323844 * x )
  else if ( i == 4 ) then
    ffun = exp ( - x )
  else if ( i == 5 ) then
    ffun = abs ( x - 0.5 )
  else if ( i == 6 ) then
    ffun = atan ( 10.0 * ( x - 0.5 ) )
  else
    ffun = 0.0E+00
  end if

  return
end
function dfundx(i,x)
!
!*******************************************************************************
!
!! DFUNDX returns the derivative of one of six different functions of X.
!
!
  real del
  real dfun
  real dfundx
  integer i
  real pi
  real pi2
  real x
!
  if ( i == 1 ) then
    dfundx=2*x
  else if ( i == 2 ) then
    dfundx=(x-1)**2 + x*2*(x-1)
  else if ( i == 3 ) then
    pi=3.14159265358979323844
    pi2=2.0*pi
    dfundx=pi2*cos(pi2*x)
  else if ( i == 4 ) then
    dfundx=-exp(-x)
  else if ( i == 5 ) then
    dfundx=-1
    if(x.ge.0.5)dfundx=1
  else if ( i == 6 ) then
    del=0.001
    dfun=atan(10.0*(x+del-0.5))-atan(10.0*(x-del-0.5))
    dfundx=dfun/(2.0*del)
  else
    dfundx = 0.0E+00
  end if

  return
end
function d2fundx(i,x)
!
!*******************************************************************************
!
!! D2FUNDX returns the second derivative of one of six different functions of X.
!
!
  real d2
  real d2fundx
  real del
  integer i
  real, parameter :: pi = 3.1415926535897932384
  real pi2
  real x
!
  if ( i == 1 ) then
    d2fundx=2
  else if ( i == 2 ) then
    d2fundx=2*(x-1) + 2*(x-1) + x*2
  else if ( i == 3 ) then
    pi2=2.0*pi
    d2fundx=-pi2**2*sin(pi2*x)
  else if ( i == 4 ) then
    d2fundx=exp(-x)
  else if ( i == 5 ) then
    d2fundx=0.0
    if(x.eq.0.5)d2fundx=1e6
  else if ( i == 6 ) then
    del=0.001
    d2  =atan(10.0*(x+del-0.5))-2*atan(10.0*(x-0.5))+atan(10.0*(x-del-0.5))
    d2fundx=d2/del**2
  else
    d2fundx = 0.0
  end if

  return
end
subroutine gnuplot
!
!*******************************************************************************
!
!! GNUPLOT writes out a file of commands to the GNUPLOT graphics package.
!
  integer lgp
!
  lgp = 12

  open ( unit = lgp, file = 'rbf_gnu_commands.txt', status = 'replace' ) 

  write ( lgp, * ) 'set title "mq-lsh method interpolation errors"'
  write ( lgp, * ) 'set xlabel "x"'
  write ( lgp, * ) 'set ylabel "f, df/dx error"'
  write ( lgp, * ) 'plot "mqle.dat" u 2:5 w l 1,"mqle.dat" u 2:8 w l -1'
  write ( lgp, * ) 'pause -1'
  write ( lgp, * ) 'set title "mq-ls : f - original and interpolated function"'
  write ( lgp, * ) 'set ylabel "f"'
  write ( lgp, * ) 'plot "mqle.dat" u 2:3 w l 1,"mqle.dat" u 2:4 w l 3,\\'
  write ( lgp, * ) '"mqle2.dat" u 2:3 w p 3'
  write ( lgp, * ) 'pause -1'
  write ( lgp, * ) 'set title "mq-ls : df/dx - original and interpolated function"'
  write ( lgp, * ) 'set ylabel "df/dx"'
  write ( lgp, * ) 'plot "mqle.dat" u 2:6 w l 1,"mqle.dat" u 2:7 w l 5,\\'
  write ( lgp, * ) '"mqle2.dat" u 2:6 w p 3'
  write ( lgp, * ) 'pause -1'
  write ( lgp, * ) 'q'

  close ( unit = lgp )

  write ( *, * ) 'starting gnuplot. hit return for new figure'

  call system ( 'gnuplot rbf_gnu_commands.txt' )

  return
end
subroutine hlsrec ( na, x, fvec, iflag )
!
!*******************************************************************************
!
!! HLSREC calculates the functions at X.
!
!
  integer, parameter :: nmax = 501
  integer, parameter :: nxtm = 1001
!
  integer na
!
  real aaa(nmax,nmax)
  real asa(nmax,nmax)
  real b(nmax)
  real c(nmax)
  real cond
  real d2fdx2
  real d2fundx
  real dfdx
  real dfundx
  real ffun
  real fun
  real fvec(na)
  real g
  real gx
  real gxc
  real g2x
  real hxt
  integer i
  integer iflag
  integer ifun
  integer ipr
  integer iprec
  integer ipvt(nmax)
  integer j
  integer k
  integer m
  integer mode
  integer n
  integer, save :: nit = 0
  integer nn
  integer nqq
  integer nxt
  real resid
  real uxk
  real w(nmax)
  real x(na)
  real xk(nmax)
  real xt(nxtm)
  real xx2
  real y(nmax)
  real z
!
  common /i_common/ ifun, m, mode, nn, nqq, nxt
  common /r_common/ aaa, asa, b, hxt, xk, xt
!
!  These BLANKING STATEMENT FUNCTIONS SHOULD NEVER NEVER BE USED!
!
! declaration operator-functions
!
  fun(z)    = ffun(ifun,z)
  dfdx(z)   = dfundx(ifun,z)
  d2fdx2(z) = d2fundx(ifun,z)

  g(j,z)  = sqrt((z-xk(j))**2+c(j)**2)

  gx(j,z)  = (z-xk(j))/sqrt((z-xk(j))**2+c(j)**2)
!
! gxc(j,z)  = -c(j)*(z-xk(j))/sqrt(((z-xk(j))**2+c(j)**2)**3)
!
  gxc(j,z)  = (z-xk(j))/sqrt(((z-xk(j))**2+c(j)**2)**3)

  g2x(j,z) =      1.0/sqrt((z-xk(j))**2+c(j)**2) &
     - (z-xk(j))**2/sqrt(((z-xk(j))**2+c(j)**2)**3)


  iprec = 1
  ipr = 0
  n = m
  nn = n
  nit = nit+1
!
!  Form a matrix using the current x = [x,c]
!
  c(1:n) = x(1:n)

  do j = 1, n
    do i = 1, n
      aaa(i,j) = g(j,xk(i))
    end do
  end do

  do i = 1, n
    y(i) = fun ( xk(i) )
  end do
!
!  Preconditioning
!
  if ( iprec /= 0 )then

    do i = 1, n
      b(i) = aaa(i,i)
    end do

    do j = 1, n
      do i = 1, n
        aaa(i,j) = aaa(i,j) / b(i)
      end do
    end do

    y(1:n) = y(1:n) / b(1:n)

  end if

  call decomp ( nmax, n, aaa, cond, ipvt, w )

! if ( nit<=20) then
!   write(45,1300) n,nxt,cond
! end if

1300  format(' hsolvc mqls: knots',i5,' points',i5, ' cond  = ',1pe15.5)

  call solve ( nmax, n, aaa, y, ipvt )
!
!  c_j related matrix block
!
  do i = 1, n

    fvec(i) = 0.0

    do k = 1, n
      uxk = 0.0
      do j = 1, n
        uxk = uxk + y(j) * gx(j,xk(k))
      end do

      fvec(i) = fvec(i) + y(i) * gxc(i,xk(k)) * (uxk-dfdx(xk(k)))

    end do

  end do

  resid = 0.0

  xx2 = sqrt ( sum ( x(1:n)**2 ) )
  resid = sqrt ( sum ( fvec(1:n)**2 ) )

  if ( nit <= 10 .or. mod ( nit, 10 ) == 0 ) then
    write(6,100)nit,na,m,xx2,resid,cond
  end if

100   format('hlsrec: nit,nn,m = ',3i4,' |x| =',1pe11.3,' resid=',1pe10.2 &
       ,' cond = ',1pe10.2)

  return
end
subroutine decomp ( ndim, n, a, cond, ipvt, work )
!
!*******************************************************************************
!
!! DECOMP computes a PLU decomposition of a matrix.
!
!
  integer ndim
  integer n
!
  real a(ndim,n)
  real anorm
  real cond
  real ek
  integer i
  integer ipvt(n)
  integer j
  integer k
  integer kb
  integer m
  real t
  real work(n)
  real ynorm
  real znorm
!
  ipvt(n) = 1

  if ( n == 1 ) then
    if ( a(1,1) /= 0.0 ) then
      cond = 1.0
    else
      cond = huge ( cond )
    end if
    return
  end if

  anorm = 0.0
  do j = 1, n
    t = sum ( abs ( a(1:n,j) ) )
    anorm = max ( anorm, t )
  end do

  do k = 1, n-1

    m = k
    do i = k+1,n
      if(abs(a(i,k)) > abs(a(m,k)))m = i
    end do

    ipvt(k) = m
    if ( m /= k ) ipvt(n) = -ipvt(n)

    call r_swap ( a(m,k), a(k,k) )

    if ( a(k,k) /= 0.0 ) then

      do i = k+1,n
        a(i,k) = -a(i,k) / a(k,k)
      end do

      do j = k+1,n

        call r_swap ( a(m,j), a(k,j) )

        if ( a(k,j) /= 0.0 ) then
          do i = k+1, n
            a(i,j) = a(i,j) + a(i,k) * a(k,j)
          end do
        end if

     end do

    end if

  end do

  do k = 1,n

    t = 0.0
    do i = 1, k-1
      t = t + a(i,k) * work(i)
    end do

    ek = 1.0
    if ( t < 0.0 ) ek = -1.0

    if ( a(k,k) == 0.0 ) then
      cond = huge ( cond )
      return
    end if

    work(k) = -(ek+t) / a(k,k)
  end do

  do k = n-1, 1, -1

    t = 0.0
    do i = k+1, n
      t = t + a(i,k) * work(k)
    end do

    work(k) = t
    m = ipvt(k)

    if ( m /= k ) then
      call r_swap ( work(m), work(k) )
    end if

  end do

  ynorm = sum ( abs ( work(1:n) ) )

  call solve ( ndim, n, a, work, ipvt )

  znorm = sum ( abs ( work(1:n) ) )

  cond = anorm*znorm/ynorm

  cond = max ( cond, 1.0 )

  return
end
subroutine solve ( ndim, n, a, b, ipvt )
!
!*******************************************************************************
!
!! SOLVE solves a linear system whose matrix has been decomposed by DECOMP.
!
!
  integer ndim
  integer n
!
  real a(ndim,n)
  real b(n)
  integer i
  integer ipvt(n)
  integer k
  integer m
!
  if ( n == 1 ) then
    b(1) = b(1) / a(1,1)
    return
  end if

  do k = 1, n-1
    m = ipvt(k)
    call r_swap ( b(m), b(k) )
    do i = k+1, n
      b(i) = b(i) + a(i,k) * b(k)
    end do
  end do

  do k = n, 2, -1

    b(k) = b(k) / a(k,k)

    do i = 1, k-1
      b(i) = b(i) - a(i,k) * b(k)
    end do

  end do

  b(1) = b(1) / a(1,1)

  return
end