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,
-
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.
-
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.
-
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
-
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,
-
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(
1:npnt)
= getDisEuclid(point(
1:ndim,
1:npnt), method
= method)
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.
- 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 ⛓
10 type(display_type) :: disp
12 real(RKG) ,
allocatable :: point(:)
13 real(RKG) ,
allocatable :: ref(:)
17 call disp%show(
"point = [real(RKG) :: 1, 3, 2]")
18 point
= [
real(RKG) ::
1,
3,
2]
21 call disp%show(
"[getDisEuclid(point), getDisEuclid(point(1), point(2), point(3))]")
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
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.")
35 call disp%show(
"z = (1._RKG, -1._RKG) * huge(0._RKG) / 10")
36 z
= (
1._RKG,
-1._RKG)
* huge(
0._RKG)
/ 10
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.")
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(
"!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
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]
56 call disp%show(
"[getDisEuclid(point - ref), getDisEuclid(point, ref)]")
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(
"!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%")
67 real(RKG),
allocatable :: point(:,:), ref(:,:), distance(:,:)
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])
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])
77 call disp%show(
"distance = getDisEuclid(point, ref)")
81 call disp%show(
"distance = getDisEuclid(point, ref, euclid)")
85 call disp%show(
"distance = getDisEuclid(point, ref, euclidu)")
89 call disp%show(
"distance = getDisEuclid(point, ref, euclidsq)")
This is a generic method of the derived type display_type with pass attribute.
This is a generic method of the derived type display_type with pass attribute.
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...
type(display_type) disp
This is a scalar module variable an object of type display_type for general display.
This module defines the relevant Fortran kind type-parameters frequently used in the ParaMonte librar...
integer, parameter LK
The default logical kind in the ParaMonte library: kind(.true.) in Fortran, kind(....
integer, parameter IK
The default integer kind in the ParaMonte library: int32 in Fortran, c_int32_t in C-Fortran Interoper...
integer, parameter SK
The default character kind in the ParaMonte library: kind("a") in Fortran, c_char in C-Fortran Intero...
integer, parameter RKH
The scalar integer constant of intrinsic default kind, representing the highest-precision real kind t...
integer, parameter RKS
The single-precision real kind in Fortran mode. On most platforms, this is an 32-bit real kind.
Generate and return an object of type display_type.
Example Unix compile command via Intel ifort compiler ⛓
3ifort -fpp -standard-semantics -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
Example Windows Batch compile command via Intel ifort compiler ⛓
2set PATH=..\..\..\lib;%PATH%
3ifort /fpp /standard-semantics /O3 /I:..\..\..\include main.F90 ..\..\..\lib\libparamonte*.lib /exe:main.exe
Example Unix / MinGW compile command via GNU gfortran compiler ⛓
3gfortran -cpp -ffree-line-length-none -O3 -Wl,-rpath,../../../lib -I../../../inc main.F90 ../../../lib/libparamonte* -o main.exe
Example output ⛓
2point
= [
real(RKG) ::
1,
3,
2]
4+1.00000000,
+3.00000000,
+2.00000000
6+3.74165726,
+3.74165726
9point
= [
real(RKG) ::
1,
1,
1]
* huge(
0._RKG)
/ 10
11+0.340282347E+38,
+0.340282347E+38,
+0.340282347E+38
13+0.589386287E+38,
+0.589386287E+38
16z
= (
1._RKG,
-1._RKG)
* huge(
0._RKG)
/ 10
18+0.340282347E+38,
+0.340282347E+38,
+0.340282347E+38
20+0.481231910E+38,
+0.481231910E+38,
+0.481231910E+38
28point
= [
real(RKG) ::
1,
0,
5,
2]
29ref
= [
real(RKG) ::
2,
1,
3,
0]
31+1.00000000,
+0.00000000,
+5.00000000,
+2.00000000
33+3.16227770,
+3.16227770
41point
= reshape([
real(RKG) ::
0,
0,
0,
1,
0,
0,
0,
1,
0,
0,
0,
1,
1,
1,
1], shape
= [
3,
5])
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])
48-1.00000000,
+1.00000000
49-1.00000000,
+1.00000000
50-1.00000000,
+1.00000000
53+1.73205078,
+2.44948983,
+2.44948983,
+2.44948983,
+3.46410155
54+1.73205078,
+1.41421354,
+1.41421354,
+1.41421354,
+0.00000000
57+1.73205078,
+2.44948983,
+2.44948983,
+2.44948983,
+3.46410155
58+1.73205078,
+1.41421354,
+1.41421354,
+1.41421354,
+0.00000000
61+1.73205078,
+2.44948983,
+2.44948983,
+2.44948983,
+3.46410155
62+1.73205078,
+1.41421354,
+1.41421354,
+1.41421354,
+0.00000000
65+3.00000000,
+6.00000000,
+6.00000000,
+6.00000000,
+12.0000000
66+3.00000000,
+2.00000000,
+2.00000000,
+2.00000000,
+0.00000000
- 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.
-
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.
-
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.
- Copyright
- Computational Data Science Lab
- Author:
- Amir Shahmoradi, September 1, 2017, 12:00 AM, Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin
Definition at line 484 of file pm_distanceEuclid.F90.
1.9.3