From e3a08947bd72ec0c89a16169a3724b781e93849e Mon Sep 17 00:00:00 2001 From: vijay gopal chilkuri Date: Tue, 25 May 2021 17:48:25 +0530 Subject: [PATCH] Working on rescaled distances. #15 --- org/qmckl_distance.org | 267 +++++++++++++++++++++++++++++++++++++++++ org/qmckl_nucleus.org | 170 ++++++++++++++++++++++++++ 2 files changed, 437 insertions(+) diff --git a/org/qmckl_distance.org b/org/qmckl_distance.org index 4ea0e2b..8c5c191 100644 --- a/org/qmckl_distance.org +++ b/org/qmckl_distance.org @@ -775,6 +775,273 @@ end function test_qmckl_dist qmckl_exit_code test_qmckl_dist(qmckl_context context); assert(test_qmckl_dist(context) == QMCKL_SUCCESS); #+end_src + +* Rescaled Distance + +** ~qmckl_distance_rescaled~ + :PROPERTIES: + :Name: qmckl_distance_rescaled + :CRetType: qmckl_exit_code + :FRetType: qmckl_exit_code + :END: + + ~qmckl_distance_rescaled~ computes the matrix of the rescaled distances between all + pairs of points in two sets, one point within each set: + + \[ + C_{ij} = TODO + \] + + If the input array is normal (~'N'~), the xyz coordinates are in + the leading dimension: ~[n][3]~ in C and ~(3,n)~ in Fortran. + + #+NAME: qmckl_distance_rescaled_args + | qmckl_context | context | in | Global state | + | char | transa | in | Array ~A~ is ~'N'~: Normal, ~'T'~: Transposed | + | char | transb | in | Array ~B~ is ~'N'~: Normal, ~'T'~: Transposed | + | int64_t | m | in | Number of points in the first set | + | int64_t | n | in | Number of points in the second set | + | double | A[][lda] | in | Array containing the $m \times 3$ matrix $A$ | + | int64_t | lda | in | Leading dimension of array ~A~ | + | double | B[][ldb] | in | Array containing the $n \times 3$ matrix $B$ | + | int64_t | ldb | in | Leading dimension of array ~B~ | + | double | C[n][ldc] | out | Array containing the $m \times n$ matrix $C$ | + | int64_t | ldc | in | Leading dimension of array ~C~ | + +*** Requirements + + - ~context~ is not ~QMCKL_NULL_CONTEXT~ + - ~m > 0~ + - ~n > 0~ + - ~lda >= 3~ if ~transa == 'N'~ + - ~lda >= m~ if ~transa == 'T'~ + - ~ldb >= 3~ if ~transb == 'N'~ + - ~ldb >= n~ if ~transb == 'T'~ + - ~ldc >= m~ + - ~A~ is allocated with at least $3 \times m \times 8$ bytes + - ~B~ is allocated with at least $3 \times n \times 8$ bytes + - ~C~ is allocated with at least $m \times n \times 8$ bytes + +*** C header + + #+CALL: generate_c_header(table=qmckl_distance_rescaled_args,rettyp=get_value("CRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src c :tangle (eval h_func) :comments org + qmckl_exit_code qmckl_distance_rescaled ( + const qmckl_context context, + const char transa, + const char transb, + const int64_t m, + const int64_t n, + const double* A, + const int64_t lda, + const double* B, + const int64_t ldb, + double* const C, + const int64_t ldc ); + #+end_src + +*** Source + #+begin_src f90 :tangle (eval f) +integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, & + A, LDA, B, LDB, C, LDC) & + result(info) + use qmckl + implicit none + integer(qmckl_context) , intent(in) :: context + character , intent(in) :: transa, transb + integer*8 , intent(in) :: m, n + integer*8 , intent(in) :: lda + real*8 , intent(in) :: A(lda,*) + integer*8 , intent(in) :: ldb + real*8 , intent(in) :: B(ldb,*) + integer*8 , intent(in) :: ldc + real*8 , intent(out) :: C(ldc,*) + + integer*8 :: i,j + real*8 :: x, y, z + integer :: transab + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (m <= 0_8) then + info = QMCKL_INVALID_ARG_4 + return + endif + + if (n <= 0_8) then + info = QMCKL_INVALID_ARG_5 + return + endif + + if (transa == 'N' .or. transa == 'n') then + transab = 0 + else if (transa == 'T' .or. transa == 't') then + transab = 1 + else + transab = -100 + endif + + if (transb == 'N' .or. transb == 'n') then + continue + else if (transa == 'T' .or. transa == 't') then + transab = transab + 2 + else + transab = -100 + endif + + if (transab < 0) then + info = QMCKL_INVALID_ARG_1 + return + endif + + if (iand(transab,1) == 0 .and. LDA < 3) then + info = QMCKL_INVALID_ARG_7 + return + endif + + if (iand(transab,1) == 1 .and. LDA < m) then + info = QMCKL_INVALID_ARG_7 + return + endif + + if (iand(transab,2) == 0 .and. LDA < 3) then + info = QMCKL_INVALID_ARG_7 + return + endif + + if (iand(transab,2) == 2 .and. LDA < m) then + info = QMCKL_INVALID_ARG_7 + return + endif + + + select case (transab) + + case(0) + + do j=1,n + do i=1,m + x = A(1,i) - B(1,j) + y = A(2,i) - B(2,j) + z = A(3,i) - B(3,j) + C(i,j) = x*x + y*y + z*z + end do + C(:,j) = dsqrt(C(:,j)) + end do + + case(1) + + do j=1,n + do i=1,m + x = A(i,1) - B(1,j) + y = A(i,2) - B(2,j) + z = A(i,3) - B(3,j) + C(i,j) = x*x + y*y + z*z + end do + C(:,j) = dsqrt(C(:,j)) + end do + + case(2) + + do j=1,n + do i=1,m + x = A(1,i) - B(j,1) + y = A(2,i) - B(j,2) + z = A(3,i) - B(j,3) + C(i,j) = x*x + y*y + z*z + end do + C(:,j) = dsqrt(C(:,j)) + end do + + case(3) + + do j=1,n + do i=1,m + x = A(i,1) - B(j,1) + y = A(i,2) - B(j,2) + z = A(i,3) - B(j,3) + C(i,j) = x*x + y*y + z*z + end do + C(:,j) = dsqrt(C(:,j)) + end do + + end select + +end function qmckl_distance_rescaled_f + #+end_src + +*** Performance + + This function is more efficient when ~A~ and ~B~ are transposed. + +** C interface :noexport: + + #+CALL: generate_c_interface(table=qmckl_distance_rescaled_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_distance_rescaled & + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + character , intent(in) , value :: transa + character , intent(in) , value :: transb + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(out) :: C(ldc,n) + integer (c_int64_t) , intent(in) , value :: ldc + + integer(c_int32_t), external :: qmckl_distance_rescaled_f + info = qmckl_distance_rescaled_f & + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) + + end function qmckl_distance_rescaled + #+end_src + + #+CALL: generate_f_interface(table=qmckl_distance_rescaled_args,rettyp=get_value("FRetType"),fname=get_value("Name")) + + #+RESULTS: + #+begin_src f90 :tangle (eval fh_func) :comments org :exports none + interface + integer(c_int32_t) function qmckl_distance_rescaled & + (context, transa, transb, m, n, A, lda, B, ldb, C, ldc) & + bind(C) + use, intrinsic :: iso_c_binding + import + implicit none + + integer (c_int64_t) , intent(in) , value :: context + character , intent(in) , value :: transa + character , intent(in) , value :: transb + integer (c_int64_t) , intent(in) , value :: m + integer (c_int64_t) , intent(in) , value :: n + real (c_double ) , intent(in) :: A(lda,*) + integer (c_int64_t) , intent(in) , value :: lda + real (c_double ) , intent(in) :: B(ldb,*) + integer (c_int64_t) , intent(in) , value :: ldb + real (c_double ) , intent(out) :: C(ldc,n) + integer (c_int64_t) , intent(in) , value :: ldc + + end function qmckl_distance_rescaled + end interface + #+end_src + +*** Test :noexport: * End of files :noexport: #+begin_src c :comments link :tangle (eval c_test) diff --git a/org/qmckl_nucleus.org b/org/qmckl_nucleus.org index 834535b..b09cf6f 100644 --- a/org/qmckl_nucleus.org +++ b/org/qmckl_nucleus.org @@ -79,9 +79,11 @@ typedef struct qmckl_nucleus_struct { int64_t num; int64_t repulsion_date; int64_t nn_distance_date; + int64_t nn_distance_rescaled_date; double* coord; double* charge; double* nn_distance; + double* nn_distance_rescaled; double repulsion; int32_t uninitialized; bool provided; @@ -639,6 +641,174 @@ assert(fabs(distance[1]-2.070304721365169) < 1.e-12); #+end_src + +** Nucleus-nucleus rescaled distances + +*** Get + + #+begin_src c :comments org :tangle (eval h_func) :noweb yes +qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_get_nucleus_nn_distance_rescaled(qmckl_context context, double* distance_rescaled) +{ + /* Check input parameters */ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + qmckl_exit_code rc = qmckl_provide_nn_distance_rescaled(context); + if (rc != QMCKL_SUCCESS) return rc; + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + size_t sze = ctx->nucleus.num * ctx->nucleus.num; + memcpy(distance_rescaled, ctx->nucleus.nn_distance_rescaled, sze * sizeof(double)); + + return QMCKL_SUCCESS; +} + #+end_src + +*** Provide :noexport: + + #+begin_src c :comments org :tangle (eval h_private_func) :noweb yes :exports none +qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context); + #+end_src + + #+begin_src c :comments org :tangle (eval c) :noweb yes :exports none +qmckl_exit_code qmckl_provide_nn_distance_rescaled(qmckl_context context) +{ + /* Check input parameters */ + if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) { + return (char) 0; + } + + qmckl_context_struct* const ctx = (qmckl_context_struct* const) context; + assert (ctx != NULL); + + if (!ctx->nucleus.provided) return QMCKL_NOT_PROVIDED; + + /* Allocate array */ + if (ctx->nucleus.nn_distance_rescaled == NULL) { + + qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero; + mem_info.size = ctx->nucleus.num * ctx->nucleus.num * sizeof(double); + double* nn_distance_rescaled = (double*) qmckl_malloc(context, mem_info); + + if (nn_distance_rescaled == NULL) { + return qmckl_failwith( context, + QMCKL_ALLOCATION_FAILED, + "qmckl_nn_distance_rescaled", + NULL); + } + ctx->nucleus.nn_distance_rescaled = nn_distance_rescaled; + } + + qmckl_exit_code rc = + qmckl_compute_nn_distance_rescaled(context, + ctx->nucleus.num, + ctx->nucleus.coord, + ctx->nucleus.nn_distance_rescaled); + if (rc != QMCKL_SUCCESS) { + return rc; + } + + ctx->nucleus.nn_distance_rescaled_date = ctx->date; + + return QMCKL_SUCCESS; +} + #+end_src + +*** Compute + + #+NAME: qmckl_nn_distance_rescaled_args + | qmckl_context | context | in | Global state | + | int64_t | nucl_num | in | Number of nuclei | + | double | coord[3][nucl_num] | in | Nuclear coordinates (au) | + | double | nn_distance_rescaled[nucl_num][nucl_num] | out | Nucleus-nucleus rescaled distances (au) | + + #+begin_src f90 :comments org :tangle (eval f) :noweb yes +integer function qmckl_compute_nn_distance_rescaled_f(context, nucl_num, coord, nn_distance_rescaled) & + result(info) + use qmckl + implicit none + integer(qmckl_context), intent(in) :: context + integer*8 , intent(in) :: nucl_num + double precision , intent(in) :: coord(nucl_num,3) + double precision , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num) + + integer*8 :: k + + info = QMCKL_SUCCESS + + if (context == QMCKL_NULL_CONTEXT) then + info = QMCKL_INVALID_CONTEXT + return + endif + + if (nucl_num <= 0) then + info = QMCKL_INVALID_ARG_2 + return + endif + + info = qmckl_distance_rescaled(context, 'T', 'T', nucl_num, nucl_num, & + coord, nucl_num, & + coord, nucl_num, & + nn_distance_rescaled, nucl_num) + +end function qmckl_compute_nn_distance_rescaled_f + #+end_src + + #+begin_src c :tangle (eval h_private_func) :comments org :exports none +qmckl_exit_code qmckl_compute_nn_distance_rescaled ( + const qmckl_context context, + const int64_t nucl_num, + const double* coord, + double* const nn_distance_rescaled ); + #+end_src + + + #+CALL: generate_c_interface(table=qmckl_nn_distance_rescaled_args,rettyp="qmckl_exit_code",fname="qmckl_compute_nn_distance") + + #+RESULTS: + #+begin_src f90 :tangle (eval f) :comments org :exports none + integer(c_int32_t) function qmckl_compute_nn_distance_rescaled & + (context, nucl_num, coord, nn_distance_rescaled) & + bind(C) result(info) + + use, intrinsic :: iso_c_binding + implicit none + + integer (c_int64_t) , intent(in) , value :: context + integer (c_int64_t) , intent(in) , value :: nucl_num + real (c_double ) , intent(in) :: coord(nucl_num,3) + real (c_double ) , intent(out) :: nn_distance_rescaled(nucl_num,nucl_num) + + integer(c_int32_t), external :: qmckl_compute_nn_distance_rescaled_f + info = qmckl_compute_nn_distance_rescaled_f & + (context, nucl_num, coord, nn_distance_rescaled) + + end function qmckl_compute_nn_distance_rescaled + #+end_src + +*** Test + + #+begin_src c :tangle (eval c_test) +/* Reference input data */ +/* TODO */ + +//assert(qmckl_nucleus_provided(context)); +// +//double distance[nucl_num*nucl_num]; +//rc = qmckl_get_nucleus_nn_distance(context, distance); +//assert(distance[0] == 0.); +//assert(distance[1] == distance[nucl_num]); +//assert(fabs(distance[1]-2.070304721365169) < 1.e-12); + + #+end_src + ** Nuclear repulsion energy \[