rnatm

rnatm_voxel_partition.c (17384B)

---

 /* Copyright (C) 2022, 2023, 2025 Centre National de la Recherche Scientifique
  * Copyright (C) 2022, 2023, 2025 Institut Pierre-Simon Laplace
  * Copyright (C) 2022, 2023, 2025 Institut de Physique du Globe de Paris
  * Copyright (C) 2022, 2023, 2025 |Méso|Star> (contact@meso-star.com)
  * Copyright (C) 2022, 2023, 2025 Observatoire de Paris
  * Copyright (C) 2022, 2023, 2025 Université de Reims Champagne-Ardenne
  * Copyright (C) 2022, 2023, 2025 Université de Versaille Saint-Quentin
  * Copyright (C) 2022, 2023, 2025 Université Paul Sabatier
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #include "rnatm_c.h"
 #include "rnatm_log.h"
 #include "rnatm_voxel_partition.h"
 
 #include <rsys/condition.h>
 #include <rsys/mem_allocator.h>
 #include <rsys/mutex.h>
 #include <rsys/ref_count.h>
 
 /* A tile stores N*M voxels where N is the number of voxels and M is the voxel
  * width as defined when creating the pool. This width makes it possible to
  * store in the same tile the radiative coefficients of a voxel for several
  * spectral data. For each voxel, the tile actually stores M voxels stored one
  * after the other */
  struct tile {
   struct list_node node;
   float voxels[1]; /* Flexible array member */
 };
 
 struct partition {
   struct list_node node;
   size_t id; /* Unique identifier of the partition */
   struct tile* tile; /* Set of voxels */
   struct pool* pool;
 
   /* Number of references on the partition. Its initial value is set when the
    * partition is commited. Once fetched, a reference is given to the caller
    * that it then released by calling partition_free */
   ATOMIC ref;
 };
 
 struct pool {
   /* Allocated partitions and tiles */
   struct partition* parts;
   char* tiles;
   size_t tile_sz; /* Size in bytes of a tile */
 
   struct list_node parts_free; /* List of free partitions */
   /* List of committed partition sorted in ascending order wrt partition id */
   struct list_node parts_commit;
 
   struct list_node tiles_free; /* List of available tiles of voxels */
 
   /* Dummy voxel returned by a partition when no voxel is reserved for it */
   float empty_voxel[NFLOATS_PER_VOXEL];
 
   struct mutex* mutex;
   struct cond* cond_new;
   struct cond* cond_fetch;
   struct cond* cond_tile;
 
   size_t next_part_id; /* Identifier of the next partition */
   size_t partition_definition; /* #voxels along the 3 axis */
   size_t partition_nvoxels; /* Overall number of voxels in a partition */
   size_t voxel_width; /* Number of items per voxel data */
 
   struct mem_allocator* allocator;
   ATOMIC error; /* Is the pool not valid? */
   ref_T ref;
 };
 
 /*******************************************************************************
  * Helper functions
  ******************************************************************************/
 static INLINE res_T
 check_pool_create_args(const struct pool_create_args* args)
 {
   if(!args
   || !args->npartitions
   || !args->voxel_width
   || args->npartitions < args->npreallocated_partitions
   || !IS_POW2(args->partition_definition)) {
     return RES_BAD_ARG;
   }
 
   return RES_OK;
 }
 
 static INLINE void
 tile_init(struct tile* tile)
 {
   ASSERT(tile);
   list_init(&tile->node);
 }
 
 static INLINE void
 tile_release(struct tile* tile)
 {
   ASSERT(tile && is_list_empty(&tile->node));
   (void)tile;
 }
 
 static void
 partition_init
   (struct pool* pool,
    struct partition* partition)
 {
   ASSERT(pool && partition);
   list_init(&partition->node);
   partition->id = SIZE_MAX;
   partition->tile = NULL;
   partition->pool = pool;
   partition->ref = 0;
 }
 
 static void
 partition_release(struct partition* partition)
 {
   ASSERT(partition && is_list_empty(&partition->node));
   if(partition->tile) {
     ASSERT(is_list_empty(&partition->tile->node));
     list_add(&partition->pool->tiles_free, &partition->tile->node);
   }
 }
 
 static struct partition*
 get_free_partition(struct pool* pool)
 {
   struct list_node* node = NULL;
   struct partition* partition = NULL;
   ASSERT(pool);
 
   mutex_lock(pool->mutex);
 
   /* Waits for a free partition */
   while(is_list_empty(&pool->parts_free) && !pool->error) {
     cond_wait(pool->cond_new, pool->mutex);
   }
 
   if(pool->error) {
     /* An error occurs */
     partition = NULL;
     mutex_unlock(pool->mutex);
 
   } else {
     /* Retrieve the next available partition */
     node = list_head(&pool->parts_free);
     list_del(node);
     mutex_unlock(pool->mutex);
 
     partition = CONTAINER_OF(node, struct partition, node);
     partition->id = SIZE_MAX;
   }
   return partition;
 }
 
 static res_T
 reserve_partition_voxels(struct partition* partition)
 {
   struct pool* pool = NULL;
   res_T res = RES_OK;
   ASSERT(partition && !partition->tile);
 
   pool = partition->pool;
 
   mutex_lock(partition->pool->mutex);
 
   /* Waits for a free tile */
   while(is_list_empty(&pool->tiles_free) && !pool->error) {
     cond_wait(pool->cond_tile, pool->mutex);
   }
 
   if(pool->error) {
     /* An error occurs */
     mutex_unlock(pool->mutex);
     res = RES_UNKNOWN_ERR;
 
   } else {
     struct list_node* node = NULL;
 
     /* Get an available voxel tile */
     node = list_head(&pool->tiles_free);
     list_del(node);
     mutex_unlock(pool->mutex);
 
     partition->tile = CONTAINER_OF(node, struct tile, node);
   }
   return res;
 }
 
 static void
 release_pool(ref_T* ref)
 {
   struct pool* pool = CONTAINER_OF(ref, struct pool, ref);
   struct list_node* node = NULL;
   struct list_node* tmp_node = NULL;
   ASSERT(ref);
 
   if(pool->mutex) mutex_destroy(pool->mutex);
   if(pool->cond_new) cond_destroy(pool->cond_new);
   if(pool->cond_fetch) cond_destroy(pool->cond_fetch);
   if(pool->cond_tile) cond_destroy(pool->cond_tile);
 
   LIST_FOR_EACH_SAFE(node, tmp_node, &pool->parts_free) {
     struct partition* partition = CONTAINER_OF(node, struct partition, node);
     list_del(node);
     partition_release(partition);
   }
 
   LIST_FOR_EACH_SAFE(node, tmp_node, &pool->parts_commit) {
     struct partition* partition = CONTAINER_OF(node, struct partition, node);
     list_del(node);
     partition_release(partition);
   }
 
   LIST_FOR_EACH_SAFE(node, tmp_node, &pool->tiles_free) {
     struct tile* tile = CONTAINER_OF(node, struct tile, node);
     list_del(node);
     tile_release(tile);
   }
 
   ASSERT(is_list_empty(&pool->parts_free));
   ASSERT(is_list_empty(&pool->parts_commit));
   ASSERT(is_list_empty(&pool->tiles_free));
 
   MEM_RM(pool->allocator, pool->parts);
   MEM_RM(pool->allocator, pool->tiles);
   MEM_RM(pool->allocator, pool);
 }
 
 /*******************************************************************************
  * Partition of voxels
  ******************************************************************************/
 void
 partition_free(struct partition* partition)
 {
   struct pool* pool = NULL;
   int free_tile = 0;
   ASSERT(partition);
 
   pool = partition->pool;
 
   /* The partition reference counter can be zero before being decremented: it is
    * initialized to the voxel width once the partition has been commited. So, if
    * the partition is released before any commit (for example, when it doesn't
    * overlap the atmosphere grid), it can be negative after it has been
    * decremented */
   if(ATOMIC_DECR(&partition->ref) > 0)
     return; /* The partition is still referenced */
 
   mutex_lock(pool->mutex);
   list_move_tail(&partition->node, &pool->parts_free); /* Free the partition */
   partition->ref = 0; /* Reset to 0 rather than letting a negative value */
   if(partition->tile) { /* Free the reserved tile */
     list_move_tail(&partition->tile->node, &pool->tiles_free);
     partition->tile = NULL;
     free_tile = 1;
   }
   mutex_unlock(pool->mutex);
 
   /* Notify a thread waiting for a free partition that we just registered one */
   cond_signal(pool->cond_new);
 
   if(free_tile) {
     /* Notify a partition waiting for a free tile that we just registered one */
     cond_signal(pool->cond_tile);
   }
 }
 
 void
 partition_empty(struct partition* partition)
 {
   ASSERT(partition && partition->tile);
 
   mutex_lock(partition->pool->mutex);
   list_move_tail(&partition->tile->node, &partition->pool->tiles_free);
   partition->tile = NULL;
   mutex_unlock(partition->pool->mutex);
 
   /* Notify a partition waiting for a free tile that we just registered one */
   cond_signal(partition->pool->cond_tile);
 }
 
 void
 partition_commit(struct partition* partition, const size_t refs_count)
 {
   struct list_node* node = NULL;
   struct pool* pool = NULL;
   ASSERT(partition);
   ASSERT(refs_count && refs_count <= partition->pool->voxel_width);
 
   pool = partition->pool;
 
   /* Setup the number of partition references */
   partition->ref = (ATOMIC)refs_count;
 
   /* Committed partitions are sorted in ascending order of their id. We are
    * therefore looking for the partition whose id is less than the id of the
    * partition to be committed in order to add the latter in the right place */
   mutex_lock(pool->mutex);
   LIST_FOR_EACH_REVERSE(node, &pool->parts_commit) {
     struct partition* partition2 = CONTAINER_OF(node, struct partition, node);
     if(partition2->id < partition->id) break;
   }
   list_add(node, &partition->node);
   mutex_unlock(pool->mutex);
 
   /* Notify the threads waiting for a valid partition that we just registered
    * one. Note that a partition can register several items when voxel_width > 1
    * and therefore several threads can wait for the same partition.
    * Consequently, we broadcast the signal to all threads that are blocked on
    * fetch condition */
   cond_broadcast(pool->cond_fetch);
 }
 
 size_t
 partition_get_id(const struct partition* partition)
 {
   ASSERT(partition);
   return partition->id;
 }
 
 size_t
 partition_get_definition(const struct partition* partition)
 {
   ASSERT(partition);
   return partition->pool->partition_definition;
 }
 
 float*
 partition_get_voxel
   (struct partition* part,
    const size_t ivoxel,
    const size_t iitem)
 {
   ASSERT(part && ivoxel < part->pool->partition_nvoxels && part->tile != NULL);
   ASSERT(iitem < part->pool->voxel_width);
   return part->tile->voxels
     + NFLOATS_PER_VOXEL*(ivoxel*part->pool->voxel_width + iitem);
 }
 
 const float*
 partition_cget_voxel
   (struct partition* part,
    const size_t ivoxel,
    const size_t iitem)
 {
   ASSERT(part && ivoxel < part->pool->partition_nvoxels);
   ASSERT(iitem < part->pool->voxel_width);
   if(part->tile == NULL) {
     return part->pool->empty_voxel;
   } else {
     return partition_get_voxel(part, ivoxel, iitem);
   }
 }
 
 void
 partition_clear_voxels(struct partition* partition)
 {
   size_t ivoxel = 0;
   ASSERT(partition);
 
   if(partition->tile == NULL) return; /* Nothing to do */
 
   FOR_EACH(ivoxel, 0, partition->pool->partition_nvoxels) {
     size_t iitem = 0;
     FOR_EACH(iitem, 0, partition->pool->voxel_width) {
       float* voxel = partition_get_voxel(partition, ivoxel, iitem);
       voxel_clear(voxel);
     }
   }
 }
 
 void
 partition_accum(struct partition* dst, struct partition* src)
 {
   size_t ivoxel = 0;
   ASSERT(dst->pool->partition_nvoxels == src->pool->partition_nvoxels);
   ASSERT(dst->pool->voxel_width == src->pool->voxel_width);
   ASSERT(src->tile && dst->tile); /* Partition cannot be empty */
 
   FOR_EACH(ivoxel, 0, src->pool->partition_nvoxels) {
     size_t iitem = 0;
     FOR_EACH(iitem, 0, src->pool->voxel_width) {
       const float* src_voxel = partition_cget_voxel(src, ivoxel, iitem);
       float* dst_voxel = partition_get_voxel(dst, ivoxel, iitem);
       voxel_accum(dst_voxel, src_voxel);
     }
   }
 }
 
 /*******************************************************************************
  * Pool of partitions
  ******************************************************************************/
 res_T
 pool_create
   (const struct pool_create_args* args,
    struct pool** out_pool)
 {
   struct mem_allocator* allocator = NULL;
   struct pool* pool = NULL;
   size_t ipartition = 0;
   size_t nvoxels = 0;
   size_t tile_sz = 0;
 
   res_T res = RES_OK;
   ASSERT(check_pool_create_args(args) == RES_OK && out_pool);
 
   allocator = args->allocator ? args->allocator : &mem_default_allocator;
   pool = MEM_CALLOC(allocator, 1, sizeof(*pool));
   if(!pool) {
     res = RES_MEM_ERR;
     goto error;
   }
   pool->allocator = allocator;
   pool->partition_definition = args->partition_definition;
   pool->voxel_width = args->voxel_width;
   pool->partition_nvoxels =
     pool->partition_definition
   * pool->partition_definition
   * pool->partition_definition;
   ref_init(&pool->ref);
   list_init(&pool->parts_free);
   list_init(&pool->parts_commit);
   list_init(&pool->tiles_free);
   voxel_clear(pool->empty_voxel);
 
   /* Create the mutex and condition variables used to synchronize threads */
   pool->mutex = mutex_create();
   if(!pool->mutex) { res = RES_UNKNOWN_ERR; goto error; }
   pool->cond_new = cond_create();
   if(!pool->cond_new) { res = RES_UNKNOWN_ERR; goto error; }
   pool->cond_fetch = cond_create();
   if(!pool->cond_fetch) { res = RES_UNKNOWN_ERR; goto error; }
   pool->cond_tile = cond_create();
   if(!pool->cond_tile) { res = RES_UNKNOWN_ERR; goto error; }
 
   /* Allocate the partitions */
   pool->parts = MEM_CALLOC(allocator, args->npartitions, sizeof(*pool->parts));
   if(!pool->parts) { res = RES_MEM_ERR; goto error; }
 
   /* Allocate the tiles */
   nvoxels = pool->partition_nvoxels;
   tile_sz =
     sizeof(struct tile)
   - sizeof(float)/*Dummy member*/
   + sizeof(float[NFLOATS_PER_VOXEL]) * nvoxels * pool->voxel_width;
   tile_sz = ALIGN_SIZE(tile_sz, ALIGNOF(struct tile));
   pool->tiles = MEM_CALLOC(allocator, args->npreallocated_partitions, tile_sz);
   if(!pool->tiles) { res = RES_MEM_ERR; goto error; }
 
   /* Setup the free partitions */
   FOR_EACH(ipartition, 0, args->npartitions) {
     struct partition* partition = pool->parts + ipartition;
     partition_init(pool, partition);
     list_add(&pool->parts_free, &partition->node);
   }
 
   /* Setup the free tiles */
   FOR_EACH(ipartition, 0, args->npreallocated_partitions) {
     struct tile* tile = (struct tile*)(pool->tiles + ipartition*tile_sz);
     tile_init(tile);
     list_add(&pool->tiles_free, &tile->node);
   }
 
 exit:
   *out_pool = pool;
   return res;
 error:
   if(pool) { pool_ref_put(pool); pool = NULL; }
   goto exit;
 }
 
 void
 pool_ref_get(struct pool* pool)
 {
   ASSERT(pool);
   ref_get(&pool->ref);
 }
 
 void
 pool_ref_put(struct pool* pool)
 {
   ASSERT(pool);
   ref_put(&pool->ref, release_pool);
 }
 
 size_t
 pool_get_partition_definition(const struct pool* pool)
 {
   ASSERT(pool);
   return pool->partition_definition;
 }
 
 size_t
 pool_get_voxel_width(const struct pool* pool)
 {
   ASSERT(pool);
   return pool->voxel_width;
 }
 
 struct partition*
 pool_next_partition(struct pool* pool)
 {
   struct partition* partition;
   ASSERT(pool);
 
   partition = pool_dummy_partition(pool);
   if(!partition) return NULL;
 
   mutex_lock(pool->mutex);
   partition->id = pool->next_part_id;
   pool->next_part_id += 1;
   mutex_unlock(pool->mutex);
 
   return partition;
 }
 
 struct partition*
 pool_dummy_partition(struct pool* pool)
 {
   struct partition* partition;
   res_T res = RES_OK;
   ASSERT(pool);
 
   partition = get_free_partition(pool);
   res = reserve_partition_voxels(partition);
   if(res != RES_OK) {
     partition_free(partition);
     return NULL;
   }
   return partition;
 }
 
 struct partition*
 pool_fetch_partition(struct pool* pool, const size_t ipartition)
 {
   struct partition* found_partition = NULL;
   struct list_node* node = NULL;
   ASSERT(pool);
 
   mutex_lock(pool->mutex);
   while(!pool->error) {
 
     /* Search for the partition that matches the submitted id */
     LIST_FOR_EACH(node, &pool->parts_commit) {
       struct partition* partition = CONTAINER_OF(node, struct partition, node);
 
       /* The partition to fetch has been found */
       if(partition->id == ipartition) {
         found_partition = partition;
         break;
       }
 
       /* Partitions are sorted in ascending order. Stop the linear search if the
        * current id is greater than the submitted one */
       if(partition->id > ipartition) break;
     }
 
     if(!found_partition) {
       /* The partition is still not committed. The thread is waiting for it */
       cond_wait(pool->cond_fetch, pool->mutex);
 
     } else {
       /* Do not remove the found partition from the list of committed
        * partitions since the same partition can be fetched several times when
        * the width of the voxels is greater than 1 */
       /* list_del(&found_partition->node); */
       break;
     }
   }
   mutex_unlock(pool->mutex);
 
   return found_partition;
 }
 
 void
 pool_invalidate(struct pool* pool)
 {
   ASSERT(pool);
 
   /* Notifies that an error has occurred */
   mutex_lock(pool->mutex);
   pool->error = 1;
   mutex_unlock(pool->mutex);
 
   /* Wakes up all the waiting threads */
   cond_broadcast(pool->cond_new);
   cond_broadcast(pool->cond_fetch);
   cond_broadcast(pool->cond_tile);
 }
 
 void
 pool_reset(struct pool* pool)
 {
   ASSERT(pool);
   mutex_lock(pool->mutex);
   pool->next_part_id = 0;
   mutex_unlock(pool->mutex);
 }