Lenstra-Lenstra-Lovász lattice basis reduction.

Former-commit-id: 3ada9f1ebf783d07c31671fd81cb6f7d6fe6187b

qpms/CMakeLists.txt

@@ -13,7 +13,9 @@ include_directories(${DIRS})
 add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c ewald.c
 	ewaldsf.c pointgroups.c latticegens.c
 	lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c
-	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c)
+	bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c
+	lll.c
+	)
 use_c99()
 
 set(LIBS ${LIBS} ${GSL_LIBRARIES} ${GSLCBLAS_LIBRARIES})

qpms/__init__.py

@@ -7,14 +7,14 @@ from sys import platform as __platform
 import warnings as __warnings
 
 try:
-    from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc
+    from .qpms_c import PointGroup, FinitePointGroup, FinitePointGroupElement, Particle, scatsystem_set_nthreads, ScatteringSystem, ScatteringMatrix, pitau, set_gsl_pythonic_error_handling, pgsl_ignore_error, gamma_inc, lll_reduce
 except ImportError as ex:
     if __platform == "linux" or __platform == "linux2":
         if 'LD_LIBRARY_PATH' not in __os.environ.keys():
             __warnings.warn("Environment variable LD_LIBRARY_PATH has not been set. Make it point to a directory where you installed libqpms and run python again")
         else:
             __warnings.warn("Does your LD_LIBRARY_PATH include a directory where you installed libqpms? Check and run python again."
-                'Currently, I see LD_LIBRARY_PATH="%s"' % __os.environ['LD_LIBRARY_PATH'])        
+                '\nCurrently, I see LD_LIBRARY_PATH="%s"' % __os.environ['LD_LIBRARY_PATH'])        
     raise ex
 from .qpms_p import *
 from .cyquaternions import CQuat, IRot3

qpms/lattices.h

@@ -63,6 +63,18 @@ static inline point2d point2d_fromxy(const double x, const double y) {
 	return p;
 }
 
+/// Lattice basis reduction.
+/** This is currenty a bit naïve implementation of
+ * Lenstra-Lenstra-Lovász algorithm.
+ *
+ * The reduction happens in-place, i.e. the basis vectors in \a b are
+ * replaced with the reduced basis.
+ */
+int qpms_reduce_lattice_basis(double *b, ///< Array of dimension [bsize][ndim].
+		const size_t bsize, ///< Number of the basis vectors (dimensionality of the lattice).
+		const size_t ndim ///< Dimension of the space into which the lattice is embedded.
+		);
+
 /// Generic lattice point generator type.
 /**
  * A bit of OOP-in-C brainfuck here.

qpms/lll.c

@@ -0,0 +1,99 @@
+#include <cblas.h>
+#include <string.h>
+#include <math.h>
+#include <qpms_error.h>
+
+static inline size_t mu_index(size_t k, size_t j) {
+  return k * (k - 1) / 2 + j;
+}
+
+/// Gram-Schmidt orthogonalisation.
+/** Does not return the actual orthogonal basis (as it is not needed
+ * for the LLL algorithm as such)  but rather only
+ * the mu(i,j) coefficients and squared norms of the orthogonal vectors
+ */
+static void gram_schmidt(
+    double *mu, ///< Array of \f[ \mu_{k,j} = \frac{\vect{v}_i \cdot\vect{v}_i^*}{|\vect v_j^*|^2}\f] of length mu_index(bsize, 0), 
+    double *vstar_sqnorm, ///< Array of \f$ \abs{\vect v_i^*}^2 \f$ of length bsize.
+    const double *v, ///< Vectors to orthogonalise, size [bsize][ndim], 
+    const size_t bsize, ///< Size of the basis ( = dimensionality of the lattice)
+    const size_t ndim ///< Dimensionality of the space.
+    )
+{
+  double v_star[bsize][ndim];
+  for (size_t i = 0; i < bsize; ++i) {
+    memcpy(v_star[i], v+i*ndim, ndim*sizeof(double));
+    double parallel_part[ndim /*???*/];
+    memset(parallel_part, 0, sizeof(parallel_part));
+    for (size_t j = 0; j < i; ++j) {
+      double mu_numerator = cblas_ddot(ndim, v + i*ndim, 1, v_star[j], 1);
+      mu[mu_index(i, j)] = mu_numerator / vstar_sqnorm[j];
+      cblas_daxpy(ndim, mu[mu_index(i, j)], v_star[j], 1, parallel_part, 1);
+    }
+    cblas_daxpy(ndim, -1, parallel_part, 1, v_star[i], 1);
+    vstar_sqnorm[i] = cblas_ddot(ndim, v_star[i], 1, v_star[i], 1);
+  }
+}
+
+static inline double fsq(double x) { return x * x; };
+
+// A naïve implementation of Lenstra-Lenstra-Lovász algorithm.
+int qpms_reduce_lattice_basis(double *b, const size_t bsize, const size_t ndim)
+{
+  QPMS_ENSURE(bsize <= ndim, 
+      "The basis must have less elements (%zd) than the space dimension (%zd).",
+      bsize, ndim);
+  double mu[mu_index(bsize,0)];
+  double bstar_sqnorm[bsize];
+  gram_schmidt(mu, bstar_sqnorm, b, bsize, ndim);
+  size_t k = 1;
+  while (k < bsize) {
+    // Step 1 of LLL, achieve mu(k, k-1) <= 0.5
+    if (fabs(mu[mu_index(k, k-1)]) > 0.5) {
+      double r = round(mu[mu_index(k, k-1)]);
+      // "normalize" b(k), replacing it with b(k) - r b(k-1)
+      cblas_daxpy(ndim, -r, b+(k-1)*ndim, 1, b+k*ndim, 1);
+      // update mu to correspond to the new b(k)
+      for(size_t j = 0; j < bsize; ++j)
+        mu[mu_index(k, j)] -= r*mu[mu_index(k-1, j)];
+    }
+    // Step 2
+    if (k > 0 && // Case 1
+        bstar_sqnorm[k] < (0.75 - fsq(mu[mu_index(k, k-1)])) * bstar_sqnorm[k-1]) {
+      // swap b(k) and b(k-1)
+      cblas_dswap(ndim, &b[k*ndim], 1, &b[(k-1)*ndim], 1);
+      double B_k = bstar_sqnorm[k];
+      double mu_kkm1_old = mu[mu_index(k, k-1)];
+      double C = B_k + fsq(mu_kkm1_old) * bstar_sqnorm[k-1];
+      mu[mu_index(k, k-1)] *= bstar_sqnorm[k-1] / C;
+      bstar_sqnorm[k] *= bstar_sqnorm[k-1] / C;
+      bstar_sqnorm[k-1] = C;
+      for(size_t j = k+1; j < bsize; ++j) {
+        double m = mu[mu_index(j, k-1)];
+        mu[mu_index(j, k-1)] = m*m + mu[mu_index(j, k)] * B_k / C;
+        mu[mu_index(j, k)] = m - mu[mu_index(j, k)] * mu_kkm1_old;
+      }
+      for(size_t j = 0; j < k-1; ++j) {
+        double m = mu[mu_index(k-1, j)];
+        mu[mu_index(k-1, j)] = mu[mu_index(k, j)];
+        mu[mu_index(k, j)] = m;
+      }
+      --k;
+    } else { // Case 2
+      size_t l = k;
+      while(l > 0) {
+        --l;
+        if(fabs(mu[mu_index(k, l)] > 0.5)) {
+          double r = round(mu[mu_index(k, l)]);
+          cblas_daxpy(ndim, -r, b+l*ndim, 1, b+k*ndim, 1);
+          for (size_t j = 0; j < l; ++j)
+            mu[mu_index(k, j)] -= r * mu[mu_index(l, j)];
+          mu[mu_index(k, l)] -= r;
+          l = k;
+        }
+      }
+      ++k;
+    }
+  }
+  return 0;
+}

qpms/qpms_c.pyx

@@ -609,3 +609,30 @@ def gamma_inc_CF(double a, cdouble x):
     with pgsl_ignore_error(15): #15 is underflow
         cx_gamma_inc_CF_e(a, x, &res)
     return (res.val, res.err)
+
+def lll_reduce(basis):
+    """
+    Lattice basis reduction with the Lenstra-Lenstra-Lovász algorithm.
+
+    basis is array_like with dimensions (n, d), where
+    n is the size of the basis (dimensionality of the lattice)
+    and d is the dimensionality of the space into which the lattice
+    is embedded.
+    """
+    basis = np.array(basis, copy=True, order='C')
+    if len(basis.shape) != 2:
+        raise ValueError("Expected two-dimensional array (got %d-dimensional)"
+                % len(basis.shape))
+    cdef size_t n, d
+    n, d = basis.shape
+    if n > d:
+        raise ValueError("Real space dimensionality (%d) cannot be smaller than"
+                "the dimensionality of the lattice (%d) embedded into it."
+                % (d, n))
+    cdef double [:,:] basis_view = basis
+    if 0 != qpms_reduce_lattice_basis(&basis_view[0,0], n, d):
+        raise RuntimeError("Something weird happened")
+    return basis
+
+
+

qpms/qpms_cdefs.pxd

@@ -209,6 +209,7 @@ cdef extern from "lattices.h":
         PGEN_1D_INC_TOWARDS_ORIGIN
     PGen PGen_1D_new_minMaxR(double period, double offset, double minR, bint inc_minR,
             double maxR, bint inc_maxR, PGen_1D_incrementDirection incdir)
+    int qpms_reduce_lattice_basis(double *b, size_t bsize, size_t ndim)
 
 
 cdef extern from "quaternions.h":