Point

1. Context

1 Context

The following data stored in the context:

Variable	Type	Description
`num`	`int64_t`	Total number of points
`coord`	`qmckl_matrix`	`num` × 3 matrix3

We consider that the matrix is stored 'transposed' and 'normal' corresponds to the 3 × num matrix.

1.1 Data structure

typedef struct qmckl_point_struct {
  int64_t      num;
  qmckl_matrix coord;
} qmckl_point_struct;

qmckl_exit_code qmckl_init_point(qmckl_context context);

qmckl_exit_code qmckl_init_point(qmckl_context context) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return false;
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);

  qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
  mem_info.size = sizeof(qmckl_point_struct);
  ctx->point = (qmckl_point_struct*) qmckl_malloc(context, mem_info);
  if (ctx->point == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_init_point",
                           NULL);
  }
  memset(ctx->point, 0, sizeof(qmckl_point_struct)); 

  return QMCKL_SUCCESS;
}

1.2 Access functions

Access functions return QMCKL_SUCCESS when the data has been successfully retrieved. They return QMCKL_INVALID_CONTEXT when the context is not a valid context. If the function returns successfully, the variable pointed by the pointer given in argument contains the requested data. Otherwise, this variable is untouched.

1.2.1 Number of points

Returns the number of points stored in the context.

interface
  integer(c_int32_t) function qmckl_get_point_num(context, num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none

    integer (c_int64_t) , intent(in)  , value :: context                                  
    integer (c_int64_t) , intent(out)         :: num 
  end function
end interface

1.2.2 Point coordinates

Returns the point coordinates as sequences of (x,y,z). The pointer is assumed to point on a memory block of size size_max ≥ 3 * point_num.

interface
  integer(c_int32_t) function qmckl_get_point(context, transp, coord, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none

    integer (c_int64_t) , intent(in)  , value :: context
    character(c_char)   , intent(in)  , value :: transp
    real    (c_double ) , intent(out)         :: coord(*)
    integer (c_int64_t) , intent(in)          :: size_max
  end function
end interface

1.3 Initialization functions

When the data is set in the context, if the arrays are large enough, we overwrite the data contained in them.

To set the data relative to the points in the context, one of the following functions need to be called.

qmckl_exit_code qmckl_set_point (qmckl_context context,
                                 const char transp,
                                 const double* coord,
                                 const int64_t size_max);

Copy a sequence of (x,y,z) into the context.

qmckl_exit_code
qmckl_set_point (qmckl_context context,
                 const char transp,
                 const double* coord,
                 const int64_t num)
{

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return QMCKL_NULL_CONTEXT;
  }

  if (transp != 'N' && transp != 'T') {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_2,
                           "qmckl_set_point",
                           "transp should be 'N' or 'T'");
  }

  if (coord == NULL) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_ARG_3,
                           "qmckl_set_point",
                           "coord is a NULL pointer");
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct* const) context;
  assert (ctx != NULL);
  assert (ctx->point != NULL);

  qmckl_exit_code rc;
  if (ctx->point->num < num) {

    if (ctx->point->coord.data != NULL) {
      rc = qmckl_matrix_free(context, ctx->point->coord);
      return rc;
    }

    ctx->point->coord = qmckl_matrix_alloc(context, num, 3);
    if (ctx->point->coord.data == NULL) {
      return qmckl_failwith( context,
                             QMCKL_ALLOCATION_FAILED,
                             "qmckl_set_point",
                             NULL);
    }

  };

  ctx->point->num = num;

  if (transp == 'T') {
    memcpy(ctx->point->coord.data, coord, 3*num*sizeof(double));
  } else {
    for (int64_t i=0 ; i<num ; ++i) {
      qmckl_mat(ctx->point->coord, i, 0) = coord[3*i  ];
      qmckl_mat(ctx->point->coord, i, 1) = coord[3*i+1];
      qmckl_mat(ctx->point->coord, i, 2) = coord[3*i+2];
    }
  }

  return QMCKL_SUCCESS;

}

interface
  integer(c_int32_t) function qmckl_set_point(context, &
       transp, coord, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none

    integer (c_int64_t) , intent(in)  , value :: context
    character(c_char)   , intent(in)  , value :: transp
    real    (c_double ) , intent(in)          :: coord(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function
end interface

1.4 Test

/* Reference input data */
int64_t point_num = chbrclf_elec_num;
const double* coord     = &(chbrclf_elec_coord[0][0][0]);

/* --- */

qmckl_exit_code rc;
double coord2[point_num*3];
double coord3[point_num*3];


rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_NOT_PROVIDED);

rc = qmckl_set_point (context, 'N', coord, point_num);
assert(rc == QMCKL_SUCCESS);

int64_t n;
rc = qmckl_get_point_num (context, &n);
assert(rc == QMCKL_SUCCESS);
assert(n == point_num);

rc = qmckl_get_point (context, 'N', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i<3*point_num ; ++i) {
  assert( coord[i] == coord2[i] );
}

rc = qmckl_get_point (context, 'T', coord2, (point_num*3));
assert(rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i<point_num ; ++i) {
  assert( coord[3*i+0] == coord2[i] );
  assert( coord[3*i+1] == coord2[i+point_num] );
  assert( coord[3*i+2] == coord2[i+point_num*2] );
}

rc = qmckl_set_point (context, 'T', coord2, point_num);
assert(rc == QMCKL_SUCCESS);

rc = qmckl_get_point (context, 'N', coord3, (point_num*3));
assert(rc == QMCKL_SUCCESS);

for (int64_t i=0 ; i<3*point_num ; ++i) {
  assert( coord[i] == coord3[i] );
}