Scamp on the Mac

configure.ac

@@ -50,6 +50,7 @@ AC_SUBST(DATE3, "$date3")
 
 # Include macros
 sinclude(acx_atlas.m4)
+sinclude(acx_accelerate.m4)
 sinclude(acx_openblas.m4)
 sinclude(acx_curl.m4)
 sinclude(acx_fftw.m4)
@@ -101,6 +102,15 @@ AC_ARG_ENABLE(mkl,
 	AC_MSG_RESULT([yes]),
 	AC_MSG_RESULT([no]))
 
+# Provide special options for Apple Accelerate
+AC_MSG_CHECKING([whether we should use Appel's Accelerate])
+AC_ARG_ENABLE(accelerate,
+  [AS_HELP_STRING([--enable-accelerate],
+  [Use Apple's Accelerate for linear algebra (off by default)])],
+        CC="cc"
+  AC_MSG_RESULT([yes]),
+  AC_MSG_RESULT([no]))
+
 # Checks for programs.
 AC_LANG(C)
 
@@ -285,32 +295,36 @@ else
     ],
     AC_MSG_ERROR([$FFTW_ERROR Exiting.])
   )
-
   if test "x$enable_openblas" = "xyes"; then
 ######## Handle the OpenBLAS library (linear algebra: BLAS + LAPACKe) ########
     ACX_OPENBLAS($with_openblas_libdir, $with_openblas_incdir, $use_pthreads, no,
-      [
-        AM_CFLAGS="$AM_CFLAGS $OPENBLAS_CFLAGS "
-        AM_LDFLAGS="$AM_LDFLAGS $OPENBLAS_LDFLAGS "
-        LIBS="$OPENBLAS_LIBS $LIBS"
-        if test "$OPENBLAS_WARN" != ""; then
-          AC_MSG_WARN([$OPENBLAS_WARN])
-        fi
-      ],
-      AC_MSG_ERROR([$OPENBLAS_ERROR Exiting.])
+    [
+    AM_CFLAGS="$AM_CFLAGS $OPENBLAS_CFLAGS "
+    AM_LDFLAGS="$AM_LDFLAGS $OPENBLAS_LDFLAGS "
+    LIBS="$OPENBLAS_LIBS $LIBS"
+    if test "$OPENBLAS_WARN" != ""; then
+      AC_MSG_WARN([$OPENBLAS_WARN])
+    fi
+    ],
+    AC_MSG_ERROR([$OPENBLAS_ERROR Exiting.])
     )
-
   else
+################ handle the Apple Accelerate framework (linear algebra) ######
+    if test x$enable_accelerate = xyes; then
+      ACX_ACCEL()
+      LIBS="$ACCEL_LIBS $LIBS"
+    else
 ######### handle the ATLAS library (linear algebra: BLAS + cLAPACK) ##########
-    ACX_ATLAS($with_atlas_libdir, $with_atlas_incdir, $use_pthreads,
-      [
+      ACX_ATLAS($with_atlas_libdir, $with_atlas_incdir, $use_pthreads,
+        [
         [LIBS="$ATLAS_LIBS $LIBS"]
         if test "$ATLAS_WARN" != ""; then
-	  AC_MSG_WARN([$ATLAS_WARN])
-	fi
-      ],
-      AC_MSG_ERROR([$ATLAS_ERROR Exiting.])
-    )
+          AC_MSG_WARN([$ATLAS_WARN])
+        fi
+        ],
+        AC_MSG_ERROR([$ATLAS_ERROR Exiting.])
+      )
+    fi
   fi
 fi
 

m4/acx_accelerate.m4

@@ -0,0 +1,83 @@
+dnl
+dnl				acx_accelerate.m4
+dnl
+dnl Set up options for using the Mac OSX Accelerate Framework.
+dnl
+dnl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+dnl
+dnl	This file part of:	AstrOmatic software
+dnl
+dnl	Copyright:		(C) 2002-2020 IAP/CNRS/SorbonneU
+dnl
+dnl	License:		GNU General Public License
+dnl
+dnl	AstrOmatic software is free software: you can redistribute it and/or
+dnl	modify it under the terms of the GNU General Public License as
+dnl	published by the Free Software Foundation, either version 3 of the
+dnl	License, or (at your option) any later version.
+dnl	AstrOmatic software is distributed in the hope that it will be useful,
+dnl	but WITHOUT ANY WARRANTY; without even the implied warranty of
+dnl	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+dnl	GNU General Public License for more details.
+dnl	You should have received a copy of the GNU General Public License
+dnl	along with AstrOmatic software.
+dnl	If not, see <http://www.gnu.org/licenses/>.
+dnl
+dnl	Last modified:		12/08/2020
+dnl
+dnl %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+dnl
+dnl @synopsis ACX_ACCEL([MKL_DIR, ILP64_FLAG, STATIC_FLAG, CONV_LIBS])
+dnl
+dnl This macro sets the MKL_CFLAGS, MKL_LDFLAGS and MKL_LIBS variables to
+dnl for compiling and linking with INTEL's MKL. A coma-separated list of
+dnl convenience libraries may be included in the linked group for static linking.
+dnl You may wish to use these variables in your default CFLAGS:
+dnl
+dnl        CFLAGS="$CFLAGS $MKL_CFLAGS"
+dnl
+dnl You may wish to use these variables in your default LDFLAGS:
+dnl
+dnl        LDFLAGS="$LDFLAGS $MKL_LDLAGS"
+dnl
+dnl You may wish to use these variables in your default LIBS:
+dnl
+dnl        LIBS="$LIBS $MKL_LIBS"
+dnl
+
+AC_DEFUN([ACX_ACCEL], [
+AC_REQUIRE([AC_CANONICAL_HOST])
+
+
+dnl --------------------------
+dnl Exit if host is not MacOSX
+dnl --------------------------
+case $host_os in
+  darwin* ) ;;
+  *)
+    ACCEL_WARN="Accelerate only available on Mac OSX"
+    AC_SUBST(ACCEL_WARN)
+    exit
+esac
+
+acx_accelerate_ok=no
+AC_CHECK_HEADERS([Accelerate/Accelerate.h], [acx_accelerate_ok=yes])
+AC_SUBST(ACCEL_LIBS, "-framework Accelerate")
+AC_SUBST(MKL_LDFLAGS, "")
+
+dnl --------------------
+dnl Set internal flags
+dnl --------------------
+
+AC_DEFINE(HAVE_ACCELERATE,1, [Define if you have the Accelerate libraries.])
+AC_DEFINE(HAVE_LAPACK,1, [Define if you have the LAPACK libraries.])
+
+dnl --------------------
+dnl Set include files
+dnl --------------------
+
+AC_DEFINE(ACCELERATE_H, "Accelerate/Accelerate.h", [Accelerate header filename.])
+AC_DEFINE(LAPACK_H, "Accelerate/Accelerate.h", [LAPACK header filename.])
+
+])dnl ACX_ACCEL
+

src/astrsolve.c

@@ -46,6 +46,10 @@
 #include ATLAS_LAPACK_H
 #endif
 
+#ifdef HAVE_ACCELERATE
+#include ACCELERATE_H
+#endif
+
 #ifdef HAVE_LAPACKE
 #include LAPACKE_H
 #define MATSTORAGE_PACKED 1
@@ -528,6 +532,12 @@ void	astrsolve_fgroups(fgroupstruct **fgroups, int nfgroup)
 #endif
         free(lap_ipiv);
 
+#elif defined(HAVE_ACCELERATE)
+    int info,n_beta=1;
+    char* upper_or_lower="U";
+    if (dposv_(upper_or_lower, &ncoefftot,
+        &n_beta, alpha, &ncoefftot, beta, &ncoefftot, &info) != 0)
+      warning("Not a positive definite matrix", " in astrometry solver");
 
 #else
 

src/match.c

@@ -46,6 +46,10 @@
 #include ATLAS_LAPACK_H
 #endif
 
+#ifdef HAVE_ACCELERATE
+#include ACCELERATE_H
+#endif
+
 #ifdef HAVE_LAPACKE
 #include LAPACKE_H
 #endif
@@ -1315,6 +1319,17 @@ void match_refine(setstruct *set, setstruct *refset, double matchresol,
     {
         LAPACKE_dpotrs(LAPACK_COL_MAJOR, 'U', 3, 1, alpha, 3, blng, 3);
         LAPACKE_dpotrs(LAPACK_COL_MAJOR, 'U', 3, 1, alpha, 3, blat, 3);
+#elif defined(HAVE_ACCELERATE)
+  int info;
+  int N=3;
+  int LDA=3;
+  int LDB=3;
+  int NRHS=1;
+  char* upper_or_lower="U";
+  if (dpotrf_(upper_or_lower, &N, alpha, &LDA, &info) == 0)
+    {
+    dpotrs_(upper_or_lower, &N, &NRHS, alpha, &LDA, blng, &LDB, &info);
+    dpotrs_(upper_or_lower, &N, &NRHS, alpha, &LDA, blat, &LDB, &info);
 #else
         if (clapack_dpotrf(CblasRowMajor, CblasUpper, 3, alpha, 3) == 0)
         {

src/mosaic.c

@@ -58,10 +58,15 @@
 #include ATLAS_LAPACK_H
 #endif
 
+#ifdef HAVE_ACCELERATE
+#include ACCELERATE_H
+#endif
+
 #ifdef HAVE_LAPACKE
 #include LAPACKE_H
 #endif
 
+
 /*------------------- global variables for multithreading -------------------*/
 #ifdef USE_THREADS
  pthread_t		*mosaic_thread;
@@ -248,6 +253,10 @@ void	adjust_set(fieldstruct **fields, int nfield, int s)
 /* Derive the CRPIXs */
 #if defined(HAVE_LAPACKE)
   LAPACKE_dgesv(LAPACK_ROW_MAJOR, naxis, 1, cd, naxis, ipiv, x, 1);
+#elif defined(HAVE_ACCELERATE)
+  int NRHS=1;
+  int info;
+  dgesv_(&naxis, &NRHS, cd, &naxis, ipiv, x, &naxis, &info);
 #else
   clapack_dgesv(CblasRowMajor, naxis, 1, cd, naxis, ipiv, x, naxis);
 #endif

src/photsolve.c

@@ -50,6 +50,10 @@
 #include ATLAS_LAPACK_H
 #endif
 
+#ifdef HAVE_ACCELERATE
+#include ACCELERATE_H
+#endif
+
 #ifdef HAVE_LAPACKE
 #include LAPACKE_H
 #endif
@@ -378,7 +382,11 @@ void	photsolve_fgroups(fgroupstruct **fgroups, int nfgroup)
 #if defined(HAVE_LAPACKE)
       LAPACKE_dposv(LAPACK_COL_MAJOR, 'L',
 		ncoefftot, 1, alpha, ncoefftot, beta, ncoefftot);
-#else
+#elif defined(HAVE_ACCELERATE)
+      int info,n_beta=1;
+      char* upper_or_lower="U";
+      dposv_(upper_or_lower, &ncoefftot, &n_beta, alpha, &ncoefftot, beta, &ncoefftot, &info);
+ #else
       clapack_dposv(CblasRowMajor, CblasUpper,
 		ncoefftot, 1, alpha, ncoefftot, beta, ncoefftot);
 #endif

src/proper.c

@@ -57,6 +57,10 @@
 #include ATLAS_LAPACK_H
 #endif
 
+#ifdef HAVE_ACCELERATE
+#include ACCELERATE_H
+#endif
+
 #ifdef HAVE_LAPACKE
 #include LAPACKE_H
 #endif
@@ -383,6 +387,10 @@ void	astrprop_fgroup(fgroupstruct *fgroup)
       nfreemin = astrprop_solve(fgroup,samp,wcsec,alpha,beta,wis, &chi2min);
 #if defined(HAVE_LAPACKE)
     LAPACKE_dpotri(LAPACK_COL_MAJOR, 'L',ncoeff, alpha, ncoeff);
+#elif defined(HAVE_ACCELERATE)
+    int info,n_beta=1;
+    char* upper_or_lower="U";
+    dpotri_(upper_or_lower, &ncoeff, alpha, &ncoeff, &info);
 #else
     clapack_dpotri(CblasRowMajor, CblasUpper, ncoeff, alpha, ncoeff);
 #endif
@@ -605,6 +613,10 @@ static int	astrprop_solve(fgroupstruct *fgroup, samplestruct *samp,
   memcpy(a,alpha,ncoeff*ncoeff*sizeof(double));
 #if defined(HAVE_LAPACKE)
   LAPACKE_dposv(LAPACK_COL_MAJOR, 'L', ncoeff, 1, alpha, ncoeff, beta, ncoeff);
+#elif defined(HAVE_ACCELERATE)
+  int info,n_beta=1;
+  char* upper_or_lower="U";
+  dposv_(upper_or_lower, &ncoeff, &n_beta, alpha, &ncoeff, beta, &ncoeff, &info);
 #else
   clapack_dposv(CblasRowMajor,CblasUpper,ncoeff,1, alpha,ncoeff,beta,ncoeff);
 #endif