Consolidate vector mathematics into a single C++ file native/c_vectormath.cpp. Make some of these external C++ declarations available in other modules by using a Cython pxd file, vectormath.pxd.

Co-Authored-By: David Ignacio Cortee <[email protected]>

Author: rpep

Makefile

@@ -20,7 +20,8 @@ CPPFLAGS = -fPIC -g -Iinclude/ \
            -lfftw3_omp -lfftw3 \
            -lsundials_cvodes -lsundials_nvecserial -lsundials_nvecopenmp \
            -lblas -llapack \
-           -fopenmp
+           -fopenmp \
+           -std=c++14
 
 LDFLAGS = -shared
 SOURCES = $(shell echo native/src/*.cpp)

fidimag/atomistic/a_clib.pyx

@@ -1,9 +1,11 @@
 # distutils: language = c++
-
 import numpy
 cimport numpy as np
+cimport fidimag.common.vectormath as vmath
 np.import_array()
 
+
+
 cdef extern from "a_random.h":
     ctypedef struct mt19937_state:
         pass
@@ -97,8 +99,6 @@ cdef extern from "a_clib.h":
 
     # -------------------------------------------------------------------------
 
-    void normalise(double *m, int *pins, int n)
-
     # used for sllg
     void llg_rhs_dw_c(double *m, double *h, double *dm, double *T, double *alpha,
                       double *mu_s_inv, int *pins, double *eta, int n,
@@ -272,7 +272,7 @@ def compute_demag_full(double [:] spin,
 
 def normalise_spin(np.ndarray[double, ndim=1, mode="c"] spin,
               np.ndarray[int, ndim=1, mode="c"] pins, n):
-    normalise(&spin[0], &pins[0], n)
+    vmath.normalise(&spin[0], &pins[0], n)
 
 def compute_llg_rhs_dw(double [:] dm,
                        double [:] spin,

fidimag/common/neb_method/nebm_clib.pyx

@@ -1,5 +1,7 @@
 # distutils: language = c++
 
+cimport fidimag.common.vectormath as vmath
+
 cdef extern from "c_nebm_lib.h":
     void compute_tangents_C(double *tangents,
                             double *y,
@@ -23,7 +25,7 @@ cdef extern from "c_nebm_lib.h":
                                    int n_images,
                                    int n_dofs_image)
 
-    void normalise(double * a, int n)
+
     void normalise_images_C(double * y, int n_images, int n_dofs_image)
     void normalise_spins_C(double * y, int n_images, int n_dofs_image)
     double compute_distance_cartesian(double * A, double * B, int n_dofs_image,
@@ -112,7 +114,7 @@ def compute_effective_force(double [:] G,
                               )
 
 def normalise_clib(double [:] a, n):
-    normalise(&a[0], n)
+    vmath.normalise(&a[0], n)
 
 def project_images(double [:]  vector,
                    double [:]  y,

fidimag/common/sundials/cvode.pyx

@@ -7,6 +7,7 @@ cimport openmp
 np.import_array()  # don't remove or you'll segfault
 from libc.string cimport memcpy
 import sys
+from fidimag.common cimport vectormath
 
 
 cdef extern from "sundials/sundials_types.h":
@@ -16,7 +17,7 @@ cdef extern from "sundials/sundials_types.h":
 cdef extern from "sundials/sundials_nvector.h":
     cdef struct _generic_N_Vector:
         void *content
-        
+
     ctypedef _generic_N_Vector *N_Vector
     N_Vector N_VNew_Serial(long int vec_length)
     N_Vector N_VNew_OpenMP(long int vec_length, int num_threads)
@@ -170,16 +171,10 @@ cdef extern from "sundials/sundials_iterative.h":
     int PREC_LEFT
     int PREC_RIGHT
     int PREC_BOTH
-    
+ 
     int MODIFIED_GS
     int CLASSICAL_GS
 
-
-
-cdef extern from "a_clib.h":
-    void normalise(double * m, int nxyz)
-
-
 cdef struct cv_userdata:
     void * rhs_fun
     void * y

fidimag/common/vectormath.pxd

@@ -0,0 +1,3 @@
+cdef extern from "c_vectormath.h":
+     void normalise(double *m, int *pins, int n)
+     void normalise(double * a, int n)

native/include/a_clib.h

@@ -4,27 +4,11 @@
 #include <fftw3.h>
 #include <math.h>
 #include<omp.h>
-
 #include "a_random.h"
+#include "c_vectormath.h"
 
 #define WIDE_PI 3.1415926535897932384626433832795L
 
-// ----------------------------------------------------------------------------
-/* 3 components for the cross product calculations */
-
-inline double cross_x(double a0, double a1, double a2, double b0, double b1,
-                      double b2) {
-  return a1 * b2 - a2 * b1;
-}
-inline double cross_y(double a0, double a1, double a2, double b0, double b1,
-                      double b2) {
-  return a2 * b0 - a0 * b2;
-}
-inline double cross_z(double a0, double a1, double a2, double b0, double b1,
-                      double b2) {
-  return a0 * b1 - a1 * b0;
-}
-
 // ----------------------------------------------------------------------------
 // From exch.c
 
@@ -101,8 +85,6 @@ void compute_px_py_c(double *spin, int nx, int ny, int nz, double *px,
 // ----------------------------------------------------------------------------
 // From sllg.c
 
-void normalise(double *m, int *pins, int n);
-
 void llg_rhs_dw_c(double * m, double * h, double * dm, double * T, double * alpha,
                   double * mu_s_inv, int *pins, double * eta, int n, double gamma,
                   double dt);

native/include/c_clib.h

@@ -4,42 +4,6 @@
 #include<omp.h>
 #define WIDE_PI 3.1415926535897932384626433832795L
 
-// ----------------------------------------------------------------------------
-
-/* 3 components for the cross product calculations */
-inline double cross_x(double a0, double a1, double a2, double b0, double b1,
-                      double b2) {
-  return a1 * b2 - a2 * b1;
-}
-inline double cross_y(double a0, double a1, double a2, double b0, double b1,
-                      double b2) {
-  return a2 * b0 - a0 * b2;
-}
-inline double cross_z(double a0, double a1, double a2, double b0, double b1,
-                      double b2) {
-  return a0 * b1 - a1 * b0;
-}
-
-inline void normalise(double *m, int *pins, int n){
-	int i, j, k;
-	double mm;
-	for (int id = 0; id < n; id++) {
-        i = 3*id;
-        j = i + 1;
-        k = j + 1;
-
-        if (pins[id]>0) continue;
-
-        mm = sqrt(m[i] * m[i] + m[j] * m[j] + m[k] * m[k]);
-        if(mm > 0) {
-            mm = 1 / mm;
-            m[i] *= mm;
-            m[j] *= mm;
-            m[k] *= mm;
-        }
-    }   
-}
-
 // ----------------------------------------------------------------------------
 // From: llg.c
 

native/include/c_dipolar.h

@@ -1,14 +1,11 @@
 #include<math.h>
 #include<complex>
 #include<fftw3.h>
+#include "c_vectormath.h"
 //#include<omp.h>
 
 #define WIDE_PI 3.1415926535897932384626433832795L
 
-inline double cross_x(double a0, double a1, double a2, double b0, double b1, double b2) { return a1*b2 - a2*b1; }
-inline double cross_y(double a0, double a1, double a2, double b0, double b1, double b2) { return a2*b0 - a0*b2; }
-inline double cross_z(double a0, double a1, double a2, double b0, double b1, double b2) { return a0*b1 - a1*b0; }
-
 
 enum Type_Nij {
 	Tensor_xx, Tensor_yy, Tensor_zz, Tensor_xy, Tensor_xz, Tensor_yz

native/include/c_nebm_lib.h

@@ -1,15 +1,6 @@
 #include "math.h"
 #define WIDE_PI 3.1415926535897932384626433832795L
 
-
-void cross_product(double * output, double * A, double * B);
-
-double dot_product(double * A, double * B, int n);
-
-double compute_norm(double * a, int n);
-
-void normalise(double * a, int n);
-
 void compute_tangents_C(double * ys, double * energy,
                         double * tangents, int image_num, int nodes
                         );

native/include/m_clib.h

@@ -5,10 +5,6 @@
 #define MU0 1.25663706143591728850e-6
 #define MU0_INV 795774.71545947669074
 
-double cross_x(double a0, double a1, double a2, double b0, double b1, double b2);
-double cross_y(double a0, double a1, double a2, double b0, double b1, double b2);
-double cross_z(double a0, double a1, double a2, double b0, double b1, double b2);
-
 void compute_exch_field_micro(double * m, double * field, double * energy, double * Ms_inv,
                          double A, double dx, double dy, double dz, int n, int *ngbs);