For kicks and giggles, add an (untested) C+BLAS/LAPACK

implementation of the 'ip_simple' mimetic inner product.

  Suggested by: jrn.

blas_lapack.h

@@ -0,0 +1,31 @@
+// C <- a1*op(A)*op(B) + a2*C  where  op(X) in {X, X.'}
+void dgemm_(char *transA        , char *transB         ,
+            const MAT_SIZE_T*  m, const MAT_SIZE_T* n  , const MAT_SIZE_T* k  ,
+            const double*     a1, const double*     A  , const MAT_SIZE_T* ldA,
+            const double*      B, const MAT_SIZE_T* ldB,
+            const double*     a2,       double*     C  , const MAT_SIZE_T* ldC);
+
+// C <- a1*A*A' + a2*C   *or*   C <- a1*A'*A + a2*C
+void dsyrk_(char             *uplo, char             *trans,
+            const MAT_SIZE_T *n   , const MAT_SIZE_T *k    ,
+            const double     *a1  , const double     *A    , const MAT_SIZE_T *ldA,
+            const double     *a2  ,       double     *C    , const MAT_SIZE_T *ldC);
+
+
+// B <- a*op(A)*B  *or*  B <- a*B*op(A)  where op(X) \in {X, X.', X'}
+void dtrmm_(char *, char *, char *, char *,
+            const MAT_SIZE_T *m, const MAT_SIZE_T* n  ,
+            const double     *a,
+            const double     *A, const MAT_SIZE_T* ldA,
+                  double     *B, const MAT_SIZE_T* ldB);
+
+
+void dgeqrf_(const MAT_SIZE_T *m    , const MAT_SIZE_T *n   ,
+                   double     *A    , const MAT_SIZE_T *ld  ,
+                   double     *tau  ,       double     *work,
+             const MAT_SIZE_T *lwork,       MAT_SIZE_T *info);
+
+
+void dorgqr_(const MAT_SIZE_T *m   , const MAT_SIZE_T *n    , const MAT_SIZE_T *k  ,
+                   double     *A   , const MAT_SIZE_T *ld   , const double     *tau,
+                   double     *work, const MAT_SIZE_T *lwork,       MAT_SIZE_T *info);

mimetic.c

@@ -0,0 +1,111 @@
+#include <stddef.h>
+
+#define MAT_SIZE_T size_t
+
+#include "blas_lapack.h"
+#include "mimetic.h"
+
+
+/* ------------------------------------------------------------------ */
+static double
+trace(size_t n, double *A);
+/* ------------------------------------------------------------------ */
+
+
+/* ------------------------------------------------------------------ */
+static void
+mim_ip_simple_fill_null(size_t nf, size_t d, double *C, double *A,
+                        double *Binv, double *work, MAT_SIZE_T lwork);
+/* ------------------------------------------------------------------ */
+
+
+/* ------------------------------------------------------------------ */
+void
+mim_ip_simple(MAT_SIZE_T nf, MAT_SIZE_T d, double vol,
+              double *A, double *K, double *C, double *N,
+              double *Binv, double* work, MAT_SIZE_T lwork)
+/* ------------------------------------------------------------------ */
+{
+    MAT_SIZE_T m, n, k, ld1, ld2, ldBinv;
+    double     a1, a2;
+
+    mim_ip_simple_fill_null(nf, d, C, A, Binv, work, lwork);
+
+    /* Step 4) C <- N*K */
+    m   = nf ;  n   = nf ;  k = d;
+    ld1 = nf ;  ld2 = d  ;
+    a1  = 1.0;  a2  = 0.0;
+    dgemm_("No Transpose", "No Transpose", &m, &n, &k,
+           &a1, N, &ld1, K, &ld1, &a2, C, &ld1);
+
+    /*  Step 5) Binv <- (N*K*N' + t*A*(I-Q*Q')*A) / vol */
+    a1 = 1.0                /      vol ;
+    a2 = 6.0 * trace(nf, K) / (d * vol);
+    ldBinv = nf;
+    dgemm_("No Transpose", "Transpose", &m, &m, &n,
+           &a1, C, &ld1, N, &ld1, &a2, Binv, &ldBinv);
+}
+
+
+/* ------------------------------------------------------------------ */
+static void
+mim_ip_simple_fill_null(size_t nf, size_t d, double *C, double *A,
+                        double *Binv,
+                        double *work, MAT_SIZE_T lwork)
+/* ------------------------------------------------------------------ */
+{
+    MAT_SIZE_T m, n, k, ld, info;
+
+    size_t i, j;
+    double a1, a2;
+
+    double tau[3] = { 0.0 };  /* No more than 3 spatial dimensions */
+
+    /* Step 1) Binv <- I */
+    for (i = 0; i < nf*nf; i++) { Binv[i           ] = 0.0; }
+    for (i = 0; i < nf   ; i++) { Binv[i * (nf + 1)] = 1.0; }
+
+    /* Step 2) C <- orth(A * C) */
+    for (j = 0; j < d; j++) {
+        for (i = 0; i < nf; i++) {
+            C[i + j*nf] *= A[i];
+        }
+    }
+
+    m = nf;  n = d;  ld = nf;  k = d;
+    dgeqrf_(&m, &n,     C, &ld, tau, work, &lwork, &info);
+    dorgqr_(&m, &n, &k, C, &ld, tau, work, &lwork, &info);
+
+    /* Step 3) Binv <- A * (Binv - C*C') * A */
+    a1 = -1.0;  a2 = 1.0;
+    dsyrk_("Upper Triangular", "No Transpose",
+           &m, &n, &a1, C, &ld, &a2, Binv, &ld);
+    for (j = 0; j < nf; j++) {
+        for (i = 0; i <= j; i++) {
+            Binv[i + j*nf] *= A[i] * A[j];
+        }
+    }
+
+    /* Account for DSYRK only assigning upper triangular part. */
+    for (j = 0; j < nf; j++) {
+        for (i = j + 1; i < nf; i++) {
+            A[i + j*nf] = A[j + i*nf];
+        }
+    }
+}
+
+
+/* ------------------------------------------------------------------ */
+static double
+trace(size_t n, double *A)
+/* ------------------------------------------------------------------ */
+{
+    size_t i;
+    double t = 0.0;
+
+    for (i = 0; i < n; i++) {
+        t += A[i + i*n];
+    }
+
+    return t;
+}

mimetic.h

@@ -0,0 +1,10 @@
+#ifndef MIMETIC_H_INCLUDED
+#define MIMETIC_H_INCLUDED
+
+void mim_ip_simple(MAT_SIZE_T nf, MAT_SIZE_T d, double v,
+                   double *A, double *K, double *C, double *N,
+                   double *Binv,
+                   double *work, MAT_SIZE_T lwork);
+
+
+#endif /* MIMETIC_H_INCLUDED */