C API reference (kriging_capi.F90, kriging_st_capi.f90)#

The C API is an ISO C Binding wrapper around the spatial t_kriging and space-time t_kriging_st Fortran types. Both inherit from t_kriging_base; see Fortran architecture for the shared solve framework.

Every public method is exposed as a C-callable function that takes an opaque 64-bit integer handle instead of the Fortran derived type. Handles are registry slot indices, not raw pointers.

All functions return integer(c_int) ierr: 0 = success, non-zero = error. Retrieve the error message with krige_get_last_error.

Design conventions#

Convention

Detail

Handle type

int64 (shared polymorphic registry slot index, not a raw pointer)

Boolean flags

int32: 0 = false, 1 = true

Strings

Null-terminated C char arrays; converted internally with c2fstr

Array layout

Fortran column-major (ndim, nobs) — Python transposes before calling

Array sizes

Passed explicitly alongside every pointer; no assumed-size except pointweight in krige_set_grid_block

Optional args

Python always supplies concrete values; no sentinel / has_* logic needed here


Fortran object model#

The spatial and ST CAPIs share the same handle infrastructure in kriging_capi_common.F90. That registry stores class(t_kriging_base) pointers, so a slot can hold either concrete type:

class(t_kriging_base)
  +-- type(t_kriging)       spatial API, kriging_capi.F90
  +-- type(t_kriging_st)    ST API, kriging_st_capi.f90

Each API module retrieves the base pointer and then downcasts it with select type before calling concrete-only methods. Shared methods, including persistent-factor accessors, live on t_kriging_base and are available to both spatial and ST objects.

Lifecycle#

krige_create#

int krige_create(int64_t *handle)

Allocates a new t_kriging object and returns its registry slot as an opaque handle. Python stores this handle and passes it to every subsequent call.

krige_destroy#

int krige_destroy(int64_t *handle)

Finalises and deallocates the object; zeros the handle so stale use is caught early.


Initialization#

krige_initialize#

int krige_initialize(int64_t handle,
    int ndim, int nvar, int ndrift,
    int unbias, int nsim,
    int anisotropic_search, int weight_correction,
    int use_old_weight, int store_weight,
    int cross_validation, int write_mat,
    int neglect_error, int varying_vgm, int std_ck, int verbose,
    int pf_cache,
    const char *weight_file,
    const double bounds[2],
    int seed)

Must be called once after krige_create and before any other function.

Parameter

Type

Description

ndim

int

Spatial dimensions: 2 or 3

nvar

int

Number of variables (1 = kriging, >1 = co-kriging)

ndrift

int

Number of external drift functions (0 = none)

unbias

int

1 = ordinary kriging; 0 = simple kriging

nsim

int

0 = kriging only; >0 = number of SGSIM realisations

anisotropic_search

int

0/1 use anisotropic search ellipse

weight_correction

int

0/1 clip negative weights and re-normalise

use_old_weight

int

0/1 read weights from weight_file

store_weight

int

0/1 write weights to weight_file

cross_validation

int

0/1 leave-one-out cross-validation mode

write_mat

int

0/1 write kriging matrix to CSV (debug)

neglect_error

int

0/1 write NaN instead of stopping on singular matrix

varying_vgm

int

0/1 use a different variogram per block (SVA mode)

std_ck

int

0/1 co-kriging unbiasedness formulation (only used when nvar > 1 and unbias = 1). 1 = standard co-kriging: separate per-variable constraints (Σw₁=1, Σw₂=0); matches gstat/ISATIS. 0 = Isaaks & Srivastava: single combined constraint (Σw₁+Σw₂=1) plus local-mean correction. Default: 1.

verbose

int

0/1 print progress messages

pf_cache

int

0/1 enable the persistent between-solve factorisation cache. When 1, the Cholesky factorisation of K is stored after the first krige_solve and reused on subsequent calls when observations and variogram are unchanged. Default: 0 (disabled).

weight_file

char*

Path for weight file; empty string when unused

bounds

double[2]

[lower, upper] clipping bounds on the estimate

seed

int

Random seed for SGSIM (0 = use clock)


Observations#

krige_set_obs#

int krige_set_obs(int64_t handle,
    int ivar, int nobs, int ndim_c,
    const double *coord,    // [ndim_c × nobs], Fortran column-major
    const double *value,    // [nobs]
    const double *variance, // [nobs]  measurement error variance; 0 if unknown
    double sk_mean)

Sets coordinates, values, and per-observation measurement variance for one variable. sk_mean is only used when unbias=0. Configure nmax and maxdist with krige_set_search. Duplicate coordinate tuples within the same variable are rejected. For spatial kriging, all ndim_c coordinate rows must match to count as a duplicate.

krige_set_obs_drift#

int krige_set_obs_drift(int64_t handle,
    int ivar, int ndrift_c, int nobs,
    const double *drift)    // [ndrift_c × nobs], Fortran column-major

Sets external drift values at observation locations. Must be called after krige_set_obs for the same ivar, and only when ndrift > 0.

krige_update_obs_value#

int krige_update_obs_value(int64_t handle,
    int ivar, int nobs,
    const double *value)    // [nobs]

Replaces observation values in-place without touching coordinates or the KD-tree. Use with use_old_weight to re-estimate with new data without recomputing search neighbourhoods or the LHS factorisation.

Score transforms#

int krige_transform_value_to_score(int64_t handle,
    int ivar, int n,
    const double *value,  // [n] data units
    double *score)        // [n] normal or uniform scores

int krige_transform_score_to_value(int64_t handle,
    int ivar, int n,
    const double *score,  // [n] normal or uniform scores
    double *value)        // [n] data units

Apply the active marginal transform created by krige_set_nscore or krige_set_uscore. The score type follows the active transform for ivar: normal scores for set_nscore, uniform CDF scores for set_uscore.


Variogram#

krige_reset_vgm#

int krige_reset_vgm(int64_t handle, int ivar, int jvar)

Clears all nested structures for the (ivar, jvar) pair (and its mirror (jvar, ivar) for cross-variograms) across every block. Call before krige_set_vgm when reusing an object with a different variogram model.

krige_set_vgm#

int krige_set_vgm(int64_t handle,
    int ivar, int jvar,
    const char *vtype,      // null-terminated: "sph" "exp" "gau" "pow" "lin" "hol" "bsq" "cir" "nug"
    double nugget, double sill,
    double a_major, double a_minor1, double a_minor2,
    double azimuth, double dip, double plunge)

Appends one nested structure for the (ivar, jvar) variable pair. Call multiple times to build a composite model. For co-kriging the LMC constraint b12² b11 × b22 must hold per structure.

krige_set_vgm_block#

int krige_set_vgm_block(int64_t handle,
    int ivar, int jvar, int ib,
    const char *vtype,
    double nugget, double sill,
    double a_major, double a_minor1, double a_minor2,
    double azimuth, double dip, double plunge)

Same as krige_set_vgm but targets a single block ib (1-based). Requires varying_vgm=1 in krige_initialize and krige_set_grid to have been called first (block count must be known).


Grid#

krige_set_grid#

int krige_set_grid(int64_t handle,
    int ngrid, int ndim_c,
    const double *coord,        // [ndim_c × ngrid], Fortran column-major
    const double *rangescale,   // [ngrid]  variogram range scaling; pass 1.0 when unused
    const double *localnugget)  // [ngrid]  additional per-block nugget; pass 0.0 when unused

Sets estimation targets for point kriging. Use krige_set_grid_block for block kriging or krige_set_grid_cv for cross-validation.

krige_set_grid_block#

int krige_set_grid_block(int64_t handle,
    int block_type,
    int ngrid, int ndim_c,
    const double *coord,        // [ndim_c × ngrid]
    int nblock,
    const int *nblockpnt,       // [nblock]  sub-nodes per block
    const double *pointweight,  // [sum(nblockpnt)]
    const double *blocksize,    // [ndim_c × nblock]
    const double *rangescale,   // [nblock]
    const double *localnugget)  // [nblock]

Sets estimation targets for block kriging.

block_type

Meaning

-4

Gaussian quadrature sub-nodes (auto-generated)

> 0

User-supplied sub-nodes

pointweight length equals sum(nblockpnt); Fortran derives it via sum(nblockpnt) so no separate npw argument is needed.

krige_set_grid_cv#

int krige_set_grid_cv(int64_t handle)

Configures leave-one-out cross-validation mode. Fortran derives the grid from the observation coordinates automatically. Call instead of krige_set_grid when cross_validation=1.

krige_set_grid_drift#

int krige_set_grid_drift(int64_t handle,
    int ivar,
    int ndrift_c, int nblocks,
    const double *drift)    // [ndrift_c × nblocks], Fortran column-major

Sets external drift values at block locations. Must be called after any krige_set_grid* variant, and only when ndrift > 0.

Parameter

Description

ivar

Target-variable index (1-based) whose RHS receives this drift. Pass ivar < 0 to broadcast the same drift to all target variables — the common case when external drift is the same regardless of which variable is being estimated.

ndrift_c

Number of drift functions (= ndrift from krige_initialize)

nblocks

Number of blocks (block%n, not grid%n for block kriging)

drift

Drift values [ndrift_c × nblocks], Fortran column-major

Note on ivar semantics. In krige_set_obs_drift, ivar identifies the source variable (whose observations form the F-matrix column). Here it identifies the target variable (which estimation’s RHS uses this drift). These are opposite ends of the kriging system — both use 1-based indexing and < 0 is not valid for krige_set_obs_drift.


Simulation (SGSIM)#

krige_set_sim#

int krige_set_sim(int64_t handle,
    int nblocks,
    const int    *randpath,  // [nblocks]  random visiting order
    int nsim_c, int nvar_c,
    const double *sample)    // [nsim_c × nvar_c × nblocks]

Supplies the random visiting path and pre-drawn standard-normal samples for SGSIM. Python generates both arrays before calling. Call after krige_set_grid and before krige_set_search. Only needed when nsim > 0.



Solve#

krige_prepare#

int krige_prepare(int64_t handle)

Pre-allocates result arrays and sets up the block loop. Called automatically by krige_solve; exposed separately for benchmarking.

krige_solve#

int krige_solve(int64_t handle, int nthread, int ncache)

Runs the kriging or SGSIM block loop. nthread = 0 uses the OpenMP runtime default; nthread > 0 caps the thread count for this call only. ncache = -1 keeps the object’s current hcache slot default, ncache = 0 disables the factor caches for this solve, and positive ncache sets the total shared hcache pool size for this solve only. Values from 1 through 3 are promoted to 4 because one set-associative bucket requires four slots. Results are available via the getters immediately after this returns.


Result getters#

krige_get_nblocks#

int krige_get_nblocks(int64_t handle, int *n)

Returns the number of estimation blocks.

krige_get_nsim#

int krige_get_nsim(int64_t handle, int *n)

Returns the number of simulations (1 for plain kriging).

krige_get_block_coord#

int krige_get_block_coord(int64_t handle,
    int ndim_c, int nblocks,
    double *out)    // [ndim_c × nblocks], Fortran column-major

Copies the block centroid coordinates into a caller-allocated buffer. The output is filled as out(ndim_c, nblocks) in Fortran column-major order. Python allocates a (ndim_c, nblocks) Fortran-order array and transposes to (nblocks, ndim_c) for standard row-major convention.

For SGSIM, blocks are reordered back to the original (non-randomised) index order inside krige_solve, so the coordinates returned here always correspond to krige_get_estimate_all positions at the same block index.

krige_get_estimate#

int krige_get_estimate(int64_t handle,
    int nsim_c, int nblocks,
    double *out)    // [nblocks × nsim_c], C row-major

Copies the primary-variable estimate for all simulations. Output is out[ib, isim] for block ib (0-based) and simulation isim (0-based).

krige_get_estimate_all#

int krige_get_estimate_all(int64_t handle,
    int nblocks, int nvar_c, int nsim_c,
    double *out)    // [nblocks × nvar_c × nsim_c]

Copies estimates for all variables and simulations. Output convention: out[ib, kvar, isim] — block index first, matching the (nobs, ndim) coord convention.

krige_get_variance#

int krige_get_variance(int64_t handle,
    int nblocks,
    double *out)    // [nblocks]

Copies the primary-variable kriging variance.

krige_get_variance_all#

int krige_get_variance_all(int64_t handle,
    int nblocks, int nvar_c,
    double *out)    // [nblocks × nvar_c × nvar_c]

Copies the full conditional covariance matrix at every block. out[ib, iv, jv] is the covariance between variables iv+1 and jv+1 at block ib+1; the diagonal out[ib, k, k] is variable k+1’s kriging variance.


Persistent factorisation cache#

Between-solve factorisation caching allows the Cholesky factorisation of the kriging covariance matrix K to be reused across successive krige_solve calls when observation coordinates and the variogram have not changed.

The cache is populated automatically after the first successful krige_solve and invalidated by krige_set_obs or krige_set_vgm. These two functions let Python query the cached matrices for inspection or debugging.

Internally there are three cache layers:

  • ctx%cache: one single-entry cache per worker thread. It handles immediate repeated systems and is pre-warmed from the persistent cache when available.

  • self%hcache: one bounded, shared, four-way set-associative cache for repeated systems within a single solve. Each bucket has its own lock and local LRU order, so unrelated neighbourhoods remain concurrent. Every hash candidate is verified against the full neighbour key with fcache_matches, so collisions are safe.

  • self%pf: the optional persistent cache shared by the kriging object across solves when pf_cache=1. This cache also stores the assembled matA/rhsB snapshot exposed by krige_get_factor_system.

For self%hcache, the effective pool is divided into floor(ncache / 4) buckets of four ways. Positive values below four are promoted to four; other non-multiples are rounded down to a complete bucket. A hash selects one bucket, and only its four entries are scanned. Empty ways are filled first; otherwise the bucket-local least-recently-used entry is replaced.

Lookup takes the selected bucket’s lock and copies a verified hit into the worker’s single-entry ctx%cache. Insertion uses omp_test_lock and skips publication if that bucket is contended. This avoids making completed factorizations wait merely to improve future cache hit rate. Skipped inserts and duplicate concurrent factorizations do not affect numerical results. See the architecture guide for the data layout, key construction, and synchronization details.

self%pf is saved after the parallel block loop from a thread whose current ctx%matA/ctx%rhsB still match its thread-local factors. hcache hits do not update the assembled system, so they are not used to populate matA/rhsB for inspection.

The multi-slot cache size is controlled in Fortran by factor_cache_size (default 256 total shared slots) and by a total shared-pool byte cap (MAX_HCACHE_BYTES). The byte cap accounts for the Cholesky factors, neighbour keys, and worst-case lazy SSYTRF fallback storage.

For cache-path testing, the Python solve(ncache=...) argument and C API krige_solve(..., ncache) argument override this slot count for one solve call. Use ncache=0 to disable all factorization reuse, ncache=4 for the smallest shared pool, or ncache=None in Python / ncache=-1 in C to keep the compiled default. The Makefile HCACHE variable sets that compiled default: make HCACHE=0 disables factorization reuse by default, make HCACHE=4 builds the smallest shared pool, and bare make uses the normal 256-slot default. Positive values from 1 through 3 are promoted to 4. The optional persistent self%pf cache is configured separately.

krige_get_factor_info#

int krige_get_factor_info(int64_t handle,
    int *npp_out,    // number of neighbours (rows/cols of K)
    int *p_out,      // drift + unbiasedness columns (Schur size)
    int *valid_out)  // 1 = valid; 0 = not yet computed or invalidated

Returns the dimensions and validity of the cached factorisation. Call this first to obtain npp and p before allocating buffers for krige_get_factor_matrices.

krige_get_factor_matrices#

int krige_get_factor_matrices(int64_t handle,
    int npp, int p,
    double *L_out,       // [npp * npp]       upper-tri Cholesky of K
    double *kinv_out,    // [npp * max(1,p)]   K^{-1} F
    double *schur_out)   // [max(1,p) * max(1,p)]  Cholesky of F'K^{-1}F

Copies the three persistent factor matrices into caller-allocated arrays. npp and p must match the values returned by krige_get_factor_info.

All arrays are in Fortran column-major order. The solver uses uplo='U' (LAPACK spotrf), so:

  • L_out — upper triangle is the Cholesky factor U; K = U' U. The lower triangle retains the original K values and should be ignored.

  • kinv_outK^{-1} F where F is the full drift matrix.

  • schur_out — upper triangle is the Cholesky factor of F' K^{-1} F.

krige_get_factor_system#

int krige_get_factor_system(int64_t handle,
    int npp, int p, int nvar,
    double *matA_out,    // [npp+p x npp+p] assembled LHS before factorization
    double *rhsB_out)    // [nvar x npp+p] assembled RHS before solving

Copies the assembled linear system that produced the persistent factor. This is for inspection and debugging; npp, p, and nvar must match the kriging object and the dimensions returned by krige_get_factor_info.

Invalidation rules#

Trigger

Effect on persistent cache

krige_set_obs

pf_valid = false (coordinates change K)

krige_set_vgm

pf_valid = false (variogram changes K)

krige_update_obs_value

No effect — values do not enter K

Python wrapper#

f = kriging_object.get_factor()
# f['valid']       bool
# f['npp']         int — size of K
# f['p']           int — size of Schur complement
# f['L']           ndarray (npp, npp)       — upper-triangular factor
# f['kinv_drift']  ndarray (npp, max(1,p))  — K^{-1} F
# f['schur']       ndarray (max(1,p), max(1,p))
# f['matA']        ndarray (npp+p, npp+p)   — assembled LHS
# f['rhsB']        ndarray (nvar, npp+p)    — assembled RHS

Weight store#

krige_free_weight_store#

int krige_free_weight_store(int64_t handle)

Frees the in-memory weight store.

krige_get_weight_nnear#

int krige_get_weight_nnear(int64_t handle,
    int ngroups_c, int nblock_c,
    int *out)    // [ngroups_c × nblock_c]

Copies the neighbour-count array. ngroups = nvar for kriging, ngroups = 2 * nvar for SGSIM.

krige_get_weight_inear#

int krige_get_weight_inear(int64_t handle,
    int nmax_c, int ngroups_c, int nblock_c,
    int *out)    // [nmax_c × ngroups_c × nblock_c]

Copies the 1-based observation-index array for every neighbour slot. Padded slots are zero.

krige_get_weight_data#

int krige_get_weight_data(int64_t handle,
    int nmax_c, int ngroups_c, int nvar_c, int nblock_c,
    double *out)    // [nmax_c × ngroups_c × nvar_c × nblock_c]

Copies kriging weights. Padded slots are zero.

krige_get_weight_var#

int krige_get_weight_var(int64_t handle,
    int nvar_c, int nblock_c,
    double *out)    // [nvar_c × nvar_c × nblock_c]

Copies the stored kriging variances (shape matches krige_get_variance_all transposed to Fortran order).

krige_set_weights#

int krige_set_weights(int64_t handle,
    int nmax_c, int ngroups_c, int nvar_c, int nblock_c,
    const int    *nnear_in,   // [ngroups_c × nblock_c]
    const int    *inear_in,   // [nmax_c × ngroups_c × nblock_c]
    const double *weight_in,  // [nmax_c × ngroups_c × nvar_c × nblock_c]
    const int    *order_in,   // [nblock_c]
    const double *var_in)     // [nvar_c × nvar_c × nblock_c]

Loads pre-computed weights into the store and sets use_old_weight = true so that the next krige_solve applies them directly without solving the kriging system.


Utilities#

krige_get_last_error#

int krige_get_last_error(char *buffer, int nbuf)

Copies the last error message (null-terminated) into buffer. nbuf is the buffer capacity in bytes. Returns 0 even when an error exists — the error string itself carries the information.

krige_to_str#

int64_t krige_to_str(int64_t handle)

Returns a pointer to a null-terminated Fortran character array containing a human-readable summary of the kriging object. Returns 0 on error. Do not free the returned pointer — it is owned by the Fortran object.

krige_get_max_threads / krige_get_num_threads#

void krige_get_max_threads(int *n)
void krige_get_num_threads(int *n)

Query the OpenMP thread count. Both return 1 when the library is compiled without OpenMP (--no-openmp).


Space-time API (kriging_st_capi.f90)#

All ST entry points are prefixed krige_st_ and share the same handle registry as the spatial API. Differences from the spatial API are noted below; methods not listed here have the same signature as their krige_ counterparts (e.g. krige_st_solve, krige_st_prepare, krige_st_get_estimate, …).

krige_st_set_obs#

As with krige_set_obs, duplicate observation coordinate tuples are rejected. For ST data, the duplicate key includes spatial coordinates and time.

int krige_st_set_obs(int64_t handle,
    int ivar, int nobs,
    const double *coord,    // [4 × nobs], Fortran column-major; rows 1:3 = spatial, row 4 = time
    const double *value,    // [nobs]
    const double *variance, // [nobs]
    double sk_mean)

Identical to krige_set_obs except that coord has 4 rows (3 spatial + 1 time). Configure nmax and maxdist with krige_st_set_search.

krige_st_set_grid#

int krige_st_set_grid(int64_t handle,
    int ngrid,
    const double *coord,  // [3 × ngrid], Fortran column-major (spatial only)
    const double *time)   // [ngrid]

Accepts spatial coordinates and times as separate arrays (unlike krige_set_grid which takes a single ndim-row coord array).

krige_st_set_st_model#

int krige_st_set_st_model(int64_t handle,
    const char *model,       // "sum_metric" or "product_sum"
    const char *transform,   // variogram type for f_time: "lin", "exp", "sph", …
    double at,               // joint temporal scale (same units as time coordinate)
    double time_nugget,      // temporal variogram nugget
    double time_sill,        // temporal variogram sill
    double k_ps)             // product-sum k (ignored for sum_metric)

Sets global ST model parameters shared by all variogram entries.

krige_st_set_vgm#

int krige_st_set_vgm(int64_t handle,
    int ivar,
    int jvar,
    const char *vtype,
    double nugget,
    double sill,
    double a_major,
    double a_minor1,
    double a_minor2,
    double azimuth,
    double dip,
    double plunge,
    int is_product) // 0 = additive, 1 = multiply the preceding spatial structure

Adds one spatial marginal structure. Consecutive structures with is_product=1 are multiplied in covariance space, matching krige_set_vgm.

krige_st_set_vgm_temporal#

int krige_st_set_vgm_temporal(int64_t handle,
    int ivar,
    int jvar,
    const char *vtype,
    double nugget,
    double sill,
    double at_k,
    int is_product) // 0 = additive, 1 = multiply the preceding temporal structure

Adds one temporal marginal structure. at_k is in the same units as the observation time coordinate. Consecutive structures with is_product=1 are multiplied in covariance space, matching krige_set_vgm.


Handle registry (internals)#

Python receives a 1-based slot index (int64) rather than a raw Fortran pointer. This avoids passing non-C-interoperable derived-type addresses through ctypes.

Function

Description

store_obj

Finds the first free slot; grows the registry (doubles) if full

get_obj

Validates the slot index and recovers the typed Fortran pointer

release_obj

Nullifies the slot on krige_destroy; does not compact the array

The registry starts at 16 slots and doubles on demand, keeping old slot numbers stable for existing Python handles.