diff --git a/README.md b/README.md index 3bfe1a5..ec28fed 100644 --- a/README.md +++ b/README.md @@ -7,10 +7,34 @@ Delaunay triangulation. In addition to the original Fortran source code, this repository contains a wrapper for Python 3.6+ and C bindings. Command line drivers are also provided with the original Fortran code. -## Organization and Usage +## Usage -The physical organization is as follows. Note that each of the following -directories could be independently downloaded. +`DELAUNAYSPARSE` contains several modes of operation. + +In the original ACM TOMS release, two Fortran driver subroutines were provided: + * `DELAUNAYSPARSES` runs the serial driver to identify the vertices + of the simplex/simplices containing one or more interpolation points. + Can also (optionally) be set to compute and return the value of the + Delaunay interpolant. + * `DELAUNAYSPARSEP` runs the parallel driver to identify the vertices + of the simplex/simplices containing one or more interpolation points. + Can also (optionally) be set to compute and return the value of the + Delaunay interpolant (must set the `OMP_NUM_THREADS` environment + variable). + +Additionally, two command-line drivers are provided, which read input +from files: + * `delsparses` (uses the serial driver), and + * `delsparsep` (uses the parallel driver). + +In this repository, two additional interfaces are provided for calling +from C/C++ (`c_binding`) and Python 3 (`python`). + +Further detailed user information is documented in the `USAGE` document. + +## Organization + +The physical organization is as follows. * `src` contains the original unmodified Fortran source code, as published in ACM TOMS Algorithm 1012. This includes 2 command line drivers @@ -30,6 +54,10 @@ directories could be independently downloaded. A test file `test_install.c` can be used for usage examples. This directory's internal README also contains best practices when calling Fortran from C/C++. + * `docs` contains the html source for generating the DelaunaySparse website. + * `USAGE` provides additional detailed user information. + * DelaunaySparse is shared under the MIT Software License, in the `LICENSE` + file. ## Citation diff --git a/USAGE.md b/USAGE.md new file mode 100644 index 0000000..d09f3c5 --- /dev/null +++ b/USAGE.md @@ -0,0 +1,549 @@ +# Usage Information for DELAUNAYSPARSE. + +DELAUNAYSPARSE solves the multivariate interpolation problem: + +Given a set of `N` points `PTS` and a set of `M` interpolation points +`Q` in `R^D`, for each interpolation point `Q_i` in `Q`, identify the set +of `D+1` data points in `PTS` that are the vertices of a Delaunay simplex +containing `Q_i`. + +These vertices can be used to calculate the Delaunay interpolant using +a piecewise linear model. + +For more information on the underlying algorithm, see + +Chang et al. 2018. A polynomial time algorithm for multivariate interpolation +in arbitrary dimension via the Delaunay triangulation. In Proc. ACMSE 2018 +Conference. + +For more information on this software, see + +Chang et al. 2020. Algorithm 1012: DELAUNAYSPARSE: Interpolation via a sparse +subset of the Delaunay triangulation in medium to high dimensions. +ACM Trans. Math. Softw. 46(4). Article No. 38. + +DELAUNAYSPARSE contains a Fortran module + * `delsparse`; + +as well as C bindings + * `delsparsec`; + +two command-line drivers + * `delsparses` and + * `delsparsep`; + +and a Python 3 wrapper + * `delsparsepy`. + +These interfaces are described in the following sections. + +## Fortran interface + +DELAUNAYSPARSE is written in Fortran 2003, and this is its native interface. +The Fortran interface contains two drivers: + * `DELAUNAYSPARSES` (serial driver) and + * `DELAUNAYSPARSEP` (OpenMP parallel driver). + +### DELAUNAYSPARSES + +The interface for DELAUNAYSPARSES is + +``` +SUBROUTINE DELAUNAYSPARSES( D, N, PTS, M, Q, SIMPS, WEIGHTS, IERR, & + INTERP_IN, INTERP_OUT, EPS, EXTRAP, RNORM, & + IBUDGET, CHAIN, EXACT ) +``` + +Each of the above parameters is described below. + + +On input: + + * `D` is the dimension of the space for `PTS` and `Q`. + + * `N` is the number of data points in `PTS`. + + * `PTS(1:D,1:N)` is a real valued matrix with `N` columns, each containing the + coordinates of a single data point in R^D. + + * `M` is the number of interpolation points in `Q`. + + * `Q(1:D,1:M)` is a real valued matrix with `M` columns, each containing the + coordinates of a single interpolation point in R^D. + + +On output: + + * `PTS` and `Q` have been rescaled and shifted. All the data points in `PTS` + are now contained in the unit hyperball in R^D, and the points in `Q` + have been shifted and scaled accordingly in relation to `PTS`. + + * `SIMPS(1:D+1,1:M)` contains the `D+1` integer indices (corresponding to + columns in `PTS`) for the `D+1` vertices of the Delaunay simplex + containing each interpolation point in `Q`. + + * `WEIGHTS(1:D+1,1:M)` contains the `D+1` real-valued weights for expressing + each point in `Q` as a convex combination of the `D+1` corresponding vertices + in `SIMPS`. + + * `IERR(1:M)` contains integer valued error flags associated with the + computation of each of the `M` interpolation points in `Q`. The error + codes are: + + Codes 0, 1, 2 are expected to occur during normal execution. + + - 00 : Succesful interpolation. + - 01 : Succesful extrapolation (up to the allowed extrapolation distance). + - 02 : This point was outside the allowed extrapolation distance; the + corresponding entries in SIMPS and WEIGHTS contain zero values. + + Error codes 10--28 indicate that one or more inputs contain illegal + values or are incompatible with each other. + + - 10 : The dimension `D` must be positive. + - 11 : Too few data points to construct a triangulation (i.e., `N < D+1`). + - 12 : No interpolation points given (i.e., `M < 1`). + - 13 : The first dimension of `PTS` does not agree with the dimension `D`. + - 14 : The second dimension of `PTS` does not agree with the number of + points `N`. + - 15 : The first dimension of `Q` does not agree with the dimension `D`. + - 16 : The second dimension of `Q` does not agree with the number of + interpolation points `M`. + - 17 : The first dimension of the output array `SIMPS` does not match the + number of vertices needed for a `D`-simplex (`D+1`). + - 18 : The second dimension of the output array `SIMPS` does not match the + number of interpolation points `M`. + - 19 : The first dimension of the output array `WEIGHTS` does not match the + number of vertices for a a `D`-simplex (`D+1`). + - 20 : The second dimension of the output array `WEIGHTS` does not match + the number of interpolation points `M`. + - 21 : The size of the error array `IERR` does not match the number of + interpolation points `M`. + - 22 : `INTERP_IN` cannot be present without `INTERP_OUT` or vice versa. + - 23 : The first dimension of `INTERP_IN` does not match the first + dimension of `INTERP_OUT`. + - 24 : The second dimension of `INTERP_IN` does not match the number of + data points `PTS`. + - 25 : The second dimension of `INTERP_OUT` does not match the number of + interpolation points `M`. + - 26 : The budget supplied in `IBUDGET` does not contain a positive + integer. + - 27 : The extrapolation distance supplied in `EXTRAP` cannot be negative. + - 28 : The size of the `RNORM` output array does not match the number of + interpolation points `M`. + + The errors 30, 31 typically indicate that DELAUNAYSPARSE has been given + an unclean dataset. These errors can be fixed by preprocessing your + data (remove duplicate points and apply PCA or other dimension reduction + technique). + + - 30 : Two or more points in the data set `PTS` are too close together with + respect to the working precision (`EPS`), which would result in a + numerically degenerate simplex. + - 31 : All the data points in `PTS` lie in some lower dimensional linear + manifold (up to the working precision), and no valid triangulation + exists. + + The error code 40 occurs when another earlier error prevented this point + from ever being evaluated. + + - 40 : An error caused `DELAUNAYSPARSES` to terminate before this value + could be computed. Note: The corresponding entries in `SIMPS` and + `WEIGHTS` may contain garbage values. + + The error code 50 corresponds to allocation of the internal WORK array. + Check your systems internal memory settings and limits, in relation + to the problem size and DELAUNAYSPARSE's space requirements (see TOMS + Alg. paper for more details on DELAUNAYSPARSE's space requirements). + + - 50 : A memory allocation error occurred while allocating the work array + `WORK`. + + The errors 60, 61 should not occur with the default settings. If one of + these errors is observed, then it is likely that either the value of + the optional inputs `IBUDGET` or `EPS` has been adjusted in a way that is + unwise, or there may be another issue with the problem settings, which + is manifesting in an unusual way. + + - 60 : The budget was exceeded before the algorithm converged on this + value. If the dimension is high, try increasing `IBUDGET`. This + error can also be caused by a working precision `EPS` that is too + small for the conditioning of the problem. + - 61 : A value that was judged appropriate later caused LAPACK to + encounter a singularity. Try increasing the value of `EPS`. + + The errors 70--72 were caused by the DWNNLS library from SLATEC, which + is only used during extrapolation. Note that there is a known issue + with this library, when it is linked against included public-domain + copy of BLAS/LAPACK, instead of an installed version + (i.e., `-lblas` `-llapack`). + + - 70 : Allocation error for the extrapolation work arrays. + - 71 : The SLATEC subroutine `DWNNLS` failed to converge during the + projection of an extrapolation point onto the convex hull. + - 72 : The SLATEC subroutine `DWNNLS` has reported a usage error. + + The errors 72, 80--83 should never occur, and likely indicate a + compiler bug or hardware failure. + + - 80 : The LAPACK subroutine `DGEQP3` has reported an illegal value. + - 81 : The LAPACK subroutine `DGETRF` has reported an illegal value. + - 82 : The LAPACK subroutine `DGETRS` has reported an illegal value. + - 83 : The LAPACK subroutine `DORMQR` has reported an illegal value. + + +Optional arguments: + + * `INTERP_IN(1:IR,1:N)` contains real valued response vectors for each of + the data points in `PTS` on input. The first dimension of `INTERP_IN` is + inferred to be the dimension of these response vectors, and the + second dimension must match `N`. If present, the response values will + be computed for each interpolation point in `Q`, and stored in `INTERP_OUT`, + which therefore must also be present. If both `INTERP_IN` and `INTERP_OUT` + are omitted, only the containing simplices and convex combination + weights are returned. + + * `INTERP_OUT(1:IR,1:M)` contains real valued response vectors for each + interpolation point in `Q` on output. The first dimension of `INTERP_OUT` + must match the first dimension of `INTERP_IN`, and the second dimension + must match `M`. If present, the response values at each interpolation + point are computed as a convex combination of the response values + (supplied in `INTERP_IN`) at the vertices of a Delaunay simplex containing + that interpolation point. Therefore, if `INTERP_OUT` is present, then + `INTERP_IN` must also be present. If both are omitted, only the + simplices and convex combination weights are returned. + + * `EPS` contains the real working precision for the problem on input. By + default, `EPS` is assigned \sqrt{\mu} where \mu denotes the unit roundoff + for the machine. In general, any values that differ by less than `EPS` + are judged as equal, and any weights that are greater than `-EPS` are + judged as nonnegative. `EPS` cannot take a value less than the default + value of \sqrt{\mu}. If any value less than \sqrt{\mu} is supplied, + the default value will be used instead automatically. + + * `EXTRAP` contains the real maximum extrapolation distance (relative to the + diameter of `PTS`) on input. Interpolation at a point outside the convex + hull of `PTS` is done by projecting that point onto the convex hull, and + then doing normal Delaunay interpolation at that projection. + Interpolation at any point in `Q` that is more than `EXTRAP * DIAMETER(PTS)` + units outside the convex hull of `PTS` will not be done and an error code + of `2` will be returned. Note that computing the projection can be + expensive. Setting `EXTRAP=0` will cause all extrapolation points to be + ignored without ever computing a projection. By default, `EXTRAP=0.1` + (extrapolate by up to 10% of the diameter of `PTS`). + + * `RNORM(1:M)` contains the real unscaled projection (2-norm) distances from + any projection computations on output. If not present, these distances + are still computed for each extrapolation point, but are never returned. + + * `IBUDGET` on input contains the integer budget for performing flips while + iterating toward the simplex containing each interpolation point in + `Q`. This prevents `DELAUNAYSPARSES` from falling into an infinite loop when + an inappropriate value of `EPS` is given with respect to the problem + conditioning. By default, `IBUDGET=50000`. However, for extremely + high-dimensional problems and pathological inputs, the default value + may be insufficient. + + * `CHAIN` is a logical input argument that determines whether a new first + simplex should be constructed for each interpolation point + (`CHAIN=.FALSE.`), or whether the simplex walks should be "daisy-chained." + By default, `CHAIN=.FALSE.` Setting `CHAIN=.TRUE.` is generally not + recommended, unless the size of the triangulation is relatively small + or the interpolation points are known to be tightly clustered. + + * `EXACT` is a logical input argument that determines whether the exact + diameter should be computed and whether a check for duplicate data + points should be performed in advance. When `EXACT=.FALSE.`, the + diameter of `PTS` is approximated by twice the distance from the + barycenter of `PTS` to the farthest point in `PTS`, and no check is + done to find the closest pair of points, which could result in hard + to find bugs later on. When `EXACT=.TRUE.`, the exact diameter is + computed and an error is returned whenever PTS contains duplicate + values up to the precision `EPS`. By default `EXACT=.TRUE.`, but setting + `EXACT=.FALSE.` could result in significant speedup when `N` is large. + It is strongly recommended that most users leave `EXACT=.TRUE.`, as + setting `EXACT=.FALSE.` could result in input errors that are difficult + to identify. Also, the diameter approximation could be wrong by up to + a factor of two. + + +Subroutines and functions directly referenced from BLAS are + * `DDOT`, + * `DGEMV`, + * `DNRM2`, + * `DTRSM`, +and from LAPACK are + * `DGEQP3`, + * `DGETRF`, + * `DGETRS`, + * `DORMQR`. + +The SLATEC subroutine + * `DWNNLS` is also directly referenced. + +`DWNNLS` and all its SLATEC dependencies have been slightly edited to +comply with the Fortran 2008 standard, with all print statements and +references to stderr being commented out. For a reference to `DWNNLS`, +see ACM TOMS Algorithm 587 (Hanson and Haskell). +The module `REAL_PRECISION` from HOMPACK90 (ACM TOMS Algorithm 777) is +used for the real data type. The `REAL_PRECISION` module, `DELAUNAYSPARSES`, +and `DWNNLS` and its dependencies comply with the Fortran 2008 standard. + +## DELAUNAYSPARSEP + +``` +SUBROUTINE DELAUNAYSPARSEP( D, N, PTS, M, Q, SIMPS, WEIGHTS, IERR, & + INTERP_IN, INTERP_OUT, EPS, EXTRAP, RNORM, IBUDGET, CHAIN, EXACT, & + PMODE ) +``` + +Each of the above parameters is described below. + + +On input: + + * `D` is the dimension of the space for `PTS` and `Q`. + + * `N` is the number of data points in `PTS`. + + * `PTS(1:D,1:N)` is a real valued matrix with `N` columns, each containing the + coordinates of a single data point in R^D. + + * `M` is the number of interpolation points in `Q`. + + * `Q(1:D,1:M)` is a real valued matrix with `M` columns, each containing the + coordinates of a single interpolation point in R^D. + + +On output: + + * `PTS` and `Q` have been rescaled and shifted. All the data points in `PTS` + are now contained in the unit hyperball in R^D, and the points in `Q` + have been shifted and scaled accordingly in relation to `PTS`. + + * `SIMPS(1:D+1,1:M)` contains the `D+1` integer indices (corresponding to + columns in `PTS`) for the `D+1` vertices of the Delaunay simplex + containing each interpolation point in `Q`. + + * `WEIGHTS(1:D+1,1:M)` contains the `D+1` real-valued weights for expressing + each point in `Q` as a convex combination of the `D+1` corresponding vertices + in `SIMPS`. + + * `IERR(1:M)` contains integer valued error flags associated with the + computation of each of the `M` interpolation points in `Q`. The error + codes are: + + Codes 0, 1, 2 are expected to occur during normal execution. + + - 00 : Succesful interpolation. + - 01 : Succesful extrapolation (up to the allowed extrapolation distance). + - 02 : This point was outside the allowed extrapolation distance; the + corresponding entries in SIMPS and WEIGHTS contain zero values. + + Error codes 10--28 indicate that one or more inputs contain illegal + values or are incompatible with each other. + + - 10 : The dimension `D` must be positive. + - 11 : Too few data points to construct a triangulation (i.e., `N < D+1`). + - 12 : No interpolation points given (i.e., `M < 1`). + - 13 : The first dimension of `PTS` does not agree with the dimension `D`. + - 14 : The second dimension of `PTS` does not agree with the number of + points `N`. + - 15 : The first dimension of `Q` does not agree with the dimension `D`. + - 16 : The second dimension of `Q` does not agree with the number of + interpolation points `M`. + - 17 : The first dimension of the output array `SIMPS` does not match the + number of vertices needed for a `D`-simplex (`D+1`). + - 18 : The second dimension of the output array `SIMPS` does not match the + number of interpolation points `M`. + - 19 : The first dimension of the output array `WEIGHTS` does not match the + number of vertices for a a `D`-simplex (`D+1`). + - 20 : The second dimension of the output array `WEIGHTS` does not match + the number of interpolation points `M`. + - 21 : The size of the error array `IERR` does not match the number of + interpolation points `M`. + - 22 : `INTERP_IN` cannot be present without `INTERP_OUT` or vice versa. + - 23 : The first dimension of `INTERP_IN` does not match the first + dimension of `INTERP_OUT`. + - 24 : The second dimension of `INTERP_IN` does not match the number of + data points `PTS`. + - 25 : The second dimension of `INTERP_OUT` does not match the number of + interpolation points `M`. + - 26 : The budget supplied in `IBUDGET` does not contain a positive + integer. + - 27 : The extrapolation distance supplied in `EXTRAP` cannot be negative. + - 28 : The size of the `RNORM` output array does not match the number of + interpolation points `M`. + + The errors 30, 31 typically indicate that DELAUNAYSPARSE has been given + an unclean dataset. These errors can be fixed by preprocessing your + data (remove duplicate points and apply PCA or other dimension reduction + technique). + + - 30 : Two or more points in the data set `PTS` are too close together with + respect to the working precision (`EPS`), which would result in a + numerically degenerate simplex. + - 31 : All the data points in `PTS` lie in some lower dimensional linear + manifold (up to the working precision), and no valid triangulation + exists. + + The error code 40 occurs when another earlier error prevented this point + from ever being evaluated. + + - 40 : An error caused `DELAUNAYSPARSEP` to terminate before this value + could be computed. Note: The corresponding entries in `SIMPS` and + `WEIGHTS` may contain garbage values. + + The error code 50 corresponds to allocation of the internal WORK array. + Check your systems internal memory settings and limits, in relation + to the problem size and DELAUNAYSPARSE's space requirements (see TOMS + Alg. paper for more details on DELAUNAYSPARSE's space requirements). + + - 50 : A memory allocation error occurred while allocating the work array + `WORK`. + + The errors 60, 61 should not occur with the default settings. If one of + these errors is observed, then it is likely that either the value of + the optional inputs `IBUDGET` or `EPS` has been adjusted in a way that is + unwise, or there may be another issue with the problem settings, which + is manifesting in an unusual way. + + - 60 : The budget was exceeded before the algorithm converged on this + value. If the dimension is high, try increasing `IBUDGET`. This + error can also be caused by a working precision `EPS` that is too + small for the conditioning of the problem. + - 61 : A value that was judged appropriate later caused LAPACK to + encounter a singularity. Try increasing the value of `EPS`. + + The errors 70--72 were caused by the DWNNLS library from SLATEC, which + is only used during extrapolation. Note that there is a known issue + with this library, when it is linked against included public-domain + copy of BLAS/LAPACK, instead of an installed version + (i.e., `-lblas` `-llapack`). + + - 70 : Allocation error for the extrapolation work arrays. + - 71 : The SLATEC subroutine `DWNNLS` failed to converge during the + projection of an extrapolation point onto the convex hull. + - 72 : The SLATEC subroutine `DWNNLS` has reported a usage error. + + The errors 72, 80--83 should never occur, and likely indicate a + compiler bug or hardware failure. + + - 80 : The LAPACK subroutine `DGEQP3` has reported an illegal value. + - 81 : The LAPACK subroutine `DGETRF` has reported an illegal value. + - 82 : The LAPACK subroutine `DGETRS` has reported an illegal value. + - 83 : The LAPACK subroutine `DORMQR` has reported an illegal value. + + The error code 90 is unique to DELAUNAYSPARSEP. + + - 90 : The value of `PMODE` is not valid. + + +Optional arguments: + + * `INTERP_IN(1:IR,1:N)` contains real valued response vectors for each of + the data points in `PTS` on input. The first dimension of `INTERP_IN` is + inferred to be the dimension of these response vectors, and the + second dimension must match `N`. If present, the response values will + be computed for each interpolation point in `Q`, and stored in `INTERP_OUT`, + which therefore must also be present. If both `INTERP_IN` and `INTERP_OUT` + are omitted, only the containing simplices and convex combination + weights are returned. + + * `INTERP_OUT(1:IR,1:M)` contains real valued response vectors for each + interpolation point in `Q` on output. The first dimension of `INTERP_OUT` + must match the first dimension of `INTERP_IN`, and the second dimension + must match `M`. If present, the response values at each interpolation + point are computed as a convex combination of the response values + (supplied in `INTERP_IN`) at the vertices of a Delaunay simplex containing + that interpolation point. Therefore, if `INTERP_OUT` is present, then + `INTERP_IN` must also be present. If both are omitted, only the + simplices and convex combination weights are returned. + + * `EPS` contains the real working precision for the problem on input. By + default, `EPS` is assigned \sqrt{\mu} where \mu denotes the unit roundoff + for the machine. In general, any values that differ by less than `EPS` + are judged as equal, and any weights that are greater than `-EPS` are + judged as nonnegative. `EPS` cannot take a value less than the default + value of \sqrt{\mu}. If any value less than \sqrt{\mu} is supplied, + the default value will be used instead automatically. + + * `EXTRAP` contains the real maximum extrapolation distance (relative to the + diameter of `PTS`) on input. Interpolation at a point outside the convex + hull of `PTS` is done by projecting that point onto the convex hull, and + then doing normal Delaunay interpolation at that projection. + Interpolation at any point in `Q` that is more than `EXTRAP * DIAMETER(PTS)` + units outside the convex hull of `PTS` will not be done and an error code + of `2` will be returned. Note that computing the projection can be + expensive. Setting `EXTRAP=0` will cause all extrapolation points to be + ignored without ever computing a projection. By default, `EXTRAP=0.1` + (extrapolate by up to 10% of the diameter of `PTS`). + + * `RNORM(1:M)` contains the real unscaled projection (2-norm) distances from + any projection computations on output. If not present, these distances + are still computed for each extrapolation point, but are never returned. + + * `IBUDGET` on input contains the integer budget for performing flips while + iterating toward the simplex containing each interpolation point in + `Q`. This prevents `DELAUNAYSPARSEP` from falling into an infinite loop when + an inappropriate value of `EPS` is given with respect to the problem + conditioning. By default, `IBUDGET=50000`. However, for extremely + high-dimensional problems and pathological inputs, the default value + may be insufficient. + + * `CHAIN` is a logical input argument that determines whether a new first + simplex should be constructed for each interpolation point + (`CHAIN=.FALSE.`), or whether the simplex walks should be "daisy-chained." + By default, `CHAIN=.FALSE.` Setting `CHAIN=.TRUE.` is generally not + recommended, unless the size of the triangulation is relatively small + or the interpolation points are known to be tightly clustered. + + * `EXACT` is a logical input argument that determines whether the exact + diameter should be computed and whether a check for duplicate data + points should be performed in advance. When `EXACT=.FALSE.`, the + diameter of `PTS` is approximated by twice the distance from the + barycenter of `PTS` to the farthest point in `PTS`, and no check is + done to find the closest pair of points, which could result in hard + to find bugs later on. When `EXACT=.TRUE.`, the exact diameter is + computed and an error is returned whenever PTS contains duplicate + values up to the precision `EPS`. By default `EXACT=.TRUE.`, but setting + `EXACT=.FALSE.` could result in significant speedup when `N` is large. + It is strongly recommended that most users leave `EXACT=.TRUE.`, as + setting `EXACT=.FALSE.` could result in input errors that are difficult + to identify. Also, the diameter approximation could be wrong by up to + a factor of two. + + * `PMODE` is an integer specifying the level of parallelism to be exploited. + - If `PMODE = 1`, then parallelism is exploited at the level of the loop + over all interpolation points (Level 1 parallelism). + - If `PMODE = 2`, then parallelism is exploited at the level of the loops + over data points when constructing/flipping simplices (Level 2 + parallelism). + - If `PMODE = 3`, then parallelism is exploited at both levels. Note: + this implies that the total number of threads active at any time could + be up to `OMP_NUM_THREADS^2`. + By default, `PMODE` is set to `1` if there is more than 1 interpolation + point and `2` otherwise. + + +Subroutines and functions directly referenced from BLAS are + * `DDOT`, + * `DGEMV`, + * `DNRM2`, + * `DTRSM`, +and from LAPACK are + * `DGEQP3`, + * `DGETRF`, + * `DGETRS`, + * `DORMQR`. + +The SLATEC subroutine + * `DWNNLS` is also directly referenced. + +`DWNNLS` and all its SLATEC dependencies have been slightly edited to +comply with the Fortran 2008 standard, with all print statements and +references to stderr being commented out. For a reference to `DWNNLS`, +see ACM TOMS Algorithm 587 (Hanson and Haskell). +The module `REAL_PRECISION` from HOMPACK90 (ACM TOMS Algorithm 777) is +used for the real data type. The `REAL_PRECISION` module, `DELAUNAYSPARSEP`, +and `DWNNLS` and its dependencies comply with the Fortran 2008 standard. diff --git a/c_binding/Makefile b/c_binding/Makefile index bc04468..592d3ef 100644 --- a/c_binding/Makefile +++ b/c_binding/Makefile @@ -2,14 +2,15 @@ FORT = gfortran CC = gcc CFLAGS = -c OPTS = -fopenmp +LIBS = blas.f lapack.f LEGACY = -std=legacy -all: test_install.o delsparse_bind_c.o delsparse.o slatec.o lapack.o blas.o - $(FORT) $(OPTS) test_install.o delsparse_bind_c.o delsparse.o slatec.o lapack.o blas.o -o test_install +all: test_install.o delsparse_bind_c.o delsparse.o slatec.o delsparse.h + $(FORT) $(OPTS) test_install.o delsparse_bind_c.o delsparse.o slatec.o $(LIBS) -o test_install ./test_install test_install.o: test_install.c - $(CC) $(CFLAGS) $(OPTS) test_install.c -o test_install.o + $(CC) $(CFLAGS) test_install.c -o test_install.o delsparse_bind_c.o: delsparse_bind_c.f90 delsparse.o $(FORT) $(CFLAGS) $(OPTS) delsparse_bind_c.f90 -o delsparse_bind_c.o @@ -20,11 +21,5 @@ delsparse.o: delsparse.f90 slatec.o : slatec.f $(FORT) $(CFLAGS) $(OPTS) $(LEGACY) slatec.f -o slatec.o -lapack.o : lapack.f - $(FORT) $(CFLAGS) $(OPTS) lapack.f -o lapack.o - -blas.o : blas.f - $(FORT) $(CFLAGS) $(OPTS) blas.f -o blas.o - clean: rm -f *.o *.mod test_install diff --git a/c_binding/delsparse.h b/c_binding/delsparse.h new file mode 100644 index 0000000..4ed0241 --- /dev/null +++ b/c_binding/delsparse.h @@ -0,0 +1,59 @@ +#ifndef DELSPARSEC +#define DELSPARSEC + +// serial subroutine: no optional arguments +extern void c_delaunaysparses(int *d, int *n, double pts[], int *m, double q[], + int simps[], double weights[], int ierr[]); + +// serial: compute interpolant values +extern void c_delaunaysparses_interp(int *d, int *n, double pts[], int *m, + double q[], int simps[], double weights[], + int ierr[], int *ir, double interp_in[], + double interp_out[]); + +// serial: optional arguments, no interpolant values +extern void c_delaunaysparses_opts(int *d, int *n, double pts[], int *m, + double q[],int simps[], double weights[], + int ierr[], double *eps, double *extrap, + double rnorm[], int *ibudget, bool *chain, + bool *exact); + +// serial: optional arguments and compute interpolant values +extern void c_delaunaysparses_interp_opts(int *d, int *n, double pts[], int *m, + double q[],int simps[], + double weights[], int ierr[], + int *ir, double interp_in[], + double interp_out[], double *eps, + double *extrap, double rnorm[], + int *ibudget, bool *chain, + bool *exact); + + +// parallel: no optional arguments +extern void c_delaunaysparsep(int *d, int *n, double pts[], int *m, double q[], + int simps[], double weights[], int ierr[]); + +// parallel: compute interpolant values +extern void c_delaunaysparsep_interp(int *d, int *n, double pts[], int *m, + double q[], int simps[], double weights[], + int ierr[], int *ir, double interp_in[], + double interp_out[]); + +// parallel: optional arguments, no interpolant values +extern void c_delaunaysparsep_opts(int *d, int *n, double pts[], int *m, + double q[],int simps[], double weights[], + int ierr[], double *eps, double *extrap, + double rnorm[], int *ibudget, bool *chain, + bool *exact, int *pmode); + +// parallel: optional arguments and compute interpolant values +extern void c_delaunaysparsep_interp_opts(int *d, int *n, double pts[], int *m, + double q[],int simps[], + double weights[], int ierr[], + int *ir, double interp_in[], + double interp_out[], double *eps, + double *extrap, double rnorm[], + int *ibudget, bool *chain, + bool *exact, int *pmode); + +#endif diff --git a/c_binding/slatec.f b/c_binding/slatec.f index 7d51578..c652a26 100644 --- a/c_binding/slatec.f +++ b/c_binding/slatec.f @@ -7,7 +7,7 @@ SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, C a covariance matrix. C***LIBRARY SLATEC C***CATEGORY K1A2A, D9 -C***TYPE REAL(KIND=R8) (LSEI-S, DLSEI-D) +C***TYPE DOUBLE PRECISION (LSEI-S, DLSEI-D) C***KEYWORDS CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING, C EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS, C QUADRATIC PROGRAMMING @@ -62,7 +62,7 @@ SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, C C The parameters for DLSEI( ) are C -C Input.. All TYPE REAL variables are REAL(KIND=R8) +C Input.. All TYPE REAL variables are DOUBLE PRECISION C C W(*,*),MDW, The array W(*,*) is doubly subscripted with C ME,MA,MG,N first dimensioning parameter equal to MDW. @@ -268,7 +268,7 @@ SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, C LIP = MG+2*N+2 C This test will not be made if IP(2).LE.0. C -C Output.. All TYPE REAL variables are REAL(KIND=R8) +C Output.. All TYPE REAL variables are DOUBLE PRECISION C C X(*),RNORME, The array X(*) contains the solution parameters C RNORML if the integer output flag MODE = 0 or 1. @@ -382,18 +382,17 @@ SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, C 900510 Convert XERRWV calls to XERMSG calls. (RWC) C 900604 DP version created from SP version. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) -C 180613 Removed prints and replaced DP --> REAL(KIND=R8). (THC) +C 180613 Removed prints and replaced DP --> DOUBLE PRECISION. (THC) C***END PROLOGUE DLSEI - USE REAL_PRECISION INTEGER IP(3), MA, MDW, ME, MG, MODE, N - REAL(KIND=R8) PRGOPT(*), RNORME, RNORML, W(MDW,*), WS(*), X(*) + DOUBLE PRECISION PRGOPT(*), RNORME, RNORML, W(MDW,*), WS(*), X(*) C EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DDOT, DH12, DLSI, DNRM2, * DSCAL, DSWAP - REAL(KIND=R8) D1MACH, DASUM, DDOT, DNRM2 + DOUBLE PRECISION D1MACH, DASUM, DDOT, DNRM2 C - REAL(KIND=R8) DRELPR, ENORM, FNORM, GAM, RB, RN, RNMAX, SIZE, + DOUBLE PRECISION DRELPR, ENORM, FNORM, GAM, RB, RN, RNMAX, SIZE, * SN, SNMAX, T, TAU, UJ, UP, VJ, XNORM, XNRME INTEGER I, IMAX, J, JP1, K, KEY, KRANKE, LAST, LCHK, LINK, M, * MAPKE1, MDEQC, MEND, MEP1, N1, N2, NEXT, NLINK, NOPT, NP1, @@ -743,7 +742,7 @@ SUBROUTINE DLSI (W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, C***SUBSIDIARY C***PURPOSE Subsidiary to DLSEI C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (LSI-S, DLSI-D) +C***TYPE DOUBLE PRECISION (LSI-S, DLSI-D) C***AUTHOR Hanson, R. J., (SNLA) C***DESCRIPTION C @@ -795,16 +794,15 @@ SUBROUTINE DLSI (W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, C 900604 DP version created from SP version. (RWC) C 920422 Changed CALL to DHFTI to include variable MA. (WRB) C***END PROLOGUE DLSI - USE REAL_PRECISION INTEGER IP(*), MA, MDW, MG, MODE, N - REAL(KIND=R8) PRGOPT(*), RNORM, W(MDW,*), WS(*), X(*) + DOUBLE PRECISION PRGOPT(*), RNORM, W(MDW,*), WS(*), X(*) C EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DDOT, DH12, DHFTI, DLPDP, * DSCAL, DSWAP - REAL(KIND=R8) D1MACH, DASUM, DDOT + DOUBLE PRECISION D1MACH, DASUM, DDOT C - REAL(KIND=R8) ANORM, DRELPR, FAC, GAM, RB, TAU, TOL, XNORM, + DOUBLE PRECISION ANORM, DRELPR, FAC, GAM, RB, TAU, TOL, XNORM, * TMP_NORM(1) INTEGER I, J, K, KEY, KRANK, KRM1, KRP1, L, LAST, LINK, M, MAP1, * MDLPDP, MINMAN, N1, N2, N3, NEXT, NP1 @@ -1079,12 +1077,12 @@ SUBROUTINE DLSI (W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, RETURN END *DECK D1MACH - REAL(KIND=R8) FUNCTION D1MACH (I) + DOUBLE PRECISION FUNCTION D1MACH (I) C***BEGIN PROLOGUE D1MACH C***PURPOSE Return floating point machine dependent constants. C***LIBRARY SLATEC C***CATEGORY R1 -C***TYPE REAL(KIND=R8) (R1MACH-S, D1MACH-D) +C***TYPE DOUBLE PRECISION (R1MACH-S, D1MACH-D) C***KEYWORDS MACHINE CONSTANTS C***AUTHOR Fox, P. A., (Bell Labs) C Hall, A. D., (Bell Labs) @@ -1151,7 +1149,6 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C comments below. (DWL) C***END PROLOGUE D1MACH C - USE REAL_PRECISION INTEGER SMALL(4) INTEGER LARGE(4) @@ -1164,7 +1161,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C for DMACH(2) is a slight lower bound. The value for DMACH(5) is C a 20-digit approximation. If one of the sets of initial data below C is preferred, do the necessary commenting and uncommenting. (DWL) - REAL(KIND=R8) DMACH(5) + DOUBLE PRECISION DMACH(5) DATA DMACH / 2.23D-308, 1.79D+308, 1.111D-16, 2.222D-16, 1 0.30102999566398119521D0 / SAVE DMACH @@ -1387,7 +1384,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C MACHINE CONSTANTS FOR THE ELXSI 6400 -C (ASSUMING REAL*8 IS THE DEFAULT REAL(KIND=R8)) +C (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION) C C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / @@ -1420,7 +1417,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE HP 2100 -C THREE WORD REAL(KIND=R8) OPTION WITH FTN4 +C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / @@ -1429,7 +1426,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 -C FOUR WORD REAL(KIND=R8) OPTION WITH FTN4 +C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / @@ -1461,7 +1458,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE IBM PC -C ASSUMES THAT ALL ARITHMETIC IS DONE IN REAL(KIND=R8) +C ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION C ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087. C C DATA SMALL(1) / 2.23D-308 / @@ -2176,7 +2173,7 @@ INTEGER FUNCTION I1MACH (I) C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE HP 2100 -C 3 WORD REAL(KIND=R8) OPTION WITH FTN4 +C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / @@ -2196,7 +2193,7 @@ INTEGER FUNCTION I1MACH (I) C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 -C 4 WORD REAL(KIND=R8) OPTION WITH FTN4 +C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / @@ -2507,11 +2504,11 @@ SUBROUTINE DH12 (MODE, LPIVOT, L1, M, U, IUE, UP, C, ICE, ICV, C***SUBSIDIARY C***PURPOSE Subsidiary to DHFTI, DLSEI and DWNNLS C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (H12-S, DH12-D) +C***TYPE DOUBLE PRECISION (H12-S, DH12-D) C***AUTHOR (UNKNOWN) C***DESCRIPTION C -C *** REAL(KIND=R8) VERSION OF H12 ****** +C *** DOUBLE PRECISION VERSION OF H12 ****** C C C.L.Lawson and R.J.Hanson, Jet Propulsion Laboratory, 1973 Jun 12 C to appear in 'Solving Least Squares Problems', Prentice-Hall, 1974 @@ -2548,13 +2545,12 @@ SUBROUTINE DH12 (MODE, LPIVOT, L1, M, U, IUE, UP, C, ICE, ICV, C 890831 Modified array declarations. (WRB) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900328 Added TYPE section. (WRB) -C 900911 Added DDOT to REAL(KIND=R8) statement. (WRB) +C 900911 Added DDOT to DOUBLE PRECISION statement. (WRB) C***END PROLOGUE DH12 - USE REAL_PRECISION INTEGER I, I2, I3, I4, ICE, ICV, INCR, IUE, J, KL1, KL2, KLP, * L1, L1M1, LPIVOT, M, MML1P2, MODE, NCV - REAL(KIND=R8) B, C, CL, CLINV, ONE, UL1M1, SM, U, UP, DDOT + DOUBLE PRECISION B, C, CL, CLINV, ONE, UL1M1, SM, U, UP, DDOT DIMENSION U(IUE,*), C(*) C BEGIN BLOCK PERMITTING ...EXITS TO 140 C***FIRST EXECUTABLE STATEMENT DH12 @@ -2654,7 +2650,7 @@ SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, C Exactly one right-hand side vector is permitted. C***LIBRARY SLATEC C***CATEGORY D9 -C***TYPE REAL(KIND=R8) (HFTI-S, DHFTI-D) +C***TYPE DOUBLE PRECISION (HFTI-S, DHFTI-D) C***KEYWORDS CURVE FITTING, LEAST SQUARES C***AUTHOR Lawson, C. L., (JPL) C Hanson, R. J., (SNLA) @@ -2711,7 +2707,7 @@ SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, C C The entire set of parameters for DHFTI are C -C INPUT.. All TYPE REAL variables are REAL(KIND=R8) +C INPUT.. All TYPE REAL variables are DOUBLE PRECISION C C A(*,*),MDA,M,N The array A(*,*) initially contains the M by N C matrix A of the least squares problem AX = B. @@ -2742,7 +2738,7 @@ SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, C C H(*),G(*),IP(*) Arrays of working space used by DHFTI. C -C OUTPUT.. All TYPE REAL variables are REAL(KIND=R8) +C OUTPUT.. All TYPE REAL variables are DOUBLE PRECISION C C A(*,*) The contents of the array A(*,*) will be C modified by the subroutine. These contents @@ -2782,11 +2778,10 @@ SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, C 901005 Replace usage of DDIFF with usage of D1MACH. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DHFTI - USE REAL_PRECISION INTEGER I, II, IOPT, IP(*), IP1, J, JB, JJ, K, KP1, KRANK, L, * LDIAG, LMAX, M, MDA, MDB, N, NB, NERR - REAL(KIND=R8) A, B, D1MACH, DZERO, FACTOR, + DOUBLE PRECISION A, B, D1MACH, DZERO, FACTOR, * G, H, HMAX, RELEPS, RNORM, SM, SM1, SZERO, TAU, TMP DIMENSION A(MDA,*),B(MDB,*),H(*),G(*),RNORM(*) SAVE RELEPS @@ -2985,7 +2980,7 @@ SUBROUTINE DLPDP (A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS, C***SUBSIDIARY C***PURPOSE Subsidiary to DLSEI C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (LPDP-S, DLPDP-D) +C***TYPE DOUBLE PRECISION (LPDP-S, DLPDP-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -3028,12 +3023,11 @@ SUBROUTINE DLPDP (A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS, C 900328 Added TYPE section. (WRB) C 910408 Updated the AUTHOR section. (WRB) C***END PROLOGUE DLPDP - USE REAL_PRECISION C INTEGER I, IS(*), IW, IX, J, L, M, MDA, MODE, MODEW, N, N1, N2, * NP1 - REAL(KIND=R8) A(MDA,*), DDOT, DNRM2, FAC, ONE, + DOUBLE PRECISION A(MDA,*), DDOT, DNRM2, FAC, ONE, * PRGOPT(*), RNORM, SC, WNORM, WS(*), X(*), YNORM, ZERO SAVE ZERO, ONE, FAC DATA ZERO,ONE /0.0D0,1.0D0/, FAC /0.1D0/ @@ -3197,7 +3191,7 @@ SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, C selected variables. C***LIBRARY SLATEC C***CATEGORY K1A2A -C***TYPE REAL(KIND=R8) (WNNLS-S, DWNNLS-D) +C***TYPE DOUBLE PRECISION (WNNLS-S, DWNNLS-D) C***KEYWORDS CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING, C EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS, C NONNEGATIVITY CONSTRAINTS, QUADRATIC PROGRAMMING @@ -3232,7 +3226,7 @@ SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, C C The parameters for DWNNLS are C -C INPUT.. All TYPE REAL variables are REAL(KIND=R8) +C INPUT.. All TYPE REAL variables are DOUBLE PRECISION C C W(*,*),MDW, The array W(*,*) is double subscripted with first C ME,MA,N,L dimensioning parameter equal to MDW. For this @@ -3392,7 +3386,7 @@ SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, C LIW = ME+MA+N C This test will not be made if IWORK(2).LE.0. C -C OUTPUT.. All TYPE REAL variables are REAL(KIND=R8) +C OUTPUT.. All TYPE REAL variables are DOUBLE PRECISION C C X(*) An array dimensioned at least N, which will C contain the N components of the solution vector @@ -3455,13 +3449,12 @@ SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, C 900510 Convert XERRWV calls to XERMSG calls, change Prologue C comments to agree with WNNLS. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) -C 180613 Removed prints and replaced DP --> REAL(KIND=R8). (THC) +C 180613 Removed prints and replaced DP --> DOUBLE PRECISION. (THC) C***END PROLOGUE DWNNLS - USE REAL_PRECISION INTEGER IWORK(*), L, L1, L2, L3, L4, L5, LIW, LW, MA, MDW, ME, * MODE, N - REAL(KIND=R8) PRGOPT(*), RNORM, W(MDW,*), WORK(*), X(*) + DOUBLE PRECISION PRGOPT(*), RNORM, W(MDW,*), WORK(*), X(*) C CHARACTER*8 XERN1 C***FIRST EXECUTABLE STATEMENT DWNNLS MODE = 0 @@ -3525,7 +3518,7 @@ SUBROUTINE DWNLSM (W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE, C***SUBSIDIARY C***PURPOSE Subsidiary to DWNNLS C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLSM-S, DWNLSM-D) +C***TYPE DOUBLE PRECISION (WNLSM-S, DWNLSM-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -3539,7 +3532,7 @@ SUBROUTINE DWNLSM (W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE, C sequence from DWNNLS for purposes of variable dimensioning). C Their contents will in general be of no interest to the user. C -C Variables of type REAL are REAL(KIND=R8). +C Variables of type REAL are DOUBLE PRECISION. C C IPIVOT(*) C An array of length N. Upon completion it contains the @@ -3590,18 +3583,17 @@ SUBROUTINE DWNLSM (W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE, C 900604 DP version created from SP version. (RWC) C 900911 Restriction on value of ALAMDA included. (WRB) C***END PROLOGUE DWNLSM - USE REAL_PRECISION INTEGER IPIVOT(*), ITYPE(*), L, MA, MDW, MME, MODE, N - REAL(KIND=R8) D(*), H(*), PRGOPT(*), RNORM, SCALE(*), TEMP(*), + DOUBLE PRECISION D(*), H(*), PRGOPT(*), RNORM, SCALE(*), TEMP(*), * W(MDW,*), WD(*), X(*), Z(*) C EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DH12, DNRM2, SLATEC_DROTM, * SLATEC_DROTMG, DSCAL, DSWAP, DWNLIT, IDAMAX, XERMSG - REAL(KIND=R8) D1MACH, DASUM, DNRM2 + DOUBLE PRECISION D1MACH, DASUM, DNRM2 INTEGER IDAMAX C - REAL(KIND=R8) ALAMDA, ALPHA, ALSQ, AMAX, BLOWUP, BNORM, + DOUBLE PRECISION ALAMDA, ALPHA, ALSQ, AMAX, BLOWUP, BNORM, * DOPE(3), DRELPR, EANORM, FAC, SM, SPARAM(5), T, TAU, WMAX, Z2, * ZZ INTEGER I, IDOPE(3), IMAX, ISOL, ITEMP, ITER, ITMAX, IWMAX, J, @@ -4177,7 +4169,7 @@ SUBROUTINE SLATEC_DROTM (N, DX, INCX, DY, INCY, DPARAM) C***PURPOSE Apply a modified Givens transformation. C***LIBRARY SLATEC (BLAS) C***CATEGORY D1A8 -C***TYPE REAL(KIND=R8) (SROTM-S, DROTM-D) +C***TYPE DOUBLE PRECISION (SROTM-S, DROTM-D) C***KEYWORDS BLAS, LINEAR ALGEBRA, MODIFIED GIVENS ROTATION, VECTOR C***AUTHOR Lawson, C. L., (JPL) C Hanson, R. J., (SNLA) @@ -4231,9 +4223,8 @@ SUBROUTINE SLATEC_DROTM (N, DX, INCX, DY, INCY, DPARAM) C 920501 Reformatted the REFERENCES section. (WRB) C 180613 Renamed SLATEC_DROTM to avoid BLAS naming conflict. (THC) C***END PROLOGUE SLATEC_DROTM - USE REAL_PRECISION - REAL(KIND=R8) DFLAG, DH12, DH22, DX, TWO, Z, DH11, DH21, + DOUBLE PRECISION DFLAG, DH12, DH22, DX, TWO, Z, DH11, DH21, 1 DPARAM, DY, W, ZERO DIMENSION DX(*), DY(*), DPARAM(5) SAVE ZERO, TWO @@ -4346,7 +4337,7 @@ SUBROUTINE SLATEC_DROTMG (DD1, DD2, DX1, DY1, DPARAM) C***PURPOSE Construct a modified Givens transformation. C***LIBRARY SLATEC (BLAS) C***CATEGORY D1B10 -C***TYPE REAL(KIND=R8) (SROTMG-S, DROTMG-D) +C***TYPE DOUBLE PRECISION (SROTMG-S, DROTMG-D) C***KEYWORDS BLAS, LINEAR ALGEBRA, MODIFIED GIVENS ROTATION, VECTOR C***AUTHOR Lawson, C. L., (JPL) C Hanson, R. J., (SNLA) @@ -4400,9 +4391,8 @@ SUBROUTINE SLATEC_DROTMG (DD1, DD2, DX1, DY1, DPARAM) C 920501 Reformatted the REFERENCES section. (WRB) C 180613 Renamed SLATEC_DROTMG to avoid BLAS naming conflict. (THC) C***END PROLOGUE SLATEC_DROTMG - USE REAL_PRECISION - REAL(KIND=R8) GAM, ONE, RGAMSQ, DD1, DD2, DH11, DH12, DH21, + DOUBLE PRECISION GAM, ONE, RGAMSQ, DD1, DD2, DH11, DH12, DH21, 1 DH22, DPARAM, DP1, DP2, DQ1, DQ2, DU, DY1, ZERO, 2 GAMSQ, DFLAG, DTEMP, DX1, TWO DIMENSION DPARAM(5) @@ -4579,7 +4569,7 @@ SUBROUTINE DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, C***SUBSIDIARY C***PURPOSE Subsidiary to DWNNLS C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLIT-S, DWNLIT-D) +C***TYPE DOUBLE PRECISION (WNLIT-S, DWNLIT-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -4605,10 +4595,9 @@ SUBROUTINE DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, C 900328 Added TYPE section. (WRB) C 900604 DP version created from SP version. . (RWC) C***END PROLOGUE DWNLIT - USE REAL_PRECISION INTEGER IDOPE(*), IPIVOT(*), ITYPE(*), L, M, MDW, N - REAL(KIND=R8) DOPE(*), H(*), RNORM, SCALE(*), W(MDW,*) + DOUBLE PRECISION DOPE(*), H(*), RNORM, SCALE(*), W(MDW,*) LOGICAL DONE C EXTERNAL DCOPY, DH12, SLATEC_DROTM, SLATEC_DROTMG, DSCAL, DSWAP, @@ -4616,7 +4605,7 @@ SUBROUTINE DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, INTEGER IDAMAX LOGICAL DWNLT2 C - REAL(KIND=R8) ALSQ, AMAX, EANORM, FACTOR, HBAR, RN, SPARAM(5), + DOUBLE PRECISION ALSQ, AMAX, EANORM, FACTOR, HBAR, RN, SPARAM(5), * T, TAU INTEGER I, I1, IMAX, IR, J, J1, JJ, JP, KRANK, L1, LB, LEND, ME, * MEND, NIV, NSOLN @@ -4869,7 +4858,7 @@ SUBROUTINE DWNLT1 (I, LEND, MEND, IR, MDW, RECALC, IMAX, HBAR, H, C***SUBSIDIARY C***PURPOSE Subsidiary to WNLIT C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLT1-S, DWNLT1-D) +C***TYPE DOUBLE PRECISION (WNLT1-S, DWNLT1-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -4885,10 +4874,9 @@ SUBROUTINE DWNLT1 (I, LEND, MEND, IR, MDW, RECALC, IMAX, HBAR, H, C 890620 Code extracted from WNLIT and made a subroutine. (RWC)) C 900604 DP version created from SP version. (RWC) C***END PROLOGUE DWNLT1 - USE REAL_PRECISION INTEGER I, IMAX, IR, LEND, MDW, MEND - REAL(KIND=R8) H(*), HBAR, SCALE(*), W(MDW,*) + DOUBLE PRECISION H(*), HBAR, SCALE(*), W(MDW,*) LOGICAL RECALC C EXTERNAL IDAMAX @@ -4934,7 +4922,7 @@ LOGICAL FUNCTION DWNLT2 (ME, MEND, IR, FACTOR, TAU, SCALE, WIC) C***SUBSIDIARY C***PURPOSE Subsidiary to WNLIT C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLT2-S, DWNLT2-D) +C***TYPE DOUBLE PRECISION (WNLT2-S, DWNLT2-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -4964,12 +4952,11 @@ LOGICAL FUNCTION DWNLT2 (ME, MEND, IR, FACTOR, TAU, SCALE, WIC) C 890620 Code extracted from WNLIT and made a subroutine. (RWC)) C 900604 DP version created from SP version. (RWC) C***END PROLOGUE DWNLT2 - USE REAL_PRECISION - REAL(KIND=R8) FACTOR, SCALE(*), TAU, WIC(*) + DOUBLE PRECISION FACTOR, SCALE(*), TAU, WIC(*) INTEGER IR, ME, MEND C - REAL(KIND=R8) RN, SN, T + DOUBLE PRECISION RN, SN, T INTEGER J C C***FIRST EXECUTABLE STATEMENT DWNLT2 @@ -4995,7 +4982,7 @@ SUBROUTINE DWNLT3 (I, IMAX, M, MDW, IPIVOT, H, W) C***SUBSIDIARY C***PURPOSE Subsidiary to WNLIT C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLT3-S, DWNLT3-D) +C***TYPE DOUBLE PRECISION (WNLT3-S, DWNLT3-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -5011,14 +4998,13 @@ SUBROUTINE DWNLT3 (I, IMAX, M, MDW, IPIVOT, H, W) C 890620 Code extracted from WNLIT and made a subroutine. (RWC)) C 900604 DP version created from SP version. (RWC) C***END PROLOGUE DWNLT3 - USE REAL_PRECISION INTEGER I, IMAX, IPIVOT(*), M, MDW - REAL(KIND=R8) H(*), W(MDW,*) + DOUBLE PRECISION H(*), W(MDW,*) C EXTERNAL DSWAP C - REAL(KIND=R8) T + DOUBLE PRECISION T INTEGER ITEMP C C***FIRST EXECUTABLE STATEMENT DWNLT3 diff --git a/c_binding/test_install.c b/c_binding/test_install.c index f34b97c..2bcaa1f 100644 --- a/c_binding/test_install.c +++ b/c_binding/test_install.c @@ -1,61 +1,7 @@ #include #include #include - -// serial subroutine: no optional arguments -extern void c_delaunaysparses(int *d, int *n, double pts[], int *m, double q[], - int simps[], double weights[], int ierr[]); - -// serial: compute interpolant values -extern void c_delaunaysparses_interp(int *d, int *n, double pts[], int *m, - double q[], int simps[], double weights[], - int ierr[], int *ir, double interp_in[], - double interp_out[]); - -// serial: optional arguments, no interpolant values -extern void c_delaunaysparses_opts(int *d, int *n, double pts[], int *m, - double q[],int simps[], double weights[], - int ierr[], double *eps, double *extrap, - double rnorm[], int *ibudget, bool *chain, - bool *exact); - -// serial: optional arguments and compute interpolant values -extern void c_delaunaysparses_interp_opts(int *d, int *n, double pts[], int *m, - double q[],int simps[], - double weights[], int ierr[], - int *ir, double interp_in[], - double interp_out[], double *eps, - double *extrap, double rnorm[], - int *ibudget, bool *chain, - bool *exact); - - -// parallel: no optional arguments -extern void c_delaunaysparsep(int *d, int *n, double pts[], int *m, double q[], - int simps[], double weights[], int ierr[]); - -// parallel: compute interpolant values -extern void c_delaunaysparsep_interp(int *d, int *n, double pts[], int *m, - double q[], int simps[], double weights[], - int ierr[], int *ir, double interp_in[], - double interp_out[]); - -// parallel: optional arguments, no interpolant values -extern void c_delaunaysparsep_opts(int *d, int *n, double pts[], int *m, - double q[],int simps[], double weights[], - int ierr[], double *eps, double *extrap, - double rnorm[], int *ibudget, bool *chain, - bool *exact, int *pmode); - -// parallel: optional arguments and compute interpolant values -extern void c_delaunaysparsep_interp_opts(int *d, int *n, double pts[], int *m, - double q[],int simps[], - double weights[], int ierr[], - int *ir, double interp_in[], - double interp_out[], double *eps, - double *extrap, double rnorm[], - int *ibudget, bool *chain, - bool *exact, int *pmode); +#include "delsparse.h" int main() { // Set the problem dimensions diff --git a/docs/LICENSE b/docs/LICENSE new file mode 100644 index 0000000..00ce8f0 --- /dev/null +++ b/docs/LICENSE @@ -0,0 +1,22 @@ +MIT License + +Copyright (c) 2020 Tyler H. Chang, Layne T. Watson, Thomas C. H. Lux, +Ali R. Butt, Kirk W. Cameron, and Yili Hong. + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/docs/body.css b/docs/body.css new file mode 100644 index 0000000..99df0e4 --- /dev/null +++ b/docs/body.css @@ -0,0 +1,249 @@ +body { +color:#404040; +font:76% Verdana,Tahoma,Arial,sans-serif; +line-height:1.2em; +margin:0 auto; +padding:0; +} + +a { +color:#387C44; +font-weight:700; +text-decoration:none; +} + +a:hover { +text-decoration:underline; +} + +a img { +border:0; +} + +p { +margin:0 0 18px 10px; +} + +blockquote { +border:1px solid #dadada; +font-size:0.9em; +margin:20px 10px; +padding:8px; +} + +h1 { +color:#387C44; +font-size:4.2em; +letter-spacing:-5px; +margin:0 0 30px 25px; +} + +h1 a { +color:#387C44; +text-transform:none; +} + +h2 { +border-bottom:4px solid #dadada; +color:#387C44; +font-size:1.4em; +letter-spacing:-1px; +margin:0 0 10px; +padding:0 2px 2px 5px; +} + +h3 { +border-bottom:1px solid #dadada; +color:#387C44; +font-size:1.2em; +font-weight:700; +margin:10px 0 8px; +padding:1px 2px 2px 3px; +} + +/*** Main wrap and header ***/ + +#wrap { +color:#404040; +margin:10px auto; +padding:0; +width:970px; +} + +#header { +margin:0; +} + +#toplinks { +font-size:0.9em; +padding:5px 2px 2px 3px; +text-align:right; +} + +#toplinks a { +color:gray; +} + +#slogan { +color:gray; +font-size:1.5em; +font-weight:700; +letter-spacing:-1px; +line-height:1.2em; +margin:15px 0 20px 35px; +} + +/*** Sidebar and menu ***/ + +#sidebar { +float:left; +line-height:1.4em; +margin:0 0 5px; +padding:1px 0 0; +width:195px; +} + +#sidebar ul { +font-size:0.9em; +list-style:none; +margin:0; +padding:0 0 15px 10px; +} + +#sidebar li { +list-style:none; +margin:0 0 4px; +padding:0; +} + +#sidebar li a { +font-size:1.2em; +font-weight:700; +padding:2px; +} + +#sidebar ul ul { +line-height:1.2em; +margin:4px 0 3px 15px; +padding:0; +} + +#sidebar ul ul li a { +font-weight:400; +} + +#sidebar h2 { +margin:3px 0 8px; +} + +/*** Main content ***/ + +#content { +float:right; +line-height:1.5em; +margin:0; +padding:0; +text-align:left; +width:750px; +} + +#contentalt { +float:left; +line-height:1.5em; +margin-right:20px; +padding:0; +text-align:left; +width:750px; +} + +#content h3,#contentalt h3 { +margin:10px 0 8px; +} + +/*** Footer ***/ + +#footer { +border-top:4px solid #dadada; +clear:both; +color:gray; +font-size:0.9em; +line-height:1.6em; +margin:0 auto; +padding:8px 0; +text-align:right; +} + +#footer p { +margin:0; +padding:0; +} + +#footer a { +color:#808080; +} + +pre { + margin: 10px 30px 10px +} + +/*** Various classes ***/ + +.box { +background:#387C44; +border:1px solid #c8c8c8; +color:#fff; +font-size:0.9em; +line-height:1.4em; +padding:10px 10px 10px 13px; +} + +.box a { +color:#f0f0f0; +} + +.left { +float:left; +margin:0 15px 4px 0; +} + +.right { +float:right; +margin:0 0 4px 15px; +} + +.readmore { +margin:-10px 10px 12px 0; +text-align:right; +} + +.timestamp { +font-size:1.2em; +margin:-5px 0 15px 10px; +} + +.timestamp a { +font-weight:normal; +} + +.green { +color:#387C44; +} + +.clear { +clear:both; +} + +.fade { +color:#c8c8c8; +} + +.gray { +color:gray; +} + +.photo { +background:#fff; +border:1px solid #bababa; +margin:6px 18px 2px 5px; +padding:2px; +} diff --git a/docs/c_delsparse.zip b/docs/c_delsparse.zip new file mode 100644 index 0000000..e06533e Binary files /dev/null and b/docs/c_delsparse.zip differ diff --git a/docs/delsparse.zip b/docs/delsparse.zip new file mode 100644 index 0000000..9e19a1b Binary files /dev/null and b/docs/delsparse.zip differ diff --git a/docs/index.html b/docs/index.html new file mode 100644 index 0000000..3f5ccfa --- /dev/null +++ b/docs/index.html @@ -0,0 +1,137 @@ + + + + + + + +DELAUNAYSPARSE + + + +
+ + + +
+

About the Software

+

+ DELAUNAYSPARSE [1] + contains both serial and parallel codes written in Fortran 2003 + (with OpenMP 4.5) for performing medium- to high-dimensional + interpolation via the Delaunay triangulation. To accommodate + the exponential growth in the size of the Delaunay triangulation in + high dimensions, DELAUNAYSPARSE computes only a sparse subset of the + complete Delaunay triangulation, as necessary for performing + interpolation at the user specified points. + DELAUNAYSPARSE implements the algorithm in [2] and + has been used for various applications including aerospace + engineering [3], HPC performance modeling [4,5], nonparametric + statistics [6], machine learning regression [7], and graph generation [8]. +

+

+ Read the User Guide » +

+ +

Publications

+

+[1] T. H. Chang, L. T. Watson, T. C. H. Lux, A. R. Butt, K. W. Cameron, +and Y. Hong, "Algorithm 1012: DELAUNAYSPARSE: Interpolation via a Sparse Subset +of the Delaunay triangulation in Medium to High Dimensions", +ACM Trans. Math. Software, 46, Article 38 (2020), 1-20. +

 + +

+[2] T. H. Chang, L. T. Watson, T. C. H. Lux, B. Li, L. Xu, A. R. Butt, +K. W. Cameron, and Y. Hong, "A polynomial time algorithm for multivariate +interpolation in arbitrary dimension via the Delaunay triangulation", +in Proc. 2018 ACMSE Conference., Association of Computing Machinery, +Richmond, KY, USA, 2018, Article No. 13. +

 + +

+[3] M. Jrad, R. K. Kapania, J. A. Schetz, and L. T. Watson, +"Self-learning, adaptive software for aerospace engineering applications: +Example of oblique shocks in supersonic flow", +in AIAA Scitech 2019 Forum, American Institute of Aeronautics and +Astronautics, Inc., San Diego, CA, USA, 2019, 1704. +

 + +

+[4] T. H. Chang, L. T. Watson, T. C. H. Lux, J. Bernard, B. Li, L. Xu, +G. Back, A. R. Butt, K. W. Cameron, and Y. Hong, +"Predicting system performance by interpolation using a high-dimensional +Delaunay triangulation", in Proc. SpringSim 2018, the 26th High Performance +Computing Symposium (HPC '18), Society for Modeling and Simulation +International, Baltimore, MD, USA, 2018, Article No. 2. +

 + +

+[5] T. C. H. Lux, L. T. Watson, T. H. Chang, J. Bernard, B. Li, L. Xu, +G. Back, A. R. Butt, K. W. Cameron, and Y. Hong, "Predictive modeling of +I/O characteristics in high performance computing systems", in +Proc. SpringSim 2018, the 26th High Performance Computing +Symposium (HPC '18), Society for Modeling and Simulation International, +Baltimore, MD, USA, 2018, Article No. 8. +

 + +

+[6] T. C. H. Lux, L. T. Watson, T. H. Chang, J. Bernard, B. Li, X. Yu, +L. Xu, A. R. Butt, K. W. Cameron, Y. Hong, and D. Yao, "Nonparametric +distribution models for predicting and managing computational performance +variability", in Proc. IEEE SoutheastCon 2018, Institute of +Electrical and Electronics Engineers, St. Petersburg, FL, USA, 2018, +7 pages. +

 + +

+[7] T. C. H. Lux, L. T. Watson, T. H. Chang, Y. Hong, and K. W. Cameron, +"Interpolation of sparse high-dimensional data", Numerical Algorithms, +(2020) 33 pages. +

 + +

+[8] T. H. Chang, "Mathematical Software for Multiobjective Optimization +Problems". +Ph.D. thesis, Virginia Polytechnic Institute and State University, +Blacksburg, VA, USA, 2020. +

 +
+ + + + + +
+ + diff --git a/docs/py_delsparse.zip b/docs/py_delsparse.zip new file mode 100644 index 0000000..f7605e3 Binary files /dev/null and b/docs/py_delsparse.zip differ diff --git a/docs/usergd.html b/docs/usergd.html new file mode 100644 index 0000000..f84766b --- /dev/null +++ b/docs/usergd.html @@ -0,0 +1,872 @@ + + + + + + + +DELAUNAYSPARSE + + + +
+ + + +
+

Overview

+

+The package +DELAUNAYSPARSE +(ACM TOMS Algorithm 1012) contains serial and parallel codes, written in +FORTRAN 2003 with OpenMP 4.5, for performing +interpolation in medium to high dimensions via a sparse subset of the +Delaunay triangulation. +In addition to the original FORTRAN source code, this site can +be used to download a Python 3.6+ wrapper and C/C++ +bindings for DELAUNAYSPARSE. +Note that each of the three downloads is self-contained, with the +Python wrapper and C bindings each containing +a subset of the FORTRAN as is needed to build and run +their respective functionalities. +Command line drivers, which accept formatted data files, are also available +by downloading the original FORTRAN source code. +

+

+The serial driver subroutine is DELAUNAYSPARSES +and the parallel driver subroutine is DELAUNAYSPARSEP. +The subroutines DELAUNAYSPARSE{S|P} use the module +REAL_PRECISION from HOMPACK90 +(ACM TOMS Algorithm 777) for specifying the real data type, +and the subroutine DWNNLS from SLATEC +to compute projections onto the convex hull during extrapolation. +The master module DELSPARSE_MOD includes the +REAL_PRECISION module and interface blocks for both +DELAUNAYSPARSES and DELAUNAYSPARSEP, as well as +an interface block for the updated SLATEC subroutine +DWNNLS, which may be of separate interest. +Comments at the beginning of the driver subroutines +DELAUNAYSPARSE{S|P} document the arguments and usage, +and examples demonstrating their usage are provided in the +sample programs samples.f90 and samplep.f90. +

+

+The physical organization of the main +TOMS source download into files is as follows. +Further details on using the Python wrapper +and C bindings downloads are given in their +respective README files, which are included in the downloads. + +

    +
  • The file delsparse.f90 contains the module + REAL_PRECISION, DELSPARSE_MOD, and + the driver subroutines DELAUNAYSPARSES, and + DELAUNAYSPARSEP.
    • +
  • The file slatec.f contains the subroutine + DWNNLS and its dependencies from the SLATEC + library. This library has been slightly modified to + comply with the modern Fortran standards. Additionally, legacy + implementations of the BLAS subroutines DROTM + and DTROMG have been included under different names to + avoid dependency issues. Depending on your compiler, you may still + receive warnings related to obsolete features.
    • +
  • The file samples.f90 contains a sample command line program + demonstrating the usage of DELAUNAYSPARSES, with optional + arguments. This program can also be used to use DELAUNAYSPARSES + on formatted data files from the command line.
    • +
  • The file samplep.f90 contains a sample command line program + demonstrating the usage of DELAUNAYSPARSEP, with optional + arguments. This program can also be used to use DELAUNAYSPARSEP + on formatted data files from the command line.
    • +
  • The file test_install.f90 contains a simple test program + that checks whether the installation of DELAUNAYSPARSE appears correct, + based on the output to a small interpolation/extrapolation problem.
    • +
  • The file sample_input2d.dat contains a sample 2-dimensional + input data set for samples.f90 and samplep.f90, + coming from the + VarSys project at Virginia Tech.
    • +
  • The file sample_input4d.dat contains a sample 4-dimensional + input data set for samples.f90 and samplep.f90, + coming from the + VarSys project at Virginia Tech.
    • +
  • The files lapack.f and blas.f contain all + LAPACK and BLAS + subroutines that are referenced (both directly and indirectly) in + DELAUNAYSPARSE.
    • +
  • A sample GNU Makefile is provided. +
+

+ +

+To check that the installation of DELAUNAYSPARSES and +DELAUNAYSPARSEP is correct, assuming that your Fortran compiler +allows mixing fixed format .f and free format .f90 +files in the same compile command, use the command +

+ +
+$FORT $OPTS delsparse.f90 slatec.f lapack.f blas.f test_install.f90 \
+  -o test_install $LIBS
+
+ +

+where $FORT is a FORTRAN 2003 compliant compiler +supporting OpenMP 4.5, $OPTS is a list of compiler +options, and $LIBS is a list of flags to link the +BLAS and LAPACK libraries, if those exist on your +system (in which case the files blas.f and lapack.f +can be omitted from the compile command). To run the parallel code, +$OPTS must include the compiler option for OpenMP. +

+ +

+Then run the tests using +

+ +
+./test_install
+
+ +

+To compile and link the sample main programs sample{s|p}.f90, use +

+ +
+$FORT $OPTS delsparse.f90 slatec.f lapack.f blas.f sample{s|p}.f90 \
+  -o sample{s|p} $LIBS
+
+ +

+similar to above. To run a sample main program, use +

+ +
+./sample{s|p} sample_input{2|4}d.dat
+
+ +

+where sample_input{2|4}d.dat could be replaced by any other +similarly formatted data file. +

+ + +

Usage Information for DELAUNAYSPARSE

+ +

+DELAUNAYSPARSE solves the multivariate interpolation problem: +

+

+Given a set of N points PTS and a set of M +interpolation points Q in R^D, for each interpolation +point Q_i in Q, identify the set +of D+1 data points in PTS that are the vertices of a +Delaunay simplex containing Q_i. +

+

+These vertices can be used to calculate the Delaunay interpolant using +a piecewise linear model. +

+

+For more information on the underlying algorithm, see +

+    @inproceedings{algorithm,
+        author={Chang, Tyler H. and Watson, Layne T. and Lux, Thomas C. H. and
+                Li, Bo and Xu, Li and Butt, Ali R. and Cameron, Kirk W. and
+                Hong, Yili},
+        title={A polynomial time algorithm for multivariate interpolation in
+               arbitrary dimension via the {D}elaunay triangulation},
+        year={2018},
+        month={mar},
+        booktitle={Proc. ACMSE 2018 Conference (ACMSE 18)},
+        publisher={ACM},
+        location={Richmond, KY},
+        doi={10.1145/3190645.3190680}
+    }
+
+

+

+For more information on this software, see +

+    @article{toms1012,
+        author={Chang, Tyler H. and Watson, Layne T. and Lux, Thomas C. H.
+                and Butt, Ali R. and Cameron, Kirk W. and Hong, Yili},
+        title={Algorithm 1012: {DELAUNAYSPARSE}: Interpolation via a sparse
+               subset of the {D}elaunay triangulation in medium to high
+               dimensions},
+        year={2020},
+        volume={46},
+        number={4},
+        articleno={38},
+        nopages={20},
+        doi={10.1145/3422818}
+    }
+
+

+

+DELAUNAYSPARSE contains a Fortran module +

    +
  • delsparse;
    • +
+as well as C bindings +
    +
  • delsparsec;
    • +
+two command-line drivers +
    +
  • delsparses and
    • +
  • delsparsep;
    • +
+ +and a Python 3 wrapper +
    +
  • python.
    • +
+

+

+These interfaces are described in the following sections. +

+ +

Fortran interface

+

+DELAUNAYSPARSE is written in Fortran 2003, and this is its native interface. +The Fortran interface contains two drivers: +

    +
  • DELAUNAYSPARSES (serial driver) and
    • +
  • DELAUNAYSPARSEP (OpenMP parallel driver).
    • +
+

+ +

DELAUNAYSPARSES

+

+The interface for DELAUNAYSPARSES is +

+SUBROUTINE DELAUNAYSPARSES( D, N, PTS, M, Q, SIMPS, WEIGHTS, IERR,     &
+                            INTERP_IN, INTERP_OUT, EPS, EXTRAP, RNORM, &
+                            IBUDGET, CHAIN, EXACT                      )
+
+

+

+Each of the above parameters is described below. +

+

+On input: +

    +
  • D is the dimension of the space for PTS and + Q.
    • +
  • N is the number of data points in PTS.
    • +
  • PTS(1:D,1:N) is a real valued matrix with N + columns, each containing the coordinates of a single data point in + R^D.
    • +
  • M is the number of interpolation points in Q.
    • +
  • Q(1:D,1:M) is a real valued matrix with M columns, + each containing the coordinates of a single interpolation point in + R^D.
    • +
+

+

+On output: +

    +
  • PTS and Q have been rescaled and shifted. + All the data points in PTS are now contained in the unit + hyperball in R^D, and the points in Q + have been shifted and scaled accordingly in relation to PTS. +
    • +
  • SIMPS(1:D+1,1:M) contains the D+1 integer + indices (corresponding to columns in PTS) for the + D+1 vertices of the Delaunay simplex containing each + interpolation point in Q.
    • +
  • WEIGHTS(1:D+1,1:M) contains the D+1 real-valued + weights for expressing each point in Q as a convex combination + of the D+1 corresponding vertices in SIMPS.
    • +
  • IERR(1:M) contains integer valued error flags associated with + the computation of each of the M interpolation points in + Q. The error codes are: +

    + Codes 0, 1, 2 are expected to occur during normal execution. +
      +
    • 00 : Succesful interpolation.
      • +
    • 01 : Succesful extrapolation (up to the allowed extrapolation + distance).
      • +
    • 02 : This point was outside the allowed extrapolation distance; the + corresponding entries in SIMPS and WEIGHTS contain zero values.
      • +
    + + Error codes 10--28 indicate that one or more inputs contain illegal + values or are incompatible with each other. +
      +
    • 10 : The dimension D must be positive.
      • +
    • 11 : Too few data points to construct a triangulation (i.e., + N < D+1).
      • +
    • 12 : No interpolation points given (i.e., M < 1).
      • +
    • 13 : The first dimension of PTS does not agree with the + dimension D.
      • +
    • 14 : The second dimension of PTS does not agree with the + number of points N.
      • +
    • 15 : The first dimension of Q does not agree with the + dimension D.
      • +
    • 16 : The second dimension of Q does not agree with the + number of interpolation points M.
      • +
    • 17 : The first dimension of the output array SIMPS does + not match the number of vertices needed for a D-simplex + (D+1).
      • +
    • 18 : The second dimension of the output array SIMPS does + not match the number of interpolation points M.
      • +
    • 19 : The first dimension of the output array WEIGHTS does + not match the number of vertices for a a D-simplex + (D+1).
      • +
    • 20 : The second dimension of the output array WEIGHTS + does not match the number of interpolation points M.
      • +
    • 21 : The size of the error array IERR does not match the + number of interpolation points M.
      • +
    • 22 : INTERP_IN cannot be present without + INTERP_OUT or vice versa.
      • +
    • 23 : The first dimension of INTERP_IN does not match the + first dimension of INTERP_OUT.
      • +
    • 24 : The second dimension of INTERP_IN does not match the + number of data points PTS.
      • +
    • 25 : The second dimension of INTERP_OUT does not match the + number of interpolation points M.
      • +
    • 26 : The budget supplied in IBUDGET does not contain a + positive integer.
      • +
    • 27 : The extrapolation distance supplied in EXTRAP cannot + be negative.
      • +
    • 28 : The size of the RNORM output array does not match the + number of interpolation points M.
      • +
    + + The errors 30, 31 typically indicate that DELAUNAYSPARSE has been given + an unclean dataset. These errors can be fixed by preprocessing your + data (remove duplicate points and apply PCA or other dimension reduction + technique). +
      +
    • 30 : Two or more points in the data set PTS are too close + together with respect to the working precision (EPS), + which would result in a numerically degenerate simplex.
      • +
    • 31 : All the data points in PTS lie in some lower + dimensional linear manifold (up to the working precision), and no valid + triangulation exists.
      • +
    + + The error code 40 occurs when another earlier error prevented this point + from ever being evaluated. +
      +
    • 40 : An error caused DELAUNAYSPARSES to terminate before + this value could be computed. Note: The corresponding entries in + SIMPS and WEIGHTS may contain garbage values. +
      • +
    + + The error code 50 corresponds to allocation of the internal WORK array. + Check your systems internal memory settings and limits, in relation + to the problem size and DELAUNAYSPARSE's space requirements (see TOMS + Alg. paper for more details on DELAUNAYSPARSE's space requirements). +
      +
    • 50 : A memory allocation error occurred while allocating the work array + WORK.
      • +
    + + The errors 60, 61 should not occur with the default settings. If one of + these errors is observed, then it is likely that either the value of + the optional inputs IBUDGET or EPS has been + adjusted in a way that is unwise, or there may be another issue with the + problem settings, which is manifesting in an unusual way. +
      +
    • 60 : The budget was exceeded before the algorithm converged on this + value. If the dimension is high, try increasing IBUDGET. + This error can also be caused by a working precision EPS + that is too small for the conditioning of the problem.
      • +
    • 61 : A value that was judged appropriate later caused LAPACK to + encounter a singularity. Try increasing the value of EPS. +
      • +
    + + The errors 70--72 were caused by the DWNNLS library from SLATEC, which + is only used during extrapolation. Note that there is a known issue + with this library, when it is linked against included public-domain + copy of BLAS/LAPACK, instead of an installed version + (i.e., -lblas -llapack). +
      +
    • 70 : Allocation error for the extrapolation work arrays.
      • +
    • 71 : The SLATEC subroutine DWNNLS failed to converge + during the projection of an extrapolation point onto the convex hull. +
      • +
    • 72 : The SLATEC subroutine DWNNLS has reported a usage + error.
      • +
    + + The errors 72, 80--83 should never occur, and likely indicate a + compiler bug or hardware failure. +
      +
    • 80 : The LAPACK subroutine DGEQP3 has reported an illegal + value.
      • +
    • 81 : The LAPACK subroutine DGETRF has reported an illegal + value.
      • +
    • 82 : The LAPACK subroutine DGETRS has reported an illegal + value.
      • +
    • 83 : The LAPACK subroutine DORMQR has reported an illegal + value.
      • +
    +
+

+

+Optional arguments: +

    +
  • INTERP_IN(1:IR,1:N) contains real valued response vectors for + each of the data points in PTS on input. The first dimension + of INTERP_IN is inferred to be the dimension of these + response vectors, and the second dimension must match N. + If present, the response values will be computed for each interpolation + point in Q, and stored in INTERP_OUT, + which therefore must also be present. If both INTERP_IN and + INTERP_OUT are omitted, only the containing simplices and + convex combination weights are returned.
    • +
  • INTERP_OUT(1:IR,1:M) contains real valued response vectors + for each interpolation point in Q on output. The first + dimension of INTERP_OUT must match the first dimension of + INTERP_IN, and the second dimension must match M. + If present, the response values at each interpolation point are computed + as a convex combination of the response values (supplied in + INTERP_IN) at the vertices of a Delaunay simplex containing + that interpolation point. Therefore, if INTERP_OUT is + present, then INTERP_IN must also be present. If both are + omitted, only the simplices and convex combination weights are returned. +
    • +
  • EPS contains the real working precision for the problem on + input. By default, EPS is assigned \sqrt{\mu} + where \mu denotes the unit roundoff for the machine. In + general, any values that differ by less than EPS are judged + as equal, and any weights that are greater than -EPS are + judged as nonnegative. EPS cannot take a value less than the + default value of \sqrt{\mu}. If any value less than + \sqrt{\mu} is supplied, the default value will be used + instead automatically.
    • +
  • EXTRAP contains the real maximum extrapolation distance + (relative to the diameter of PTS) on input. Interpolation + at a point outside the convex hull of PTS is done by + projecting that point onto the convex hull, and then doing normal + Delaunay interpolation at that projection. Interpolation at any point + in Q that is more than EXTRAP * DIAMETER(PTS) + units outside the convex hull of PTS will not be done and + an error code of 2 will be returned. Note that computing + the projection can be expensive. Setting EXTRAP=0 will + cause all extrapolation points to be ignored without ever computing a + projection. By default, EXTRAP=0.1 + (extrapolate by up to 10% of the diameter of PTS).
    • +
  • RNORM(1:M) contains the real unscaled projection (2-norm) + distances from any projection computations on output. If not present, + these distances are still computed for each extrapolation point, but are + never returned.
    • +
  • IBUDGET on input contains the integer budget for performing + flips while iterating toward the simplex containing each interpolation + point in Q. This prevents DELAUNAYSPARSES from + falling into an infinite loop when an inappropriate value of + EPS is given with respect to the problem conditioning. + By default, IBUDGET=50000. However, for extremely + high-dimensional problems and pathological inputs, the default value + may be insufficient.
    • +
  • CHAIN is a logical input argument that determines whether a + new first simplex should be constructed for each interpolation point + (CHAIN=.FALSE.), or whether the simplex walks should be + "daisy-chained." By default, CHAIN=.FALSE. Setting + CHAIN=.TRUE. is generally not recommended, unless the size + of the triangulation is relatively small or the interpolation points are + known to be tightly clustered.
    • +
  • EXACT is a logical input argument that determines whether + the exact diameter should be computed and whether a check for duplicate + data points should be performed in advance. When + EXACT=.FALSE., the diameter of PTS is + approximated by twice the distance from the barycenter of PTS + to the farthest point in PTS, and no check is done to find + the closest pair of points, which could result in hard to find bugs later + on. When EXACT=.TRUE., the exact diameter is computed and + an error is returned whenever PTS contains duplicate values up to the + precision EPS. By default EXACT=.TRUE., but + setting EXACT=.FALSE. could result in significant speedup + when N is large.
    • It is strongly recommended that most + users leave EXACT=.TRUE., as setting + EXACT=.FALSE. could result in input errors that are difficult + to identify. Also, the diameter approximation could be wrong by up to + a factor of two. +
+

+

+Subroutines and functions directly referenced from BLAS are +

    +
  • DDOT,
    • +
  • DGEMV,
    • +
  • DNRM2,
    • +
  • DTRSM,
    • +
+and from LAPACK are +
    +
  • DGEQP3,
    • +
  • DGETRF,
    • +
  • DGETRS,
    • +
  • DORMQR.
    • +
+

+

+The SLATEC subroutine +

    +
  • DWNNLS is also directly referenced.
    • +
+DWNNLS and all its SLATEC dependencies have been slightly edited +to comply with the Fortran 2008 standard, with all print statements and +references to stderr being commented out. For a reference to DWNNLS, +see ACM TOMS Algorithm 587 (Hanson and Haskell). +The module REAL_PRECISION from HOMPACK90 (ACM TOMS Algorithm 777) +is used for the real data type. The REAL_PRECISION module, +DELAUNAYSPARSES, and DWNNLS and its dependencies +comply with the Fortran 2008 standard. +

+ +

DELAUNAYSPARSEP

+

+The interface for DELAUNAYSPARSEP is +

+    SUBROUTINE DELAUNAYSPARSEP( D, N, PTS, M, Q, SIMPS, WEIGHTS, IERR,     &
+                                INTERP_IN, INTERP_OUT, EPS, EXTRAP, RNORM, &
+                                IBUDGET, CHAIN, EXACT, PMODE               )
+
+

+

+Each of the above parameters is described below. +

+

+On input: +

    +
  • D is the dimension of the space for PTS and + Q.
    • +
  • N is the number of data points in PTS.
    • +
  • PTS(1:D,1:N) is a real valued matrix with N + columns, each containing the coordinates of a single data point in + R^D.
    • +
  • M is the number of interpolation points in Q.
    • +
  • Q(1:D,1:M) is a real valued matrix with M columns, + each containing the coordinates of a single interpolation point in + R^D.
    • +
+

+

+On output: +

    +
  • PTS and Q have been rescaled and shifted. + All the data points in PTS are now contained in the unit + hyperball in R^D, and the points in Q + have been shifted and scaled accordingly in relation to PTS. +
    • +
  • SIMPS(1:D+1,1:M) contains the D+1 integer + indices (corresponding to columns in PTS) for the + D+1 vertices of the Delaunay simplex containing each + interpolation point in Q.
    • +
  • WEIGHTS(1:D+1,1:M) contains the D+1 real-valued + weights for expressing each point in Q as a convex combination + of the D+1 corresponding vertices in SIMPS.
    • +
  • IERR(1:M) contains integer valued error flags associated with + the computation of each of the M interpolation points in + Q. The error codes are: +

    + Codes 0, 1, 2 are expected to occur during normal execution. +
      +
    • 00 : Succesful interpolation.
      • +
    • 01 : Succesful extrapolation (up to the allowed extrapolation + distance).
      • +
    • 02 : This point was outside the allowed extrapolation distance; the + corresponding entries in SIMPS and WEIGHTS contain zero values.
      • +
    + + Error codes 10--28 indicate that one or more inputs contain illegal + values or are incompatible with each other. +
      +
    • 10 : The dimension D must be positive.
      • +
    • 11 : Too few data points to construct a triangulation (i.e., + N < D+1).
      • +
    • 12 : No interpolation points given (i.e., M < 1).
      • +
    • 13 : The first dimension of PTS does not agree with the + dimension D.
      • +
    • 14 : The second dimension of PTS does not agree with the + number of points N.
      • +
    • 15 : The first dimension of Q does not agree with the + dimension D.
      • +
    • 16 : The second dimension of Q does not agree with the + number of interpolation points M.
      • +
    • 17 : The first dimension of the output array SIMPS does + not match the number of vertices needed for a D-simplex + (D+1).
      • +
    • 18 : The second dimension of the output array SIMPS does + not match the number of interpolation points M.
      • +
    • 19 : The first dimension of the output array WEIGHTS does + not match the number of vertices for a a D-simplex + (D+1).
      • +
    • 20 : The second dimension of the output array WEIGHTS + does not match the number of interpolation points M.
      • +
    • 21 : The size of the error array IERR does not match the + number of interpolation points M.
      • +
    • 22 : INTERP_IN cannot be present without + INTERP_OUT or vice versa.
      • +
    • 23 : The first dimension of INTERP_IN does not match the + first dimension of INTERP_OUT.
      • +
    • 24 : The second dimension of INTERP_IN does not match the + number of data points PTS.
      • +
    • 25 : The second dimension of INTERP_OUT does not match the + number of interpolation points M.
      • +
    • 26 : The budget supplied in IBUDGET does not contain a + positive integer.
      • +
    • 27 : The extrapolation distance supplied in EXTRAP cannot + be negative.
      • +
    • 28 : The size of the RNORM output array does not match the + number of interpolation points M.
      • +
    + + The errors 30, 31 typically indicate that DELAUNAYSPARSE has been given + an unclean dataset. These errors can be fixed by preprocessing your + data (remove duplicate points and apply PCA or other dimension reduction + technique). +
      +
    • 30 : Two or more points in the data set PTS are too close + together with respect to the working precision (EPS), + which would result in a numerically degenerate simplex.
      • +
    • 31 : All the data points in PTS lie in some lower + dimensional linear manifold (up to the working precision), and no valid + triangulation exists.
      • +
    + + The error code 40 occurs when another earlier error prevented this point + from ever being evaluated. +
      +
    • 40 : An error caused DELAUNAYSPARSEP to terminate before + this value could be computed. Note: The corresponding entries in + SIMPS and WEIGHTS may contain garbage values. +
      • +
    + + The error code 50 corresponds to allocation of the internal WORK array. + Check your systems internal memory settings and limits, in relation + to the problem size and DELAUNAYSPARSE's space requirements (see TOMS + Alg. paper for more details on DELAUNAYSPARSE's space requirements). +
      +
    • 50 : A memory allocation error occurred while allocating the work array + WORK.
      • +
    + + The errors 60, 61 should not occur with the default settings. If one of + these errors is observed, then it is likely that either the value of + the optional inputs IBUDGET or EPS has been + adjusted in a way that is unwise, or there may be another issue with the + problem settings, which is manifesting in an unusual way. +
      +
    • 60 : The budget was exceeded before the algorithm converged on this + value. If the dimension is high, try increasing IBUDGET. + This error can also be caused by a working precision EPS + that is too small for the conditioning of the problem.
      • +
    • 61 : A value that was judged appropriate later caused LAPACK to + encounter a singularity. Try increasing the value of EPS. +
      • +
    + + The errors 70--72 were caused by the DWNNLS library from SLATEC, which + is only used during extrapolation. Note that there is a known issue + with this library, when it is linked against included public-domain + copy of BLAS/LAPACK, instead of an installed version + (i.e., -lblas -llapack). +
      +
    • 70 : Allocation error for the extrapolation work arrays.
      • +
    • 71 : The SLATEC subroutine DWNNLS failed to converge + during the projection of an extrapolation point onto the convex hull. +
      • +
    • 72 : The SLATEC subroutine DWNNLS has reported a usage + error.
      • +
    + + The errors 72, 80--83 should never occur, and likely indicate a + compiler bug or hardware failure. +
      +
    • 80 : The LAPACK subroutine DGEQP3 has reported an illegal + value.
      • +
    • 81 : The LAPACK subroutine DGETRF has reported an illegal + value.
      • +
    • 82 : The LAPACK subroutine DGETRS has reported an illegal + value.
      • +
    • 83 : The LAPACK subroutine DORMQR has reported an illegal + value.
      • +
    + + The error code 90 is unique to DELAUNAYSPARSEP. +
      +
    • 90 : The value of PMODE is not valid.
      • +
    +
+

+

+Optional arguments: +

    +
  • INTERP_IN(1:IR,1:N) contains real valued response vectors for + each of the data points in PTS on input. The first dimension + of INTERP_IN is inferred to be the dimension of these + response vectors, and the second dimension must match N. + If present, the response values will be computed for each interpolation + point in Q, and stored in INTERP_OUT, + which therefore must also be present. If both INTERP_IN and + INTERP_OUT are omitted, only the containing simplices and + convex combination weights are returned.
    • +
  • INTERP_OUT(1:IR,1:M) contains real valued response vectors + for each interpolation point in Q on output. The first + dimension of INTERP_OUT must match the first dimension of + INTERP_IN, and the second dimension must match M. + If present, the response values at each interpolation point are computed + as a convex combination of the response values (supplied in + INTERP_IN) at the vertices of a Delaunay simplex containing + that interpolation point. Therefore, if INTERP_OUT is + present, then INTERP_IN must also be present. If both are + omitted, only the simplices and convex combination weights are returned. +
    • +
  • EPS contains the real working precision for the problem on + input. By default, EPS is assigned \sqrt{\mu} + where \mu denotes the unit roundoff for the machine. In + general, any values that differ by less than EPS are judged + as equal, and any weights that are greater than -EPS are + judged as nonnegative. EPS cannot take a value less than the + default value of \sqrt{\mu}. If any value less than + \sqrt{\mu} is supplied, the default value will be used + instead automatically.
    • +
  • EXTRAP contains the real maximum extrapolation distance + (relative to the diameter of PTS) on input. Interpolation + at a point outside the convex hull of PTS is done by + projecting that point onto the convex hull, and then doing normal + Delaunay interpolation at that projection. Interpolation at any point + in Q that is more than EXTRAP * DIAMETER(PTS) + units outside the convex hull of PTS will not be done and + an error code of 2 will be returned. Note that computing + the projection can be expensive. Setting EXTRAP=0 will + cause all extrapolation points to be ignored without ever computing a + projection. By default, EXTRAP=0.1 + (extrapolate by up to 10% of the diameter of PTS).
    • +
  • RNORM(1:M) contains the real unscaled projection (2-norm) + distances from any projection computations on output. If not present, + these distances are still computed for each extrapolation point, but are + never returned.
    • +
  • IBUDGET on input contains the integer budget for performing + flips while iterating toward the simplex containing each interpolation + point in Q. This prevents DELAUNAYSPARSEP from + falling into an infinite loop when an inappropriate value of + EPS is given with respect to the problem conditioning. + By default, IBUDGET=50000. However, for extremely + high-dimensional problems and pathological inputs, the default value + may be insufficient.
    • +
  • CHAIN is a logical input argument that determines whether a + new first simplex should be constructed for each interpolation point + (CHAIN=.FALSE.), or whether the simplex walks should be + "daisy-chained." By default, CHAIN=.FALSE. Setting + CHAIN=.TRUE. is generally not recommended, unless the size + of the triangulation is relatively small or the interpolation points are + known to be tightly clustered.
    • +
  • EXACT is a logical input argument that determines whether + the exact diameter should be computed and whether a check for duplicate + data points should be performed in advance. When + EXACT=.FALSE., the diameter of PTS is + approximated by twice the distance from the barycenter of PTS + to the farthest point in PTS, and no check is done to find + the closest pair of points, which could result in hard to find bugs later + on. When EXACT=.TRUE., the exact diameter is computed and + an error is returned whenever PTS contains duplicate values up to the + precision EPS. By default EXACT=.TRUE., but + setting EXACT=.FALSE. could result in significant speedup + when N is large.
    • It is strongly recommended that most + users leave EXACT=.TRUE., as setting + EXACT=.FALSE. could result in input errors that are difficult + to identify. Also, the diameter approximation could be wrong by up to + a factor of two. +
  • PMODE is an integer specifying the level of parallelism to + be exploited. +
      +
    • If PMODE = 1, then parallelism is exploited at the + level of the loop over all interpolation points (Level 1 + parallelism).
      • +
    • If PMODE = 2, then parallelism is exploited at the + level of the loops over data points when constructing/flipping + simplices (Level 2 parallelism).
      • +
    • If PMODE = 3, then parallelism is exploited at both + levels. Note: this implies that the total number of threads active + at any time could be up to OMP_NUM_THREADS^2. + By default, PMODE is set to 1 if there + is more than 1 interpolation point and 2 otherwise. +
      • +
    • +
+

+

+Subroutines and functions directly referenced from BLAS are +

    +
  • DDOT,
    • +
  • DGEMV,
    • +
  • DNRM2,
    • +
  • DTRSM,
    • +
+and from LAPACK are +
    +
  • DGEQP3,
    • +
  • DGETRF,
    • +
  • DGETRS,
    • +
  • DORMQR.
    • +
+

+

+The SLATEC subroutine +

    +
  • DWNNLS is also directly referenced.
    • +
+DWNNLS and all its SLATEC dependencies have been slightly edited +to comply with the Fortran 2008 standard, with all print statements and +references to stderr being commented out. For a reference to DWNNLS, +see ACM TOMS Algorithm 587 (Hanson and Haskell). +The module REAL_PRECISION from HOMPACK90 (ACM TOMS Algorithm 777) +is used for the real data type. The REAL_PRECISION module, +DELAUNAYSPARSEP, and DWNNLS and its dependencies +comply with the Fortran 2008 standard. +

+ +
+ +

+Notes: DELAUNAYSPARSE is available free of charge via a permissive +MIT LICENSE. +

+
+ + + + + +
+ + diff --git a/python/delsparse.py b/python/delsparse.py index 3414f65..e25564d 100644 --- a/python/delsparse.py +++ b/python/delsparse.py @@ -14,11 +14,11 @@ # ^^ 'fPIC' and 'shared' are required. 'O3' is for speed and 'fopenmp' # is necessary for supporting CPU parallelism during execution. blas_lapack = "-lblas -llapack" -blas_lapack = "blas.f lapack.f" +blas_lapack = "delsparse_src/blas.f delsparse_src/lapack.f" # ^^ Use a local BLAS and LAPACK if available by commenting the second line # above. The included "blas.f" and "lapack.f" are known to cause error 71 # during extrapolation, but there is no known resolution. -ordered_dependencies = "real_precision.f90 slatec.f delsparse.f90 delsparse_bind_c.f90" +ordered_dependencies = "delsparse_src/real_precision.f90 delsparse_src/slatec.f delsparse_src/delsparse.f90 delsparse_src/delsparse_bind_c.f90" # # -------------------------------------------------------------------- diff --git a/python/blas.f b/python/delsparse_src/blas.f similarity index 100% rename from python/blas.f rename to python/delsparse_src/blas.f diff --git a/python/delsparse.f90 b/python/delsparse_src/delsparse.f90 similarity index 100% rename from python/delsparse.f90 rename to python/delsparse_src/delsparse.f90 diff --git a/python/delsparse_bind_c.f90 b/python/delsparse_src/delsparse_bind_c.f90 similarity index 100% rename from python/delsparse_bind_c.f90 rename to python/delsparse_src/delsparse_bind_c.f90 diff --git a/python/lapack.f b/python/delsparse_src/lapack.f similarity index 100% rename from python/lapack.f rename to python/delsparse_src/lapack.f diff --git a/python/real_precision.f90 b/python/delsparse_src/real_precision.f90 similarity index 100% rename from python/real_precision.f90 rename to python/delsparse_src/real_precision.f90 diff --git a/python/slatec.f b/python/delsparse_src/slatec.f similarity index 98% rename from python/slatec.f rename to python/delsparse_src/slatec.f index 7d51578..c652a26 100755 --- a/python/slatec.f +++ b/python/delsparse_src/slatec.f @@ -7,7 +7,7 @@ SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, C a covariance matrix. C***LIBRARY SLATEC C***CATEGORY K1A2A, D9 -C***TYPE REAL(KIND=R8) (LSEI-S, DLSEI-D) +C***TYPE DOUBLE PRECISION (LSEI-S, DLSEI-D) C***KEYWORDS CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING, C EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS, C QUADRATIC PROGRAMMING @@ -62,7 +62,7 @@ SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, C C The parameters for DLSEI( ) are C -C Input.. All TYPE REAL variables are REAL(KIND=R8) +C Input.. All TYPE REAL variables are DOUBLE PRECISION C C W(*,*),MDW, The array W(*,*) is doubly subscripted with C ME,MA,MG,N first dimensioning parameter equal to MDW. @@ -268,7 +268,7 @@ SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, C LIP = MG+2*N+2 C This test will not be made if IP(2).LE.0. C -C Output.. All TYPE REAL variables are REAL(KIND=R8) +C Output.. All TYPE REAL variables are DOUBLE PRECISION C C X(*),RNORME, The array X(*) contains the solution parameters C RNORML if the integer output flag MODE = 0 or 1. @@ -382,18 +382,17 @@ SUBROUTINE DLSEI (W, MDW, ME, MA, MG, N, PRGOPT, X, RNORME, C 900510 Convert XERRWV calls to XERMSG calls. (RWC) C 900604 DP version created from SP version. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) -C 180613 Removed prints and replaced DP --> REAL(KIND=R8). (THC) +C 180613 Removed prints and replaced DP --> DOUBLE PRECISION. (THC) C***END PROLOGUE DLSEI - USE REAL_PRECISION INTEGER IP(3), MA, MDW, ME, MG, MODE, N - REAL(KIND=R8) PRGOPT(*), RNORME, RNORML, W(MDW,*), WS(*), X(*) + DOUBLE PRECISION PRGOPT(*), RNORME, RNORML, W(MDW,*), WS(*), X(*) C EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DDOT, DH12, DLSI, DNRM2, * DSCAL, DSWAP - REAL(KIND=R8) D1MACH, DASUM, DDOT, DNRM2 + DOUBLE PRECISION D1MACH, DASUM, DDOT, DNRM2 C - REAL(KIND=R8) DRELPR, ENORM, FNORM, GAM, RB, RN, RNMAX, SIZE, + DOUBLE PRECISION DRELPR, ENORM, FNORM, GAM, RB, RN, RNMAX, SIZE, * SN, SNMAX, T, TAU, UJ, UP, VJ, XNORM, XNRME INTEGER I, IMAX, J, JP1, K, KEY, KRANKE, LAST, LCHK, LINK, M, * MAPKE1, MDEQC, MEND, MEP1, N1, N2, NEXT, NLINK, NOPT, NP1, @@ -743,7 +742,7 @@ SUBROUTINE DLSI (W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, C***SUBSIDIARY C***PURPOSE Subsidiary to DLSEI C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (LSI-S, DLSI-D) +C***TYPE DOUBLE PRECISION (LSI-S, DLSI-D) C***AUTHOR Hanson, R. J., (SNLA) C***DESCRIPTION C @@ -795,16 +794,15 @@ SUBROUTINE DLSI (W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, C 900604 DP version created from SP version. (RWC) C 920422 Changed CALL to DHFTI to include variable MA. (WRB) C***END PROLOGUE DLSI - USE REAL_PRECISION INTEGER IP(*), MA, MDW, MG, MODE, N - REAL(KIND=R8) PRGOPT(*), RNORM, W(MDW,*), WS(*), X(*) + DOUBLE PRECISION PRGOPT(*), RNORM, W(MDW,*), WS(*), X(*) C EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DDOT, DH12, DHFTI, DLPDP, * DSCAL, DSWAP - REAL(KIND=R8) D1MACH, DASUM, DDOT + DOUBLE PRECISION D1MACH, DASUM, DDOT C - REAL(KIND=R8) ANORM, DRELPR, FAC, GAM, RB, TAU, TOL, XNORM, + DOUBLE PRECISION ANORM, DRELPR, FAC, GAM, RB, TAU, TOL, XNORM, * TMP_NORM(1) INTEGER I, J, K, KEY, KRANK, KRM1, KRP1, L, LAST, LINK, M, MAP1, * MDLPDP, MINMAN, N1, N2, N3, NEXT, NP1 @@ -1079,12 +1077,12 @@ SUBROUTINE DLSI (W, MDW, MA, MG, N, PRGOPT, X, RNORM, MODE, WS, RETURN END *DECK D1MACH - REAL(KIND=R8) FUNCTION D1MACH (I) + DOUBLE PRECISION FUNCTION D1MACH (I) C***BEGIN PROLOGUE D1MACH C***PURPOSE Return floating point machine dependent constants. C***LIBRARY SLATEC C***CATEGORY R1 -C***TYPE REAL(KIND=R8) (R1MACH-S, D1MACH-D) +C***TYPE DOUBLE PRECISION (R1MACH-S, D1MACH-D) C***KEYWORDS MACHINE CONSTANTS C***AUTHOR Fox, P. A., (Bell Labs) C Hall, A. D., (Bell Labs) @@ -1151,7 +1149,6 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C comments below. (DWL) C***END PROLOGUE D1MACH C - USE REAL_PRECISION INTEGER SMALL(4) INTEGER LARGE(4) @@ -1164,7 +1161,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C for DMACH(2) is a slight lower bound. The value for DMACH(5) is C a 20-digit approximation. If one of the sets of initial data below C is preferred, do the necessary commenting and uncommenting. (DWL) - REAL(KIND=R8) DMACH(5) + DOUBLE PRECISION DMACH(5) DATA DMACH / 2.23D-308, 1.79D+308, 1.111D-16, 2.222D-16, 1 0.30102999566398119521D0 / SAVE DMACH @@ -1387,7 +1384,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C DATA LOG10(1), LOG10(2) / Z44133FF3, Z79FF509F / C C MACHINE CONSTANTS FOR THE ELXSI 6400 -C (ASSUMING REAL*8 IS THE DEFAULT REAL(KIND=R8)) +C (ASSUMING REAL*8 IS THE DEFAULT DOUBLE PRECISION) C C DATA SMALL(1), SMALL(2) / '00100000'X,'00000000'X / C DATA LARGE(1), LARGE(2) / '7FEFFFFF'X,'FFFFFFFF'X / @@ -1420,7 +1417,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C DATA DMACH(5) / Z'3FD34413509F79FF' / C C MACHINE CONSTANTS FOR THE HP 2100 -C THREE WORD REAL(KIND=R8) OPTION WITH FTN4 +C THREE WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2), SMALL(3) / 40000B, 0, 1 / C DATA LARGE(1), LARGE(2), LARGE(3) / 77777B, 177777B, 177776B / @@ -1429,7 +1426,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C DATA LOG10(1), LOG10(2), LOG10(3) / 46420B, 46502B, 77777B / C C MACHINE CONSTANTS FOR THE HP 2100 -C FOUR WORD REAL(KIND=R8) OPTION WITH FTN4 +C FOUR WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA SMALL(1), SMALL(2) / 40000B, 0 / C DATA SMALL(3), SMALL(4) / 0, 1 / @@ -1461,7 +1458,7 @@ REAL(KIND=R8) FUNCTION D1MACH (I) C DATA LOG10(1), LOG10(2) / Z41134413, Z509F79FF / C C MACHINE CONSTANTS FOR THE IBM PC -C ASSUMES THAT ALL ARITHMETIC IS DONE IN REAL(KIND=R8) +C ASSUMES THAT ALL ARITHMETIC IS DONE IN DOUBLE PRECISION C ON 8088, I.E., NOT IN 80 BIT FORM FOR THE 8087. C C DATA SMALL(1) / 2.23D-308 / @@ -2176,7 +2173,7 @@ INTEGER FUNCTION I1MACH (I) C DATA IMACH(16) / 1024 / C C MACHINE CONSTANTS FOR THE HP 2100 -C 3 WORD REAL(KIND=R8) OPTION WITH FTN4 +C 3 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / @@ -2196,7 +2193,7 @@ INTEGER FUNCTION I1MACH (I) C DATA IMACH(16) / 127 / C C MACHINE CONSTANTS FOR THE HP 2100 -C 4 WORD REAL(KIND=R8) OPTION WITH FTN4 +C 4 WORD DOUBLE PRECISION OPTION WITH FTN4 C C DATA IMACH( 1) / 5 / C DATA IMACH( 2) / 6 / @@ -2507,11 +2504,11 @@ SUBROUTINE DH12 (MODE, LPIVOT, L1, M, U, IUE, UP, C, ICE, ICV, C***SUBSIDIARY C***PURPOSE Subsidiary to DHFTI, DLSEI and DWNNLS C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (H12-S, DH12-D) +C***TYPE DOUBLE PRECISION (H12-S, DH12-D) C***AUTHOR (UNKNOWN) C***DESCRIPTION C -C *** REAL(KIND=R8) VERSION OF H12 ****** +C *** DOUBLE PRECISION VERSION OF H12 ****** C C C.L.Lawson and R.J.Hanson, Jet Propulsion Laboratory, 1973 Jun 12 C to appear in 'Solving Least Squares Problems', Prentice-Hall, 1974 @@ -2548,13 +2545,12 @@ SUBROUTINE DH12 (MODE, LPIVOT, L1, M, U, IUE, UP, C, ICE, ICV, C 890831 Modified array declarations. (WRB) C 891214 Prologue converted to Version 4.0 format. (BAB) C 900328 Added TYPE section. (WRB) -C 900911 Added DDOT to REAL(KIND=R8) statement. (WRB) +C 900911 Added DDOT to DOUBLE PRECISION statement. (WRB) C***END PROLOGUE DH12 - USE REAL_PRECISION INTEGER I, I2, I3, I4, ICE, ICV, INCR, IUE, J, KL1, KL2, KLP, * L1, L1M1, LPIVOT, M, MML1P2, MODE, NCV - REAL(KIND=R8) B, C, CL, CLINV, ONE, UL1M1, SM, U, UP, DDOT + DOUBLE PRECISION B, C, CL, CLINV, ONE, UL1M1, SM, U, UP, DDOT DIMENSION U(IUE,*), C(*) C BEGIN BLOCK PERMITTING ...EXITS TO 140 C***FIRST EXECUTABLE STATEMENT DH12 @@ -2654,7 +2650,7 @@ SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, C Exactly one right-hand side vector is permitted. C***LIBRARY SLATEC C***CATEGORY D9 -C***TYPE REAL(KIND=R8) (HFTI-S, DHFTI-D) +C***TYPE DOUBLE PRECISION (HFTI-S, DHFTI-D) C***KEYWORDS CURVE FITTING, LEAST SQUARES C***AUTHOR Lawson, C. L., (JPL) C Hanson, R. J., (SNLA) @@ -2711,7 +2707,7 @@ SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, C C The entire set of parameters for DHFTI are C -C INPUT.. All TYPE REAL variables are REAL(KIND=R8) +C INPUT.. All TYPE REAL variables are DOUBLE PRECISION C C A(*,*),MDA,M,N The array A(*,*) initially contains the M by N C matrix A of the least squares problem AX = B. @@ -2742,7 +2738,7 @@ SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, C C H(*),G(*),IP(*) Arrays of working space used by DHFTI. C -C OUTPUT.. All TYPE REAL variables are REAL(KIND=R8) +C OUTPUT.. All TYPE REAL variables are DOUBLE PRECISION C C A(*,*) The contents of the array A(*,*) will be C modified by the subroutine. These contents @@ -2782,11 +2778,10 @@ SUBROUTINE DHFTI (A, MDA, M, N, B, MDB, NB, TAU, KRANK, RNORM, H, C 901005 Replace usage of DDIFF with usage of D1MACH. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) C***END PROLOGUE DHFTI - USE REAL_PRECISION INTEGER I, II, IOPT, IP(*), IP1, J, JB, JJ, K, KP1, KRANK, L, * LDIAG, LMAX, M, MDA, MDB, N, NB, NERR - REAL(KIND=R8) A, B, D1MACH, DZERO, FACTOR, + DOUBLE PRECISION A, B, D1MACH, DZERO, FACTOR, * G, H, HMAX, RELEPS, RNORM, SM, SM1, SZERO, TAU, TMP DIMENSION A(MDA,*),B(MDB,*),H(*),G(*),RNORM(*) SAVE RELEPS @@ -2985,7 +2980,7 @@ SUBROUTINE DLPDP (A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS, C***SUBSIDIARY C***PURPOSE Subsidiary to DLSEI C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (LPDP-S, DLPDP-D) +C***TYPE DOUBLE PRECISION (LPDP-S, DLPDP-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -3028,12 +3023,11 @@ SUBROUTINE DLPDP (A, MDA, M, N1, N2, PRGOPT, X, WNORM, MODE, WS, C 900328 Added TYPE section. (WRB) C 910408 Updated the AUTHOR section. (WRB) C***END PROLOGUE DLPDP - USE REAL_PRECISION C INTEGER I, IS(*), IW, IX, J, L, M, MDA, MODE, MODEW, N, N1, N2, * NP1 - REAL(KIND=R8) A(MDA,*), DDOT, DNRM2, FAC, ONE, + DOUBLE PRECISION A(MDA,*), DDOT, DNRM2, FAC, ONE, * PRGOPT(*), RNORM, SC, WNORM, WS(*), X(*), YNORM, ZERO SAVE ZERO, ONE, FAC DATA ZERO,ONE /0.0D0,1.0D0/, FAC /0.1D0/ @@ -3197,7 +3191,7 @@ SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, C selected variables. C***LIBRARY SLATEC C***CATEGORY K1A2A -C***TYPE REAL(KIND=R8) (WNNLS-S, DWNNLS-D) +C***TYPE DOUBLE PRECISION (WNNLS-S, DWNNLS-D) C***KEYWORDS CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING, C EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS, C NONNEGATIVITY CONSTRAINTS, QUADRATIC PROGRAMMING @@ -3232,7 +3226,7 @@ SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, C C The parameters for DWNNLS are C -C INPUT.. All TYPE REAL variables are REAL(KIND=R8) +C INPUT.. All TYPE REAL variables are DOUBLE PRECISION C C W(*,*),MDW, The array W(*,*) is double subscripted with first C ME,MA,N,L dimensioning parameter equal to MDW. For this @@ -3392,7 +3386,7 @@ SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, C LIW = ME+MA+N C This test will not be made if IWORK(2).LE.0. C -C OUTPUT.. All TYPE REAL variables are REAL(KIND=R8) +C OUTPUT.. All TYPE REAL variables are DOUBLE PRECISION C C X(*) An array dimensioned at least N, which will C contain the N components of the solution vector @@ -3455,13 +3449,12 @@ SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, C 900510 Convert XERRWV calls to XERMSG calls, change Prologue C comments to agree with WNNLS. (RWC) C 920501 Reformatted the REFERENCES section. (WRB) -C 180613 Removed prints and replaced DP --> REAL(KIND=R8). (THC) +C 180613 Removed prints and replaced DP --> DOUBLE PRECISION. (THC) C***END PROLOGUE DWNNLS - USE REAL_PRECISION INTEGER IWORK(*), L, L1, L2, L3, L4, L5, LIW, LW, MA, MDW, ME, * MODE, N - REAL(KIND=R8) PRGOPT(*), RNORM, W(MDW,*), WORK(*), X(*) + DOUBLE PRECISION PRGOPT(*), RNORM, W(MDW,*), WORK(*), X(*) C CHARACTER*8 XERN1 C***FIRST EXECUTABLE STATEMENT DWNNLS MODE = 0 @@ -3525,7 +3518,7 @@ SUBROUTINE DWNLSM (W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE, C***SUBSIDIARY C***PURPOSE Subsidiary to DWNNLS C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLSM-S, DWNLSM-D) +C***TYPE DOUBLE PRECISION (WNLSM-S, DWNLSM-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -3539,7 +3532,7 @@ SUBROUTINE DWNLSM (W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE, C sequence from DWNNLS for purposes of variable dimensioning). C Their contents will in general be of no interest to the user. C -C Variables of type REAL are REAL(KIND=R8). +C Variables of type REAL are DOUBLE PRECISION. C C IPIVOT(*) C An array of length N. Upon completion it contains the @@ -3590,18 +3583,17 @@ SUBROUTINE DWNLSM (W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE, C 900604 DP version created from SP version. (RWC) C 900911 Restriction on value of ALAMDA included. (WRB) C***END PROLOGUE DWNLSM - USE REAL_PRECISION INTEGER IPIVOT(*), ITYPE(*), L, MA, MDW, MME, MODE, N - REAL(KIND=R8) D(*), H(*), PRGOPT(*), RNORM, SCALE(*), TEMP(*), + DOUBLE PRECISION D(*), H(*), PRGOPT(*), RNORM, SCALE(*), TEMP(*), * W(MDW,*), WD(*), X(*), Z(*) C EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DH12, DNRM2, SLATEC_DROTM, * SLATEC_DROTMG, DSCAL, DSWAP, DWNLIT, IDAMAX, XERMSG - REAL(KIND=R8) D1MACH, DASUM, DNRM2 + DOUBLE PRECISION D1MACH, DASUM, DNRM2 INTEGER IDAMAX C - REAL(KIND=R8) ALAMDA, ALPHA, ALSQ, AMAX, BLOWUP, BNORM, + DOUBLE PRECISION ALAMDA, ALPHA, ALSQ, AMAX, BLOWUP, BNORM, * DOPE(3), DRELPR, EANORM, FAC, SM, SPARAM(5), T, TAU, WMAX, Z2, * ZZ INTEGER I, IDOPE(3), IMAX, ISOL, ITEMP, ITER, ITMAX, IWMAX, J, @@ -4177,7 +4169,7 @@ SUBROUTINE SLATEC_DROTM (N, DX, INCX, DY, INCY, DPARAM) C***PURPOSE Apply a modified Givens transformation. C***LIBRARY SLATEC (BLAS) C***CATEGORY D1A8 -C***TYPE REAL(KIND=R8) (SROTM-S, DROTM-D) +C***TYPE DOUBLE PRECISION (SROTM-S, DROTM-D) C***KEYWORDS BLAS, LINEAR ALGEBRA, MODIFIED GIVENS ROTATION, VECTOR C***AUTHOR Lawson, C. L., (JPL) C Hanson, R. J., (SNLA) @@ -4231,9 +4223,8 @@ SUBROUTINE SLATEC_DROTM (N, DX, INCX, DY, INCY, DPARAM) C 920501 Reformatted the REFERENCES section. (WRB) C 180613 Renamed SLATEC_DROTM to avoid BLAS naming conflict. (THC) C***END PROLOGUE SLATEC_DROTM - USE REAL_PRECISION - REAL(KIND=R8) DFLAG, DH12, DH22, DX, TWO, Z, DH11, DH21, + DOUBLE PRECISION DFLAG, DH12, DH22, DX, TWO, Z, DH11, DH21, 1 DPARAM, DY, W, ZERO DIMENSION DX(*), DY(*), DPARAM(5) SAVE ZERO, TWO @@ -4346,7 +4337,7 @@ SUBROUTINE SLATEC_DROTMG (DD1, DD2, DX1, DY1, DPARAM) C***PURPOSE Construct a modified Givens transformation. C***LIBRARY SLATEC (BLAS) C***CATEGORY D1B10 -C***TYPE REAL(KIND=R8) (SROTMG-S, DROTMG-D) +C***TYPE DOUBLE PRECISION (SROTMG-S, DROTMG-D) C***KEYWORDS BLAS, LINEAR ALGEBRA, MODIFIED GIVENS ROTATION, VECTOR C***AUTHOR Lawson, C. L., (JPL) C Hanson, R. J., (SNLA) @@ -4400,9 +4391,8 @@ SUBROUTINE SLATEC_DROTMG (DD1, DD2, DX1, DY1, DPARAM) C 920501 Reformatted the REFERENCES section. (WRB) C 180613 Renamed SLATEC_DROTMG to avoid BLAS naming conflict. (THC) C***END PROLOGUE SLATEC_DROTMG - USE REAL_PRECISION - REAL(KIND=R8) GAM, ONE, RGAMSQ, DD1, DD2, DH11, DH12, DH21, + DOUBLE PRECISION GAM, ONE, RGAMSQ, DD1, DD2, DH11, DH12, DH21, 1 DH22, DPARAM, DP1, DP2, DQ1, DQ2, DU, DY1, ZERO, 2 GAMSQ, DFLAG, DTEMP, DX1, TWO DIMENSION DPARAM(5) @@ -4579,7 +4569,7 @@ SUBROUTINE DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, C***SUBSIDIARY C***PURPOSE Subsidiary to DWNNLS C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLIT-S, DWNLIT-D) +C***TYPE DOUBLE PRECISION (WNLIT-S, DWNLIT-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -4605,10 +4595,9 @@ SUBROUTINE DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, C 900328 Added TYPE section. (WRB) C 900604 DP version created from SP version. . (RWC) C***END PROLOGUE DWNLIT - USE REAL_PRECISION INTEGER IDOPE(*), IPIVOT(*), ITYPE(*), L, M, MDW, N - REAL(KIND=R8) DOPE(*), H(*), RNORM, SCALE(*), W(MDW,*) + DOUBLE PRECISION DOPE(*), H(*), RNORM, SCALE(*), W(MDW,*) LOGICAL DONE C EXTERNAL DCOPY, DH12, SLATEC_DROTM, SLATEC_DROTMG, DSCAL, DSWAP, @@ -4616,7 +4605,7 @@ SUBROUTINE DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, INTEGER IDAMAX LOGICAL DWNLT2 C - REAL(KIND=R8) ALSQ, AMAX, EANORM, FACTOR, HBAR, RN, SPARAM(5), + DOUBLE PRECISION ALSQ, AMAX, EANORM, FACTOR, HBAR, RN, SPARAM(5), * T, TAU INTEGER I, I1, IMAX, IR, J, J1, JJ, JP, KRANK, L1, LB, LEND, ME, * MEND, NIV, NSOLN @@ -4869,7 +4858,7 @@ SUBROUTINE DWNLT1 (I, LEND, MEND, IR, MDW, RECALC, IMAX, HBAR, H, C***SUBSIDIARY C***PURPOSE Subsidiary to WNLIT C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLT1-S, DWNLT1-D) +C***TYPE DOUBLE PRECISION (WNLT1-S, DWNLT1-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -4885,10 +4874,9 @@ SUBROUTINE DWNLT1 (I, LEND, MEND, IR, MDW, RECALC, IMAX, HBAR, H, C 890620 Code extracted from WNLIT and made a subroutine. (RWC)) C 900604 DP version created from SP version. (RWC) C***END PROLOGUE DWNLT1 - USE REAL_PRECISION INTEGER I, IMAX, IR, LEND, MDW, MEND - REAL(KIND=R8) H(*), HBAR, SCALE(*), W(MDW,*) + DOUBLE PRECISION H(*), HBAR, SCALE(*), W(MDW,*) LOGICAL RECALC C EXTERNAL IDAMAX @@ -4934,7 +4922,7 @@ LOGICAL FUNCTION DWNLT2 (ME, MEND, IR, FACTOR, TAU, SCALE, WIC) C***SUBSIDIARY C***PURPOSE Subsidiary to WNLIT C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLT2-S, DWNLT2-D) +C***TYPE DOUBLE PRECISION (WNLT2-S, DWNLT2-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -4964,12 +4952,11 @@ LOGICAL FUNCTION DWNLT2 (ME, MEND, IR, FACTOR, TAU, SCALE, WIC) C 890620 Code extracted from WNLIT and made a subroutine. (RWC)) C 900604 DP version created from SP version. (RWC) C***END PROLOGUE DWNLT2 - USE REAL_PRECISION - REAL(KIND=R8) FACTOR, SCALE(*), TAU, WIC(*) + DOUBLE PRECISION FACTOR, SCALE(*), TAU, WIC(*) INTEGER IR, ME, MEND C - REAL(KIND=R8) RN, SN, T + DOUBLE PRECISION RN, SN, T INTEGER J C C***FIRST EXECUTABLE STATEMENT DWNLT2 @@ -4995,7 +4982,7 @@ SUBROUTINE DWNLT3 (I, IMAX, M, MDW, IPIVOT, H, W) C***SUBSIDIARY C***PURPOSE Subsidiary to WNLIT C***LIBRARY SLATEC -C***TYPE REAL(KIND=R8) (WNLT3-S, DWNLT3-D) +C***TYPE DOUBLE PRECISION (WNLT3-S, DWNLT3-D) C***AUTHOR Hanson, R. J., (SNLA) C Haskell, K. H., (SNLA) C***DESCRIPTION @@ -5011,14 +4998,13 @@ SUBROUTINE DWNLT3 (I, IMAX, M, MDW, IPIVOT, H, W) C 890620 Code extracted from WNLIT and made a subroutine. (RWC)) C 900604 DP version created from SP version. (RWC) C***END PROLOGUE DWNLT3 - USE REAL_PRECISION INTEGER I, IMAX, IPIVOT(*), M, MDW - REAL(KIND=R8) H(*), W(MDW,*) + DOUBLE PRECISION H(*), W(MDW,*) C EXTERNAL DSWAP C - REAL(KIND=R8) T + DOUBLE PRECISION T INTEGER ITEMP C C***FIRST EXECUTABLE STATEMENT DWNLT3