mesh_function.F90

Go to the documentation of this file.

!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
!!
!! 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 mesh_function_oct_m
  use blas_oct_m
  use comm_oct_m
  use debug_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use lalg_basic_oct_m
  use loct_math_oct_m
  use mesh_oct_m
  use messages_oct_m
  use mpi_oct_m
  use par_vec_oct_m
  use profiling_oct_m
  use qshep_oct_m
  use quickrnd_oct_m
  use splines_oct_m
 
  implicit none
 
  private
  public ::                &
    dmf_integrate,         &
    zmf_integrate,         &
    dmf_remove_average,    &
    zmf_remove_average,    &
    dmf_dotp,              &
    zmf_dotp,              &
    dmf_nrm2,              &
    zmf_nrm2,              &
    dmf_moment,            &
    zmf_moment,            &
    dmf_random,            &
    zmf_random,            &
    dmf_interpolate_points,&
    zmf_interpolate_points,&
    mf_surface_integral,   &
    mf_line_integral,      &
    dmf_dotp_aux,          &
    zmf_dotp_aux,          &
    dmf_multipoles,        &
    zmf_multipoles,        &
    dmf_dipole,            &
    zmf_dipole,            &
    mesh_init_mesh_aux,    &
    dmf_remove_average_aux,&
    zmf_remove_average_aux,&
    dmf_dotu_aux,          &
    zmf_dotu_aux,          &
    dmf_nrm2_aux,          &
    zmf_nrm2_aux,          &
    dmf_normalize,         &
    zmf_normalize,         &
    zmf_fix_phase
 
  ! These variables are to be used by the "distdot" function, that is outside the module
  ! but inside this file.
  ! FIXME: This is very ugly, at least these values should be set by a function.
  public :: mesh_aux
  logical, public :: sp_parallel
  integer, public :: sp_np, sp_dim, sp_st1, sp_st2, sp_kp1, sp_kp2
  integer, public :: sp_distdot_mode
  type(mpi_grp_t), public :: sp_grp
 
  interface mf_surface_integral
    module procedure dmf_surface_integral_scalar, dmf_surface_integral_vector, &
      zmf_surface_integral_scalar, zmf_surface_integral_vector
  end interface mf_surface_integral
 
  interface mf_line_integral
    module procedure dmf_line_integral_scalar, dmf_line_integral_vector, &
      zmf_line_integral_scalar, zmf_line_integral_vector
  end interface mf_line_integral
 
  interface dmf_dotp
    module procedure dmf_dotp_1, dmf_dotp_2
  end interface dmf_dotp
 
  interface zmf_dotp
    module procedure zmf_dotp_1, zmf_dotp_2
  end interface zmf_dotp
 
  interface dmf_nrm2
    module procedure dmf_nrm2_1, dmf_nrm2_2
  end interface dmf_nrm2
 
  interface zmf_nrm2
    module procedure zmf_nrm2_1, zmf_nrm2_2
  end interface zmf_nrm2
 
  interface dmf_integrate
    module procedure dmf_integrate_1, dmf_integrate_2
  end interface dmf_integrate
 
  interface zmf_integrate
    module procedure zmf_integrate_1, zmf_integrate_2
  end interface zmf_integrate
 
 
  class(mesh_t), pointer :: mesh_aux => null()
 
contains
 
  subroutine mesh_init_mesh_aux(mesh)
    class(mesh_t), target, intent(in) :: mesh
 
    push_sub(mesh_init_mesh_aux)
    mesh_aux => mesh
 
    pop_sub(mesh_init_mesh_aux)
  end subroutine mesh_init_mesh_aux
 
  subroutine zmf_fix_phase(mesh, ff)
    class(mesh_t),     intent(in)     :: mesh
    complex(real64),   intent(inout)  :: ff(:)
 
    real(real64), allocatable :: ff_re(:), ff_im(:)
    real(real64) :: dot_re, dot_im, dot_reim, sol1, sol2, maxv, maxvloc
    real(real64) :: theta, ctheta, stheta
    integer :: ip, ipmax
    complex(real64)   :: zval
 
    push_sub(zmf_fix_phase)
 
    safe_allocate(ff_re(1:mesh%np))
    safe_allocate(ff_im(1:mesh%np))
 
    ff_re(:) = real(ff(:), real64)
    ff_im(:) = aimag(ff(:))
 
    dot_re = dmf_nrm2(mesh, ff_re)**2
    dot_im = dmf_nrm2(mesh, ff_im)**2
    dot_reim = dmf_dotp(mesh, ff_re, ff_im)
 
    ! We have tan(2\theta) = 2 (\int \psi_R \psi_I)/( \int \psi_R^2-\psi_I^2)
    if(dot_re+dot_im > 1e-8_real64) then
      if(m_two*dot_reim/(dot_re-dot_im) > 1e-8_real64) then
        theta = m_half * atan(m_two*dot_reim/(dot_re-dot_im))
 
        ! Get the maximum for the real part
        sol1 = cos(theta)**2 * dot_re + sin(theta)**2 * dot_im + m_two*sin(theta)*cos(theta)*dot_reim
        sol2 = sin(theta)**2 * dot_re + cos(theta)**2 * dot_re - m_two*sin(theta)*cos(theta)*dot_reim
        if(sol2>sol1) theta = theta + m_half*m_pi
 
        ctheta = cos(theta)
        stheta = sin(theta)
 
        !We do the rotation
        !$omp parallel do
        do ip = 1, mesh%np
          ff(ip) = cmplx(ctheta*ff_re(ip)-stheta*ff_im(ip), stheta*ff_re(ip)+ctheta*ff_im(ip), real64)
        end do
 
      else ! We do not know what is the direction here
        !    we set it using the maximum value of the function from rank 0
        ipmax = 0
        maxv = m_zero
        do ip = 1, mesh%np
          maxvloc = real(ff(ip)*conjg(ff(ip)), real64)
          if(maxvloc > maxv) then
            ipmax = ip
            maxv = maxvloc
          end if
        end do
        zval = ff(ipmax)
 
        !We want to use the same phase on all processors
        if(mesh%parallel_in_domains) then
          call mesh%mpi_grp%bcast(zval, 1, mpi_double_complex, 0)
        end if
 
        call lalg_scal(mesh%np, abs(zval)/zval, ff)
 
      end if
    else
      message(1) = "Internal error fixing phase: mesh function has zero norm."
      call messages_fatal(1)
    end if
 
    safe_deallocate_a(ff_re)
    safe_deallocate_a(ff_im)
 
    pop_sub(zmf_fix_phase)
  end subroutine zmf_fix_phase
 
#include "undef.F90"
#include "real.F90"
#include "mesh_function_inc.F90"
 
#include "undef.F90"
#include "complex.F90"
#include "mesh_function_inc.F90"
 
end module mesh_function_oct_m
 
#ifdef HAVE_SPARSKIT
 
real(real64) function distdot(n, x, ix, y, iy)
  use comm_oct_m
  use global_oct_m
  use, intrinsic :: iso_fortran_env
  use messages_oct_m
  use mesh_function_oct_m
  use profiling_oct_m
 
  implicit none
 
  integer,     intent(in) :: n
  real(real64),    intent(in) :: x(n)
  integer,     intent(in) :: ix
  real(real64),    intent(in) :: y(n)
  integer,     intent(in) :: iy
 
  integer :: j, ik, ist, idim, k
 
  ! SPARSKIT only calls this function with ix, iy = 1 i.e. no stride.
  assert(ix == 1)
  assert(iy == 1)
 
  select case (sp_distdot_mode)
  case (1)
    distdot = dmf_dotp_aux(x(1:n/2), y(1:n/2)) + dmf_dotp_aux(x(n/2+1:n), y(n/2+1:n))
 
  case (2)
    distdot = m_zero
    j = 1
    k = n/2+1
    do ik = sp_kp1, sp_kp2
      do ist = sp_st1, sp_st2
        do idim = 1, sp_dim
          distdot = distdot + dmf_dotp_aux(x(j: j+sp_np-1), y(j:j+sp_np-1))
          distdot = distdot + dmf_dotp_aux(x(k: k+sp_np-1), y(k:k+sp_np-1))
          j = j + sp_np
          k = k + sp_np
        end do
      end do
    end do
    if (sp_parallel) call comm_allreduce(sp_grp, distdot)
 
  case (3)
    distdot = m_zero
    j = 1
    k = n/2+1
    do ik = sp_kp1, sp_kp2
      do ist = sp_st1, sp_st2
        do idim = 1, sp_dim
          distdot = distdot + dmf_dotp_aux(x(j: j+sp_np-1), y(j:j+sp_np-1))
          distdot = distdot + dmf_dotp_aux(x(k: k+sp_np-1), y(k:k+sp_np-1))
          j = j + sp_np
          k = k + sp_np
        end do
      end do
    end do
    do ik = sp_kp1, sp_kp2
      do ist = sp_st1, sp_st2
        do idim = 1, sp_dim
          distdot = distdot + dmf_dotp_aux(x(j: j+sp_np-1), y(j:j+sp_np-1))
          distdot = distdot + dmf_dotp_aux(x(k: k+sp_np-1), y(k:k+sp_np-1))
          j = j + sp_np
          k = k + sp_np
        end do
      end do
    end do
    if (sp_parallel) call comm_allreduce(sp_grp, distdot)
 
  end select
 
end function distdot
 
#endif /* HAVE_SPARSKIT */
 
!! Local Variables:
!! mode: f90
!! coding: utf-8
!! End: