Point

1. Context

1 Context

The following data stored in the context:

Variable	Type	Description
`num`	`int64_t`	Total number of points
`date`	`uint64_t`	Last modification date of the coordinates
`coord`	`qmckl_matrix`	`num` × 3 matrix

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;
  uint64_t     date;
  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*) context;
  assert (ctx != 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. Here, num is the number of points to set.

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

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

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

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

  if (num <= 0) {
      return qmckl_failwith( context,
                             QMCKL_INVALID_ARG_3,
                             "qmckl_set_point",
                             "Number of points should be >0.");
  }

  if (size_max < 3*num) {
      return qmckl_failwith( context,
                             QMCKL_INVALID_ARG_4,
                             "qmckl_set_point",
                             "Array too small");
  }

  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*) context;
  assert (ctx != NULL);

  qmckl_exit_code rc;
  if (num != ctx->point.num) {

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

    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') {
    double *a = ctx->point.coord.data;
#ifdef HAVE_OPENMP
    #pragma omp for
#endif
    for (int64_t i=0 ; i<3*num ; ++i) {
      a[i] = coord[i];
    }
  } else {
#ifdef HAVE_OPENMP
    #pragma omp for
#endif
    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];
    }
  }

  /* Increment the date of the context */
  rc = qmckl_context_touch(context);
  assert (rc == QMCKL_SUCCESS);

  return QMCKL_SUCCESS;

}

interface
  integer(c_int32_t) function qmckl_set_point(context, &
       transp, num, 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
    integer (c_int64_t) , intent(in)  , value :: num
    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', point_num, coord, (point_num*3));
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', point_num, coord2, (point_num*3));
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] );
}