scaling_function.F90

Go to the documentation of this file.

!! Copyright (C) 2011 X. Andrade
!! Copyright (C) Luigi Genovese, Thierry Deutsch, CEA Grenoble, 2006
!!
!! This program is free software; you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 2, or (at your option)
!! any later version.
!!
!! This program is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with this program; if not, write to the Free Software
!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
!! 02110-1301, USA.
!!
 
#include <global.h>
 
module scaling_function_oct_m
  use, intrinsic :: iso_fortran_env
  use blas_oct_m
  use messages_oct_m
 
  implicit none
 
  private
 
  public ::            &
    scaling_function,  &
    scf_recursion
 
  type lazy_8_filter_t
    integer :: m = 10
    real(real64) :: ch(-10:10)
    real(real64) :: cg(-10:10)
    real(real64) :: cht(-10:10)
    real(real64) :: cgt(-10:10)
  end type lazy_8_filter_t
 
contains
 
  ! ---------------------------------------------------------
  pure function lazy_8() result(filters)
    type(lazy_8_filter_t) :: filters
    integer :: i
 
    filters%ch = 0._real64
    filters%cht = 0._real64
    filters%cg = 0._real64
    filters%cgt = 0._real64
 
    ! The normalization is chosen such that a constant function remains the
    ! same constant on each level of the transform.
    filters%ch(-7) = -5._real64/2048._real64
    filters%ch(-5) = 49._real64/2048._real64
    filters%ch(-3) = -245._real64/2048._real64
    filters%ch(-1) = 1225._real64/2048._real64
    filters%ch( 0) = 1._real64
    filters%ch( 1) = 1225._real64/2048._real64
    filters%ch( 3) = -245._real64/2048._real64
    filters%ch( 5) = 49._real64/2048._real64
    filters%ch( 7) = -5._real64/2048._real64
    filters%cht(0) = 1._real64
 
    ! g coefficients from h coefficients
    do i = -filters%m, filters%m - 1
      filters%cg(i + 1) = filters%cht(-i) * real((2*modulo(i,2)-1), real64)
      filters%cgt(i + 1) = filters%ch(-i) * real((2*modulo(i,2)-1), real64)
    end do
  end function lazy_8
 
  !!****h* BigDFT/scaling_function
  !! NAME
  !!   scaling_function
  !!
  !! FUNCTION
  !!   Calculate the values of a scaling function in real uniform grid
  !!
  !! SOURCE
  !!
  subroutine scaling_function(itype,nd,nrange,a,x)
    integer, intent(in) :: itype
    integer, intent(in) :: nd
    integer, intent(out) :: nrange
    real(real64), dimension(0:nd), intent(out) :: a,x
 
    real(real64), dimension(:), allocatable :: y
    integer :: i,nt,ni
 
    !Only itype=8,14,16,20,24,30,40,50,60,100
    select case (itype)
    case (8)
      !O.K.
    case default
      message(1) = "Only interpolating functions 8, 14, 16, 20, 24, 30, 40, 50, 60, 100."
      call messages_fatal(1)
    end select
!!$  write(unit=*,fmt="(1x,a,i0,a)") &
!!$       "Use interpolating scaling functions of ",itype," order"
 
    !Give the range of the scaling function
    !from -itype to itype
    ni=2*itype
    nrange = ni
    allocate(y(0:nd))
 
    ! plot scaling function
    x = 0.0_8
    y = 0.0_8
 
    nt=ni
    x(nt/2-1)=1._real64
    loop1: do
      nt=2*nt
      ! write(6,*) 'nd,nt',nd,nt
      select case (itype)
      case (8)
        call back_trans_8(nd,nt,x,y)
      end select
      call blas_copy(nt, y(0), 1, x(0) ,1)
      if (nt == nd) then
        exit loop1
      end if
    end do loop1
 
    !open (unit=1,file='scfunction',status='unknown')
    do i=0,nd
      a(i) = 1._real64 * i * ni / nd - (.5_real64 * ni - 1._real64)
      !write(1,*) 1._real64*i*ni/nd-(.5_real64*ni-1._real64),x(i)
    end do
    !close(1)
 
    deallocate(y)
  end subroutine scaling_function
  !!***
 
  !!****h* BigDFT/scf_recursion
  !! NAME
  !!   scf_recursion
  !!
  !! FUNCTION
  !!   Do iterations to go from p0gauss to pgauss
  !!   order interpolating scaling function
  !!
  !! SOURCE
  !!
  subroutine scf_recursion(itype,n_iter,n_range,kernel_scf,kern_1_scf)
    integer, intent(in) :: itype,n_iter,n_range
    real(kind=8), intent(inout) :: kernel_scf(-n_range:n_range)
    real(kind=8), intent(out) :: kern_1_scf(-n_range:n_range)
 
    !Only itype=8,14,16,20,24,30,40,50,60,100
    select case (itype)
    case (8)
      !O.K.
    case default
      message(1) = "Only interpolating functions 8, 14, 16, 20, 24, 30, 40, 50, 60, 100."
      call messages_fatal(1)
    end select
 
    select case (itype)
    case (8)
      call scf_recursion_8(n_iter,n_range,kernel_scf,kern_1_scf)
    end select
 
  end subroutine scf_recursion
  !!***
 
  !!****h* BigDFT/back_trans_8
  !! NAME
  !!   back_trans_8
  !!
  !! FUNCTION
  !!
  !! SOURCE
  !!
  ! backward wavelet transform
  ! nd: length of data set
  ! nt length of data in data set to be transformed
  ! m filter length (m has to be even!)
  ! x input data, y output data
  subroutine back_trans_8(nd,nt,x,y)
    integer, intent(in) :: nd,nt
    real(kind=8), intent(in) :: x(0:nd-1)
    real(kind=8), intent(out) :: y(0:nd-1)
 
    type(lazy_8_filter_t) :: filters
    integer :: i,j,ind
 
    filters = lazy_8()
 
    do i=0,nt/2-1
      y(2*i+0)=0._real64
      y(2*i+1)=0._real64
 
      do j=-filters%m/2,filters%m/2-1
 
        ! periodically wrap index if necessary
        ind=i-j
        do
          if (ind < 0) then
            ind=ind+nt/2
            cycle
          end if
          if (ind >= nt/2) then
            ind=ind-nt/2
            cycle
          end if
          exit
        end do
 
        y(2*i+0)=y(2*i+0) + filters%ch(2*j-0)*x(ind)+filters%cg(2*j-0)*x(ind+nt/2)
        y(2*i+1)=y(2*i+1) + filters%ch(2*j+1)*x(ind)+filters%cg(2*j+1)*x(ind+nt/2)
      end do
 
    end do
 
  end subroutine back_trans_8
  !!***
 
  !!****h* BigDFT/scf_recursion_8
  !! NAME
  !!   scf_recursion_8
  !!
  !! FUNCTION
  !!   Do iterations to go from p0gauss to pgauss
  !!   8th-order interpolating scaling function
  !!
  !! SOURCE
  !!
  subroutine scf_recursion_8(n_iter,n_range,kernel_scf,kern_1_scf)
    integer, intent(in) :: n_iter,n_range
    real(kind=8), intent(inout) :: kernel_scf(-n_range:n_range)
    real(kind=8), intent(out) :: kern_1_scf(-n_range:n_range)
 
    type(lazy_8_filter_t) :: filters
    real(kind=8) :: kern,kern_tot
    integer :: i_iter,i,j,ind
 
    filters = lazy_8()
 
    !Start the iteration to go from p0gauss to pgauss
    loop_iter_scf: do i_iter=1,n_iter
      kern_1_scf(:) = kernel_scf(:)
      kernel_scf(:) = 0._real64
      loop_iter_i: do i=0,n_range
        kern_tot = 0._real64
        do j=-filters%m,filters%m
          ind = 2*i-j
          if (abs(ind) > n_range) then
            kern = 0._real64
          else
            kern = kern_1_scf(ind)
          end if
          kern_tot = kern_tot + filters%ch(j)*kern
        end do
        if (abs(kern_tot) <= 1d-150) then
          !zero after (be sure because strictly == 0._real64)
          exit loop_iter_i
        else
          kernel_scf( i) = 0.5_real64*kern_tot
          kernel_scf(-i) = kernel_scf(i)
        end if
      end do loop_iter_i
    end do loop_iter_scf
  end subroutine scf_recursion_8
  !!***
 
end module scaling_function_oct_m
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: