Option to use Eigen versus Lapack to solve equations

On Windows system, it is often tedious to install Lapack/Atlas/Openblas
packages so as to simulate PSFs. Configure the project:

```
meson setup -Duse_eigen=true -Darmadillo-code:lapack=none
build-windows/
```

In order to bypass LAPACK depdendencies. Download Eigen3 the header-only
library. Implement the Bidiagonal divide and conquer SVD solver. Write a
wrapper around the Eigen logic to avoid C++ namespace pollusion.

.github/workflows/validate-build.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
       matrix:
         #os: [ubuntu-20.04, ubuntu-latest, windows-latest, macos-latest]
-        os: [ubuntu-20.04, ubuntu-latest]
+        os: [ubuntu-20.04, ubuntu-latest, windows-latest]
         use_boost: ["false", "true"]
 
     steps:
@@ -54,14 +54,17 @@ jobs:
         build/
       if: startsWith(matrix.os, 'window') == false
 
-    #- name: Resolve C++ build dependencies (msvc toolchain)
-    #  run: meson setup --vsenv build/
-    #  if: startsWith(matrix.os, 'window') == true
+    - name: Resolve C++ build dependencies (msvc toolchain)
+      run: >-
+        meson setup
+        --vsenv
+        -Darmadillo-code:lapack=none
+        -Duse_eigen=true
+        build/
+      if: startsWith(matrix.os, 'window') == true
 
     - name: Build everything
-      #run: meson compile -C build/
-      run: ninja -C build/ all
+      run: meson compile -C build/
 
     - name: Test everything
-      #run: meson test -C build/
-      run: ninja -C build/ test
+      run: meson test -C build/

3rdparty/packagefiles/armadillo-code/meson.build

@@ -5,7 +5,10 @@ project('armadillo', 'cpp',
 # Header-only library
 armadillo_inc = include_directories('include')
 
-if get_option('lapack') == 'openblas64'
+if get_option('lapack') == 'none'
+  warning('Disabled Lapack dependencies. SVD and LU algorithms will be undefined.')
+  lapack_dep = dependency('', required: false)
+elif get_option('lapack') == 'openblas64'
   lapack_dep = dependency('openblas64')
 elif get_option('lapack') == 'openblas'
   lapack_dep = dependency('openblas')

3rdparty/packagefiles/armadillo-code/meson.options

@@ -6,6 +6,7 @@ option('lapack',
         'openblas',
         'lapack64',
         'lapack',
+        'none',
     ],
     value: 'openblas',
 )

besselj_armadillo_support/meson.build

@@ -1,5 +1,6 @@
 armadillo_besselj_support_inc = include_directories('inc')
-armadillo_dep = subproject('armadillo-code').get_variable('armadillo_dep')
+armadillo_proj = subproject('armadillo-code')
+armadillo_dep = armadillo_proj.get_variable('armadillo_dep')
 
 # TODO(Antony): use "cxx.has_function" when it supports std namespace.
 if get_option('use_boost')

examples/generate-psf.cpp

@@ -2,6 +2,14 @@
 
 #include "make_psf.h"
 
+namespace {
+#ifdef LSTSQ_USE_LAPACK
+constexpr bool has_lapack = true;
+#else
+constexpr bool has_lapack = false;
+#endif
+}  // namespace
+
 int main() {
     using arma::hdf5_name;
     using microsc_psf::makePSF;
@@ -20,6 +28,7 @@ int main() {
     precision_li2017_t precision{};
     precision.num_basis = 153;
     precision.rho_samples = 1000;
+    precision.solver = has_lapack ? microsc_psf::SandersonAndCurtin2020 : microsc_psf::EigenBdcSVD;
 
     const auto psf =
         makePSF(params, {0.1_um, 0.25_um}, {256, 128}, 0.610_um, precision);

examples/meson.build

@@ -12,7 +12,10 @@ endif
 
 generate_psf_exe = executable('generate-psf',
     sources: 'generate-psf.cpp',
-    cpp_args: generate_psf_compile_args,
+    cpp_args: [
+        generate_psf_compile_args,
+        use_lapack_args,
+    ],
     dependencies: [
         microsc_psf_dep,
         hdf5_dep,

meson.options

@@ -1,2 +1,11 @@
 option('install_examples', type: 'boolean', value: false)
-option('use_boost', type: 'boolean', value: false)
+option('use_boost',
+    description: 'Compute bessel function of the first kind with Boost::Math',
+    type: 'boolean',
+    value: false,
+)
+option('use_eigen',
+    description: 'Solve linear least squares with Eigen3 the matrix library',
+    type: 'boolean',
+    value: false,
+)

microsc-psf/inc/linsolver.h

@@ -0,0 +1,23 @@
+#pragma once
+
+#include <armadillo>
+
+namespace microsc_psf {
+namespace internal {
+
+/** Solve min || A * x - b ||^2 for x , using Eigen3 solver.
+ *
+ * Typical ussage: simulate PSF on machines not having LAPACK support, e.g.
+ * Windows on ARM64 CPU.
+ *
+ * Note: This may seem contrived to have two near-identical C++ libraries
+ * serving the linear algebra code, but Armadillo has a nice Matlab-like syntax
+ * that simplifies the code review process a lot.
+ *
+ * @tparam tranpose_b True if the signal b needs a Hermitian transpose ahead of
+ * the least squares solver.
+ */
+template <bool transpose_b>
+arma::cx_mat solveWithEigen(const arma::mat& A, const arma::cx_mat& b);
+}  // namespace internal
+}  // namespace microsc_psf

microsc-psf/inc/make_psf.h

@@ -30,8 +30,14 @@ enum linear_solver_option_t : uint8_t {
     // https://arma.sourceforge.net/armadillo_solver_2020.pdf
     SandersonAndCurtin2020,
 
-    // Penrose pseudo inverse.
-    PenroseInverse
+    /** Penrose pseudo inverse. */
+    PenroseInverse,
+
+    /** Eigen's SVD solver.
+     *
+     * Reference: https://eigen.tuxfamily.org/dox/classEigen_1_1BDCSVD.html
+     */
+    EigenBdcSVD,
 };
 
 struct precision_li2017_t {

microsc-psf/meson.build

@@ -1,10 +1,44 @@
 microsc_psf_inc = include_directories('inc')
+
+if get_option('use_eigen')
+    warning('Solving linear least squares with Eigen3, not Lapack')
+
+    # Eigen3 is a direct competitor of Armadillo library. This may seem
+    # contrived to require two nearly identicial libraries to compute PSF, but
+    # Armadillo has a Matlab-like syntax that simplifies the code review process
+    # a lot.
+    eigen_solver_lib = static_library('eigen_solver',
+        sources: 'src/linsolver.cpp',
+        include_directories: include_directories('inc'),
+        dependencies: [
+            dependency('eigen', fallback: ['eigen', 'eigen_dep']),
+            armadillo_dep,
+        ],
+    )
+
+    eigen_solver_dep = declare_dependency(
+        link_with: eigen_solver_lib,
+        include_directories: include_directories('inc'),
+        compile_args: '-DLSTSQ_USE_EIGEN',
+    )
+else
+    eigen_solver_dep = []
+endif
+
+lapack_dep = armadillo_proj.get_variable('lapack_dep')
+if lapack_dep.found()
+    use_lapack_args = '-DLSTSQ_USE_LAPACK'
+else
+    use_lapack_args = []
+endif
+
 microsc_psf_lib = static_library('microsc-psf',
     sources: [
         'src/main.cpp',
     ],
     cpp_args: [
         '-D' + bessel_source,
+        use_lapack_args,
     ],
     include_directories: [
         microsc_psf_inc,
@@ -13,6 +47,7 @@ microsc_psf_lib = static_library('microsc-psf',
     dependencies: [
         armadillo_dep,
         boost_math_dep,
+        eigen_solver_dep,
     ],
 )
 

microsc-psf/src/linsolver.cpp

@@ -0,0 +1,53 @@
+#include "linsolver.h"
+
+#include <Eigen/Dense>
+#include <cassert>
+
+namespace microsc_psf {
+
+namespace internal {
+
+template <bool transpose_b>
+arma::cx_mat
+solveWithEigen(const arma::mat& A_buffer, const arma::cx_mat& b_buffer) {
+    using arma::cx_double;
+    using Eigen::ComputeThinU;
+    using Eigen::ComputeThinV;
+    using Eigen::Map;
+    using Eigen::MatrixX;
+    using Eigen::MatrixXd;
+
+    const auto m = A_buffer.n_rows;
+    const auto n = A_buffer.n_cols;
+    const auto k = (transpose_b ? b_buffer.n_rows : b_buffer.n_cols);
+
+    if constexpr (transpose_b) {
+        assert(b_buffer.n_cols == m);
+    } else {
+        assert(b_buffer.n_rows == m);
+    }
+
+    // First, map to Eigen's primary data structure.
+    const Map<const MatrixXd> A(A_buffer.memptr(), m, n);
+    const Map<const MatrixX<cx_double>> b(b_buffer.memptr(), b_buffer.n_rows, b_buffer.n_cols);
+
+    // Allocate the result buffer
+    arma::cx_mat xopt_buffer(n, k);
+
+    // Next, solve with Eigen's SVD solver
+    Map<Eigen::MatrixX<cx_double>> xopt(xopt_buffer.memptr(), n, k);
+    if constexpr (transpose_b) {
+        // Eigen may have an efficient Hermitian transposed solver. Move the transpose step into
+        // Eigen.
+        xopt = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b.transpose());
+    } else {
+        xopt = A.bdcSvd(ComputeThinU | ComputeThinV).solve(b);
+    }
+
+    return xopt_buffer;
+}
+
+template arma::cx_mat solveWithEigen<true>(const arma::mat&, const arma::cx_mat&);
+template arma::cx_mat solveWithEigen<false>(const arma::mat&, const arma::cx_mat&);
+}  // namespace internal
+}  // namespace microsc_psf

microsc-psf/src/main.cpp

@@ -4,6 +4,10 @@
 // This must come after <boost/math/special_functions/bessel.hpp>
 #include "make_psf.h"
 
+#ifdef LSTSQ_USE_EIGEN
+#include "linsolver.h"
+#endif
+
 namespace {
 
 // Return the range [0, N), excluding N.
@@ -137,7 +141,7 @@ makePSF(microscope_params_t params, scale_t<Micron> voxel, scale_t<uint32_t> vol
     {
         // Define the basis of Bessel functions.
         // Shape: number of basis function by number of rho samples.
-        auto&& J = besselj<0>(scaling_factor * Rho);
+#define J besselj<0>(scaling_factor * Rho)
 
         // Compute the approximation to the sampled pupil phase by finding the least squares
         // solution to the complex coefficients of the Fourier-Bessel expansion.
@@ -146,10 +150,35 @@ makePSF(microscope_params_t params, scale_t<Micron> voxel, scale_t<uint32_t> vol
         //
         // Note: Armadillo does not have solver for real-valued matrix and complex-valued vector.
         // Reference: https://arma.sourceforge.net/armadillo_solver_2020.pdf
-        cx_mat&& Ci = (precision.solver == SandersonAndCurtin2020)
-                          ? solve(conv_to<cx_mat>::from(J.t()), phase.t())
-                          : cx_mat(pinv(J.t()) * phase.t());
-
+        cx_mat Ci;
+        switch (precision.solver) {
+#ifdef LSTSQ_USE_LAPACK
+            case SandersonAndCurtin2020:
+                Ci = solve(conv_to<cx_mat>::from(J.t()), phase.t());
+                break;
+            case PenroseInverse:
+                Ci = pinv(J.t()) * phase.t();
+                break;
+#else
+            case SandersonAndCurtin2020:
+            case PenroseInverse:
+                throw std::invalid_argument(R"(LAPACK not found.
+Re-configure the C++ project with
+"meson configure -Darmadillo-code:lapack=..." and try again.)");
+                break;
+#endif
+
+            case EigenBdcSVD: {
+#ifndef LSTSQ_USE_EIGEN
+                throw std::invalid_argument(R"(Eigen::BdcSVD solver not implemented.
+Re-configure the C++ project with
+"meson configure -Duse_eigen=true" and try again.)");
+#else
+                constexpr bool always_transpose_phase = true;
+                Ci = microsc_psf::internal::solveWithEigen<always_transpose_phase>(J.t(), phase);
+#endif
+            }
+        }
         const cx_mat ciEle = Ele * Ci;
         PSF0 = real(ciEle % conj(ciEle));
     }