Octopus
mesh_function.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
24 use blas_oct_m
25 use comm_oct_m
26 use debug_oct_m
27 use global_oct_m
28 use, intrinsic :: iso_fortran_env
31 use mesh_oct_m
33 use mpi_oct_m
36 use qshep_oct_m
39
40 implicit none
41
42 private
43 public :: &
46 dmf_dotp, &
47 zmf_dotp, &
48 dmf_nrm2, &
49 zmf_nrm2, &
50 dmf_moment, &
51 zmf_moment, &
52 dmf_random, &
53 zmf_random, &
62 dmf_dipole, &
63 zmf_dipole, &
72
73 ! These variables are to be used by the "distdot" function, that is outside the module
74 ! but inside this file.
75 ! FIXME: This is very ugly, at least these values should be set by a function.
76 public :: mesh_aux
77 logical, public :: sp_parallel
78 integer, public :: sp_np, sp_dim, sp_st1, sp_st2, sp_kp1, sp_kp2
79 integer, public :: sp_distdot_mode
80 type(mpi_grp_t), public :: sp_grp
81
82 interface mf_surface_integral
85 end interface mf_surface_integral
86
87 interface mf_line_integral
90 end interface mf_line_integral
91
92 interface dmf_dotp
93 module procedure dmf_dotp_1, dmf_dotp_2
94 end interface dmf_dotp
95
96 interface zmf_dotp
97 module procedure zmf_dotp_1, zmf_dotp_2
98 end interface zmf_dotp
99
100 interface dmf_nrm2
101 module procedure dmf_nrm2_1, dmf_nrm2_2
102 end interface dmf_nrm2
103
104 interface zmf_nrm2
105 module procedure zmf_nrm2_1, zmf_nrm2_2
106 end interface zmf_nrm2
107
108 interface dmf_integrate
109 module procedure dmf_integrate_1, dmf_integrate_2
110 end interface dmf_integrate
111
112 interface zmf_integrate
113 module procedure zmf_integrate_1, zmf_integrate_2
114 end interface zmf_integrate
115
118 class(mesh_t), pointer :: mesh_aux => null()
119
120contains
121
125 subroutine mesh_init_mesh_aux(mesh)
126 class(mesh_t), target, intent(in) :: mesh
127
128 push_sub(mesh_init_mesh_aux)
129 mesh_aux => mesh
130
131 pop_sub(mesh_init_mesh_aux)
132 end subroutine mesh_init_mesh_aux
133
141 subroutine zmf_fix_phase(mesh, ff)
142 class(mesh_t), intent(in) :: mesh
143 complex(real64), intent(inout) :: ff(:)
144
145 real(real64), allocatable :: ff_re(:), ff_im(:)
146 real(real64) :: dot_re, dot_im, dot_reim, sol1, sol2, maxv, maxvloc
147 real(real64) :: theta, ctheta, stheta
148 integer :: ip, ipmax
149 complex(real64) :: zval
150
151 push_sub(zmf_fix_phase)
152
153 safe_allocate(ff_re(1:mesh%np))
154 safe_allocate(ff_im(1:mesh%np))
155
156 ff_re(:) = real(ff(:), real64)
157 ff_im(:) = aimag(ff(:))
158
159 dot_re = dmf_nrm2(mesh, ff_re)**2
160 dot_im = dmf_nrm2(mesh, ff_im)**2
161 dot_reim = dmf_dotp(mesh, ff_re, ff_im)
162
163 ! We have tan(2\theta) = 2 (\int \psi_R \psi_I)/( \int \psi_R^2-\psi_I^2)
164 if(dot_re+dot_im > 1e-8_real64) then
165 if(m_two*dot_reim/(dot_re-dot_im) > 1e-8_real64) then
166 theta = m_half * atan(m_two*dot_reim/(dot_re-dot_im))
167
168 ! Get the maximum for the real part
169 sol1 = cos(theta)**2 * dot_re + sin(theta)**2 * dot_im + m_two*sin(theta)*cos(theta)*dot_reim
170 sol2 = sin(theta)**2 * dot_re + cos(theta)**2 * dot_re - m_two*sin(theta)*cos(theta)*dot_reim
171 if(sol2>sol1) theta = theta + m_half*m_pi
173 ctheta = cos(theta)
174 stheta = sin(theta)
176 !We do the rotation
177 !$omp parallel do
178 do ip = 1, mesh%np
179 ff(ip) = cmplx(ctheta*ff_re(ip)-stheta*ff_im(ip), stheta*ff_re(ip)+ctheta*ff_im(ip), real64)
180 end do
181
182 else ! We do not know what is the direction here
183 ! we set it using the maximum value of the function from rank 0
184 ipmax = 0
185 maxv = m_zero
186 do ip = 1, mesh%np
187 maxvloc = real(ff(ip)*conjg(ff(ip)), real64)
188 if(maxvloc > maxv) then
189 ipmax = ip
190 maxv = maxvloc
191 end if
192 end do
193 zval = ff(ipmax)
194
195 !We want to use the same phase on all processors
196 if(mesh%parallel_in_domains) then
197 call mesh%mpi_grp%bcast(zval, 1, mpi_double_complex, 0)
198 end if
199
200 call lalg_scal(mesh%np, abs(zval)/zval, ff)
202 end if
203 else
204 message(1) = "Internal error fixing phase: mesh function has zero norm."
206 end if
207
208 safe_deallocate_a(ff_re)
209 safe_deallocate_a(ff_im)
210
212 end subroutine zmf_fix_phase
213
214#include "undef.F90"
215#include "real.F90"
216#include "mesh_function_inc.F90"
217
218#include "undef.F90"
219#include "complex.F90"
220#include "mesh_function_inc.F90"
221
222end module mesh_function_oct_m
223
224#ifdef HAVE_SPARSKIT
225
229real(real64) function distdot(n, x, ix, y, iy)
230 use comm_oct_m
231 use global_oct_m
232 use, intrinsic :: iso_fortran_env
236
237 implicit none
238
239 integer, intent(in) :: n
240 real(real64), intent(in) :: x(n)
241 integer, intent(in) :: ix
242 real(real64), intent(in) :: y(n)
243 integer, intent(in) :: iy
244
245 integer :: j, ik, ist, idim, k
246
247 ! SPARSKIT only calls this function with ix, iy = 1 i.e. no stride.
248 assert(ix == 1)
249 assert(iy == 1)
250
251 select case (sp_distdot_mode)
252 case (1)
253 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))
254
255 case (2)
257 j = 1
258 k = n/2+1
259 do ik = sp_kp1, sp_kp2
260 do ist = sp_st1, sp_st2
261 do idim = 1, sp_dim
262 distdot = distdot + dmf_dotp_aux(x(j: j+sp_np-1), y(j:j+sp_np-1))
263 distdot = distdot + dmf_dotp_aux(x(k: k+sp_np-1), y(k:k+sp_np-1))
264 j = j + sp_np
265 k = k + sp_np
266 end do
267 end do
268 end do
270
271 case (3)
273 j = 1
274 k = n/2+1
275 do ik = sp_kp1, sp_kp2
276 do ist = sp_st1, sp_st2
277 do idim = 1, sp_dim
278 distdot = distdot + dmf_dotp_aux(x(j: j+sp_np-1), y(j:j+sp_np-1))
279 distdot = distdot + dmf_dotp_aux(x(k: k+sp_np-1), y(k:k+sp_np-1))
280 j = j + sp_np
281 k = k + sp_np
282 end do
283 end do
284 end do
285 do ik = sp_kp1, sp_kp2
286 do ist = sp_st1, sp_st2
287 do idim = 1, sp_dim
288 distdot = distdot + dmf_dotp_aux(x(j: j+sp_np-1), y(j:j+sp_np-1))
289 distdot = distdot + dmf_dotp_aux(x(k: k+sp_np-1), y(k:k+sp_np-1))
290 j = j + sp_np
291 k = k + sp_np
292 end do
293 end do
294 end do
296
297 end select
298
299end function distdot
300
301#endif /* HAVE_SPARSKIT */
302
303!! Local Variables:
304!! mode: f90
305!! coding: utf-8
306!! End:
scales a vector by a constant
Definition: lalg_basic.F90:157
double sin(double __x) __attribute__((__nothrow__
double atan(double __x) __attribute__((__nothrow__
double cos(double __x) __attribute__((__nothrow__
real(real64) function distdot(n, x, ix, y, iy)
This function will be linked by SPARSKIT to perform dot products. It expects complex numbers as an ar...
This module contains interfaces for BLAS routines You should not use these routines directly....
Definition: blas.F90:118
real(real64), parameter, public m_two
Definition: global.F90:190
real(real64), parameter, public m_zero
Definition: global.F90:188
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:186
real(real64), parameter, public m_half
Definition: global.F90:194
This module defines various routines, operating on mesh functions.
integer, public sp_st2
subroutine, public zmf_random(mesh, ff, pre_shift, post_shift, seed, normalized)
This subroutine fills a function with random values.
real(real64) function zmf_nrm2_2(mesh, dim, ff, reduce)
this function returns the norm of a vector of mesh functions
subroutine, public dmf_random(mesh, ff, pre_shift, post_shift, seed, normalized)
This subroutine fills a function with random values.
subroutine, public zmf_interpolate_points(ndim, npoints_in, x_in, f_in, npoints_out, x_out, f_out)
This function receives a function f_in defined in a mesh, and returns the interpolated values of the ...
complex(real64) function, public zmf_moment(mesh, ff, idir, order)
This function calculates the "order" moment of the function ff.
real(real64) function, public dmf_dotu_aux(f1, f2)
dot product between two vectors (mesh functions) without conjugation
real(real64) function dmf_line_integral_vector(mesh, ff, line)
This subroutine calculates the line integral of a vector function on a given line.
subroutine, public zmf_normalize(mesh, dim, psi, norm)
Normalize a mesh function psi.
complex(real64) function zmf_line_integral_vector(mesh, ff, line)
This subroutine calculates the line integral of a vector function on a given line.
real(real64) function, dimension(1:dim) dmf_integrate_2(mesh, dim, ff, mask, reduce)
Integrate of a vector of functins.
complex(real64) function zmf_dotp_2(mesh, dim, f1, f2, reduce, dotu, np)
dot product for vector valued mesh functions
real(real64) function dmf_integrate_1(mesh, ff, mask, reduce)
Integrate a function on the mesh.
class(mesh_t), pointer, public mesh_aux
Globally-scoped pointer to the mesh instance.
complex(real64) function, dimension(1:dim) zmf_integrate_2(mesh, dim, ff, mask, reduce)
Integrate of a vector of functins.
complex(real64) function, public zmf_dotu_aux(f1, f2)
dot product between two vectors (mesh functions) without conjugation
logical, public sp_parallel
integer, public sp_distdot_mode
real(real64) function dmf_nrm2_2(mesh, dim, ff, reduce)
this function returns the norm of a vector of mesh functions
real(real64) function dmf_surface_integral_scalar(mesh, ff, plane)
This subroutine calculates the surface integral of a scalar function on a given plane.
integer, public sp_dim
complex(real64) function zmf_line_integral_scalar(mesh, ff, line)
This subroutine calculates the line integral of a scalar function on a given line.
subroutine, public dmf_interpolate_points(ndim, npoints_in, x_in, f_in, npoints_out, x_out, f_out)
This function receives a function f_in defined in a mesh, and returns the interpolated values of the ...
real(real64) function dmf_dotp_2(mesh, dim, f1, f2, reduce, dotu, np)
dot product for vector valued mesh functions
integer, public sp_kp2
real(real64) function, public dmf_dotp_aux(f1, f2)
dot product between two vectors (mesh functions)
real(real64) function zmf_nrm2_1(mesh, ff, reduce)
this function returns the norm of a mesh function
real(real64) function dmf_line_integral_scalar(mesh, ff, line)
This subroutine calculates the line integral of a scalar function on a given line.
real(real64) function, public zmf_nrm2_aux(ff)
calculate norm2 of a vector (mesh function)
subroutine, public dmf_multipoles(mesh, ff, lmax, multipole, mask)
This routine calculates the multipoles of a function ff.
subroutine, public zmf_multipoles(mesh, ff, lmax, multipole, mask)
This routine calculates the multipoles of a function ff.
complex(real64) function zmf_surface_integral_scalar(mesh, ff, plane)
This subroutine calculates the surface integral of a scalar function on a given plane.
integer, public sp_kp1
subroutine, public dmf_normalize(mesh, dim, psi, norm)
Normalize a mesh function psi.
type(mpi_grp_t), public sp_grp
real(real64) function dmf_dotp_1(mesh, f1, f2, reduce, dotu, np)
this function returns the dot product between two mesh functions
real(real64) function dmf_surface_integral_vector(mesh, ff, plane)
This subroutine calculates the surface integral of a vector function on a given plane.
subroutine, public zmf_fix_phase(mesh, ff)
Fix the phase of complex function.
integer, public sp_st1
real(real64) function dmf_nrm2_1(mesh, ff, reduce)
this function returns the norm of a mesh function
complex(real64) function zmf_surface_integral_vector(mesh, ff, plane)
This subroutine calculates the surface integral of a vector function on a given plane.
subroutine, public dmf_dipole(mesh, ff, dipole, mask)
This routine calculates the dipole of a function ff, for arbitrary dimensions.
subroutine, public mesh_init_mesh_aux(mesh)
Initialise a pointer to the grid/mesh, that is globally exposed, such that low level mesh operations ...
real(real64) function, public dmf_moment(mesh, ff, idir, order)
This function calculates the "order" moment of the function ff.
real(real64) function, public dmf_nrm2_aux(ff)
calculate norm2 of a vector (mesh function)
integer, public sp_np
subroutine, public zmf_dipole(mesh, ff, dipole, mask)
This routine calculates the dipole of a function ff, for arbitrary dimensions.
complex(real64) function zmf_integrate_1(mesh, ff, mask, reduce)
Integrate a function on the mesh.
complex(real64) function zmf_dotp_1(mesh, f1, f2, reduce, dotu, np)
this function returns the dot product between two mesh functions
complex(real64) function, public zmf_dotp_aux(f1, f2)
dot product between two vectors (mesh functions)
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
Some general things and nomenclature:
Definition: par_vec.F90:171