ParaMonte Fortran 2.0.0
Parallel Monte Carlo and Machine Learning Library
See the latest version documentation.
pm_distanceEuclid::getDisEuclid Interface Reference

Generate and return the (squared) Euclidean distance of a (set of) point(s) in ndim-dimensions from a reference point (possibly origin), optionally robustly without underflow or overflow.
More...

Detailed Description

Generate and return the (squared) Euclidean distance of a (set of) point(s) in ndim-dimensions from a reference point (possibly origin), optionally robustly without underflow or overflow.

Parameters
[in]x: The input scalar (or array of the same rank as other array-like arguments) of the same type and kind as the output distance containing the x component of a 3D vector whose Euclidean norm is to be computed.
(optional. It must be present if and only if the input arguments point and ref are missing and y and z are present.)
[in]y: The input scalar (or array of the same rank as other array-like arguments), of the same type and kind as x, containing the y component of a 2D or 3D vector whose Euclidean norm is to be computed.
(optional. It must be present if and only if the input arguments point and ref are missing and x an z are present.)
[in]z: The input scalar (or array of the same rank as other array-like arguments), of the same type and kind as x, containing the z component of a 3D vector whose Euclidean norm is to be computed.
(optional. It must be present if and only if the input arguments point and ref are missing and x and y are present.)
[in]point: The input contiguous vector of shape (1:ndim) or matrix of shape (1:ndim, 1:npnt) of the same type and kind as the output distance, containing a (set of npnt) point(s) in the ndim-dimensional Euclidean space whose distances with respect to the input reference point ref must be returned.
(optional. It must be present if and only if the input arguments x, y, z are missing.)
[in]ref: The input contiguous vector of shape (1:ndim) or matrix of shape (1:ndim, 1:nref) of the same type and kind as point, containing the (set of nref) reference point(s) from which the distance(s) of point must be computed.
(optional, default = [(0., i = 1, size(point, 1))]. It can be present only if the input argument point is also present.)
[in]method: The input scalar that can be,
  1. the constant euclid or an object of type euclid_type, implying that all distance calculations must be done without undue numerical overflow.
    This option is computationally the most expensive method.
  2. the constant euclidu or an object of type euclidu_type, implying that all distance calculations must be without runtime checks for numerical overflow.
    This option is computationally faster than the euclid method.
  3. the constant euclidv or an object of type euclidv_type, implying that the volumes corresponding to all Euclidean distances must be returned without runtime checks for numerical overflow.
    This option is computationally the fastest approach to computing the distances because it avoid costly sqrt() operations and runtime overflow
  4. the constant euclidsq or an object of type euclidsq_type implying that the squared values of all distance calculations must be returned without runtime checks for numerical overflow.
    This option is computationally the fastest approach to computing the distances because it avoid costly sqrt() operations and runtime overflow checks.
(optional, default = euclid)
Returns
distance : The output object of,
  1. type real of kind any supported by the processor (e.g., RK, RK32, RK64, or RK128),
containing the requested Euclidean (squared) distance(s).
The rank and shape of the output distance follows that of the interfaces illustrated below.


Possible calling interfaces

! distance with respect to origin.
distance = getDisEuclid(x, y, z) ! elemental
distance = getDisEuclid(point(1:ndim), method = method)
distance(1:npnt) = getDisEuclid(point(1:ndim, 1:npnt), method = method)
! distance with respect to custom reference.
distance = getDisEuclid(point(1:ndim), ref(1:ndim), method = method)
distance(1:nref) = getDisEuclid(point(1:ndim), ref(1:ndim, 1:nref), method = method)
distance(1:npnt) = getDisEuclid(point(1:ndim, 1:npnt), ref(1:ndim), method = method)
distance(1:npnt, 1:nref) = getDisEuclid(point(1:ndim, 1:npnt), ref(1:ndim, 1:nref), method = method)
!
Generate and return the (squared) Euclidean distance of a (set of) point(s) in ndim-dimensions from a...
This module contains procedures and generic interfaces for computing the Euclidean norm of a single p...
Warning
The condition size(point, 1) == size(ref, 1) must hold for the corresponding input arguments.
This condition is verified only if the library is built with the preprocessor macro CHECK_ENABLED=1.
The pure procedure(s) documented herein become impure when the ParaMonte library is compiled with preprocessor macro CHECK_ENABLED=1.
By default, these procedures are pure in release build and impure in debug and testing builds.
Remarks
The procedures under discussion are elemental. The elemental attribute does not apply to interfaces that take point as an input argument.
Note
The Fortran standard provides the intrinsic procedure norm2() for computing the Euclidean norm of a vector.
However, the standard does not enforce robustness of the intrinsic procedure with respect to possible underflows or overflows.
The procedures of this module ensure robustness of the distance computations.
This will inevitably lead to worse runtime performance compared to the compiler implementations of the intrinsic routine that do not respect robustness.
Use the routines of this module in place of the Fortran intrinsics if you believe there is a possibility of under/over-flow.
The procedures of this module can be used for a robust computation of abs(x) when x is a large complex value.
In such a case, calling getDisEuclid([x%re, x%im]) would be equivalent to abs(x) intrinsic operation.
However, note that the Fortran standard already offers a better intrinsic alternative to the routines of this procedure for this task, namely hypot() which is robust against overflow and underflow.
This generic interface intentionally does not have explicit procedures for 2D Euclidean distance (x, y) because the Fortran intrinsic procedure hypot() already serves the purpose.
BLAS/LAPACK equivalent:
The procedures under discussion combine, modernize, and extend the interface and functionalities of Version 3.11 of BLAS/LAPACK routine(s): dlapy3
See also
getDisEuclid
setDisEuclid
getDisMatEuclid
setDisMatEuclid
Intrinsic Fortran procedure hypot(x, y) (robust)
Intrinsic Fortran procedure norm2(x(:)) (unsafe)
BLAS/LAPACK equivalent:
The procedures under discussion combine, modernize, and extend the interface and functionalities of Version 3.11 of BLAS/LAPACK routine(s): dlapy3


Example usage

1program example
2
3 use pm_kind, only: SK, IK, LK, RKH
4 use pm_kind, only: RKG => RKS ! all processor kinds are supported.
5 use pm_io, only: display_type
7
8 implicit none
9
10 type(display_type) :: disp
11 complex(RKG) :: z
12 real(RKG) , allocatable :: point(:)
13 real(RKG) , allocatable :: ref(:)
14 disp = display_type(file = "main.out.F90")
15
16 call disp%skip()
17 call disp%show("point = [real(RKG) :: 1, 3, 2]")
18 point = [real(RKG) :: 1, 3, 2]
19 call disp%show("point")
20 call disp%show( point )
21 call disp%show("[getDisEuclid(point), getDisEuclid(point(1), point(2), point(3))]") ! norm2(point),
22 call disp%show( [getDisEuclid(point), getDisEuclid(point(1), point(2), point(3))] ) ! norm2(point),
23 call disp%skip()
24
25 call disp%skip()
26 call disp%show("point = [real(RKG) :: 1, 1, 1] * huge(0._RKG) / 10")
27 point = [real(RKG) :: 1, 1, 1] * huge(0._RKG) / 10
28 call disp%show("point")
29 call disp%show( point )
30 call disp%show("[getDisEuclid(point), getDisEuclid(point(1), point(2), point(3))] !, sqrt(dot_product(point, point)) norm2(point), ! Note that GNU gfortran `norm2()` respects robustness, while Intel ifort does not.")
31 call disp%show( [getDisEuclid(point), getDisEuclid(point(1), point(2), point(3))] ) !, sqrt(dot_product(point, point))
32 call disp%skip()
33
34 call disp%skip()
35 call disp%show("z = (1._RKG, -1._RKG) * huge(0._RKG) / 10")
36 z = (1._RKG, -1._RKG) * huge(0._RKG) / 10
37 call disp%show("point")
38 call disp%show( point )
39 call disp%show("[hypot(z%re, z%im), abs(z), getDisEuclid([z%re, z%im])] ! norm2([z%re, z%im]), Note that GNU gfortran `norm2()` respects robustness, while Intel ifort does not.")
40 call disp%show( [hypot(z%re, z%im), abs(z), getDisEuclid([z%re, z%im])] )
41 call disp%skip()
42
43 call disp%skip()
44 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
45 call disp%show("! Compute the distance of a point with respect to a reference in arbitrary dimensions without undue overflow/underflow.")
46 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
47 call disp%skip()
48
49 call disp%skip()
50 call disp%show("point = [real(RKG) :: 1, 0, 5, 2]")
51 point = [real(RKG) :: 1, 0, 5, 2]
52 call disp%show("ref = [real(RKG) :: 2, 1, 3, 0]")
53 ref = [real(RKG) :: 2, 1, 3, 0]
54 call disp%show("point")
55 call disp%show( point )
56 call disp%show("[getDisEuclid(point - ref), getDisEuclid(point, ref)]") ! norm2(point - ref),
57 call disp%show( [getDisEuclid(point - ref), getDisEuclid(point, ref)] ) ! norm2(point - ref),
58 call disp%skip()
59
60 call disp%skip()
61 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
62 call disp%show("! Compute the asymmetric matrix of (squared) distances of a set of points from a set of reference points with or without undue overflow/underflow.")
63 call disp%show("!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
64 call disp%skip()
65
66 block
67 real(RKG), allocatable :: point(:,:), ref(:,:), distance(:,:)
68 call disp%skip()
69 call disp%show("point = reshape([real(RKG) :: 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1], shape = [3, 5])")
70 point = reshape([real(RKG) :: 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1], shape = [3, 5])
71 call disp%show("point")
72 call disp%show( point )
73 call disp%show("ref = reshape([real(RKG) :: -1, -1, -1, 1, 1, 1], shape = [3, 2]) ! two reference points.")
74 ref = reshape([real(RKG) :: -1, -1, -1, 1, 1, 1], shape = [3, 2])
75 call disp%show("ref")
76 call disp%show( ref )
77 call disp%show("distance = getDisEuclid(point, ref)")
78 distance = getDisEuclid(point, ref)
79 call disp%show("distance")
80 call disp%show( distance )
81 call disp%show("distance = getDisEuclid(point, ref, euclid)")
82 distance = getDisEuclid(point, ref, euclid)
83 call disp%show("distance")
84 call disp%show( distance )
85 call disp%show("distance = getDisEuclid(point, ref, euclidu)")
86 distance = getDisEuclid(point, ref, euclidu)
87 call disp%show("distance")
88 call disp%show( distance )
89 call disp%show("distance = getDisEuclid(point, ref, euclidsq)")
90 distance = getDisEuclid(point, ref, euclidsq)
91 call disp%show("distance")
92 call disp%show( distance )
93 call disp%skip()
94 end block
95
96end program example
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11726
This is a generic method of the derived type display_type with pass attribute.
Definition: pm_io.F90:11508
type(euclidu_type), parameter euclidu
This is a scalar parameter object of type euclidu_typethat is exclusively used to request unsafe meth...
type(euclid_type), parameter euclid
This is a scalar parameter object of type euclid_type that is exclusively used to request safe method...
type(euclidsq_type), parameter euclidsq
This is a scalar parameter object of type euclidsq_typethat is exclusively used to request computing ...
This module contains classes and procedures for input/output (IO) or generic display operations on st...
Definition: pm_io.F90:252
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
Definition: pm_io.F90:11393
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
Definition: pm_kind.F90:268
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
Definition: pm_kind.F90:541
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
Definition: pm_kind.F90:540
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
Definition: pm_kind.F90:539
integer, parameter RKH
The scalar integer constant of intrinsic default kind, representing the highest-precision real kind t...
Definition: pm_kind.F90:858
integer, parameter RKS
The single-precision real kind in Fortran mode. On most platforms, this is an 32-bit real kind.
Definition: pm_kind.F90:567
Generate and return an object of type display_type.
Definition: pm_io.F90:10282

Example Unix compile command via Intel ifort compiler
1#!/usr/bin/env sh
2rm main.exe
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example Windows Batch compile command via Intel ifort compiler
1del main.exe
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
4main.exe

Example Unix / MinGW compile command via GNU gfortran compiler
1#!/usr/bin/env sh
2rm main.exe
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
4./main.exe

Example output
1
2point = [real(RKG) :: 1, 3, 2]
3point
4+1.00000000, +3.00000000, +2.00000000
5[getDisEuclid(point), getDisEuclid(point(1), point(2), point(3))]
6+3.74165726, +3.74165726
7
8
9point = [real(RKG) :: 1, 1, 1] * huge(0._RKG) / 10
10point
11+0.340282347E+38, +0.340282347E+38, +0.340282347E+38
12[getDisEuclid(point), getDisEuclid(point(1), point(2), point(3))] !, sqrt(dot_product(point, point)) norm2(point), ! Note that GNU gfortran `norm2()` respects robustness, while Intel ifort does not.
13+0.589386287E+38, +0.589386287E+38
14
15
16z = (1._RKG, -1._RKG) * huge(0._RKG) / 10
17point
18+0.340282347E+38, +0.340282347E+38, +0.340282347E+38
19[hypot(z%re, z%im), abs(z), getDisEuclid([z%re, z%im])] ! norm2([z%re, z%im]), Note that GNU gfortran `norm2()` respects robustness, while Intel ifort does not.
20+0.481231910E+38, +0.481231910E+38, +0.481231910E+38
21
22
23!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24! Compute the distance of a point with respect to a reference in arbitrary dimensions without undue overflow/underflow.
25!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26
27
28point = [real(RKG) :: 1, 0, 5, 2]
29ref = [real(RKG) :: 2, 1, 3, 0]
30point
31+1.00000000, +0.00000000, +5.00000000, +2.00000000
32[getDisEuclid(point - ref), getDisEuclid(point, ref)]
33+3.16227770, +3.16227770
34
35
36!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37! Compute the asymmetric matrix of (squared) distances of a set of points from a set of reference points with or without undue overflow/underflow.
38!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39
40
41point = reshape([real(RKG) :: 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1], shape = [3, 5])
42point
43+0.00000000, +1.00000000, +0.00000000, +0.00000000, +1.00000000
44+0.00000000, +0.00000000, +1.00000000, +0.00000000, +1.00000000
45+0.00000000, +0.00000000, +0.00000000, +1.00000000, +1.00000000
46ref = reshape([real(RKG) :: -1, -1, -1, 1, 1, 1], shape = [3, 2]) ! two reference points.
47ref
48-1.00000000, +1.00000000
49-1.00000000, +1.00000000
50-1.00000000, +1.00000000
51distance = getDisEuclid(point, ref)
52distance
53+1.73205078, +1.73205078
54+2.44948983, +1.41421354
55+2.44948983, +1.41421354
56+2.44948983, +1.41421354
57+3.46410155, +0.00000000
58distance = getDisEuclid(point, ref, euclid)
59distance
60+1.73205078, +1.73205078
61+2.44948983, +1.41421354
62+2.44948983, +1.41421354
63+2.44948983, +1.41421354
64+3.46410155, +0.00000000
65distance = getDisEuclid(point, ref, euclidu)
66distance
67+1.73205078, +1.73205078
68+2.44948983, +1.41421354
69+2.44948983, +1.41421354
70+2.44948983, +1.41421354
71+3.46410155, +0.00000000
72distance = getDisEuclid(point, ref, euclidsq)
73distance
74+3.00000000, +3.00000000
75+6.00000000, +2.00000000
76+6.00000000, +2.00000000
77+6.00000000, +2.00000000
78+12.0000000, +0.00000000
79
80
Test:
test_pm_distanceEuclid
Todo:
Normal Priority: A benchmark comparison with the equivalent compiler implementations would be informative.


Final Remarks


If you believe this algorithm or its documentation can be improved, we appreciate your contribution and help to edit this page's documentation and source file on GitHub.
For details on the naming abbreviations, see this page.
For details on the naming conventions, see this page.
This software is distributed under the MIT license with additional terms outlined below.

  1. If you use any parts or concepts from this library to any extent, please acknowledge the usage by citing the relevant publications of the ParaMonte library.
  2. If you regenerate any parts/ideas from this library in a programming environment other than those currently supported by this ParaMonte library (i.e., other than C, C++, Fortran, MATLAB, Python, R), please also ask the end users to cite this original ParaMonte library.

This software is available to the public under a highly permissive license.
Help us justify its continued development and maintenance by acknowledging its benefit to society, distributing it, and contributing to it.

Author:
Amir Shahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas Austin

Definition at line 484 of file pm_distanceEuclid.F90.


The documentation for this interface was generated from the following file: