Move some linear algebra/matrix operations to internal API (#1165)

* Move unoptimized matrix operations to internal API

* Minor docstring formatting in ceed-basis.c

* Update Julia/Python/Fortran APIs for removed public functions

* Fix include for tests

* Fixes for removed Python API functions

Author: sebastiangrimberg

examples/python/tutorial-3-basis.ipynb

@@ -312,111 +312,6 @@
     "  # Check that (1' * G * u - u' * (G' * 1)) is numerically zero\n",
     "  print('1T * G * u - uT * (GT * 1) =', np.abs(sum1 - sum2))"
    ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Advanced topics"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "* In the following example, we demonstrate the QR factorization of a basis matrix.\n",
-    "The representation is similar to LAPACK's [`dgeqrf`](https://www.netlib.org/lapack/explore-html/dd/d9a/group__double_g_ecomputational_ga3766ea903391b5cf9008132f7440ec7b.html#ga3766ea903391b5cf9008132f7440ec7b), with elementary reflectors in the lower triangular block, scaled by `tau`."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "qr = np.array([1, -1, 4, 1, 4, -2, 1, 4, 2, 1, -1, 0], dtype=\"float64\")\n",
-    "tau = np.empty(3, dtype=\"float64\")\n",
-    "\n",
-    "qr, tau = ceed.qr_factorization(qr, tau, 4, 3)\n",
-    "\n",
-    "print('qr =')\n",
-    "print(qr.reshape(4, 3))\n",
-    "\n",
-    "print('tau =')\n",
-    "print(tau)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "* In the following example, we demonstrate the symmetric Schur decomposition of a basis matrix"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "A = np.array([0.19996678, 0.0745459, -0.07448852, 0.0332866,\n",
-    "              0.0745459, 1., 0.16666509, -0.07448852,\n",
-    "              -0.07448852, 0.16666509, 1., 0.0745459,\n",
-    "              0.0332866, -0.07448852, 0.0745459, 0.19996678], dtype=\"float64\")\n",
-    "\n",
-    "lam = ceed.symmetric_schur_decomposition(A, 4)\n",
-    "\n",
-    "print(\"Q =\")\n",
-    "for i in range(4):\n",
-    "  for j in range(4):\n",
-    "    if A[j+4*i] <= 1E-14 and A[j+4*i] >= -1E-14:\n",
-    "       A[j+4*i] = 0\n",
-    "    print(\"%12.8f\"%A[j+4*i])\n",
-    "\n",
-    "print(\"lambda =\")\n",
-    "for i in range(4):\n",
-    "  if lam[i] <= 1E-14 and lam[i] >= -1E-14:\n",
-    "    lam[i] = 0\n",
-    "  print(\"%12.8f\"%lam[i])"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "* In the following example, we demonstrate the simultaneous diagonalization of a basis matrix"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "M = np.array([0.19996678, 0.0745459, -0.07448852, 0.0332866,\n",
-    "              0.0745459, 1., 0.16666509, -0.07448852,\n",
-    "              -0.07448852, 0.16666509, 1., 0.0745459,\n",
-    "              0.0332866, -0.07448852, 0.0745459, 0.19996678], dtype=\"float64\")\n",
-    "K = np.array([3.03344425, -3.41501767, 0.49824435, -0.11667092,\n",
-    "              -3.41501767, 5.83354662, -2.9167733, 0.49824435,\n",
-    "              0.49824435, -2.9167733, 5.83354662, -3.41501767,\n",
-    "              -0.11667092, 0.49824435, -3.41501767, 3.03344425], dtype=\"float64\")\n",
-    "\n",
-    "x, lam = ceed.simultaneous_diagonalization(K, M, 4)\n",
-    "\n",
-    "print(\"x =\")\n",
-    "for i in range(4):\n",
-    "  for j in range(4):\n",
-    "    if x[j+4*i] <= 1E-14 and x[j+4*i] >= -1E-14:\n",
-    "      x[j+4*i] = 0\n",
-    "    print(\"%12.8f\"%x[j+4*i])\n",
-    "\n",
-    "print(\"lambda =\")\n",
-    "for i in range(4):\n",
-    "  if lam[i] <= 1E-14 and lam[i] >= -1E-14:\n",
-    "    lam[i] = 0\n",
-    "  print(\"%12.8f\"%lam[i])"
-   ]
   }
  ],
  "metadata": {

include/ceed/backend.h

@@ -181,8 +181,6 @@ typedef enum {
 CEED_EXTERN const char *const CeedFESpaces[];
 
 CEED_EXTERN int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *colo_grad_1d);
-CEED_EXTERN int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m, CeedInt n,
-                                      CeedInt k, CeedInt row, CeedInt col);
 CEED_EXTERN int CeedBasisIsTensor(CeedBasis basis, bool *is_tensor);
 CEED_EXTERN int CeedBasisGetData(CeedBasis basis, void *data);
 CEED_EXTERN int CeedBasisSetData(CeedBasis basis, void *data);
@@ -282,5 +280,10 @@ CEED_EXTERN int CeedOperatorSetSetupDone(CeedOperator op);
 
 CEED_INTERN int CeedMatrixMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, CeedInt n,
                                          CeedInt kk);
+CEED_EXTERN int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, CeedInt n);
+CEED_EXTERN int CeedHouseholderApplyQ(CeedScalar *mat_A, const CeedScalar *mat_Q, const CeedScalar *tau, CeedTransposeMode t_mode, CeedInt m,
+                                      CeedInt n, CeedInt k, CeedInt row, CeedInt col);
+CEED_EXTERN int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n);
+CEED_EXTERN int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *x, CeedScalar *lambda, CeedInt n);
 
 #endif

include/ceed/ceed.h

@@ -388,9 +388,6 @@ CEED_EXTERN int CeedBasisDestroy(CeedBasis *basis);
 
 CEED_EXTERN int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d);
 CEED_EXTERN int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, CeedScalar *q_weight_1d);
-CEED_EXTERN int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, CeedInt m, CeedInt n);
-CEED_EXTERN int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, CeedScalar *lambda, CeedInt n);
-CEED_EXTERN int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, CeedScalar *mat_B, CeedScalar *x, CeedScalar *lambda, CeedInt n);
 
 /** Handle for the user provided CeedQFunction callback function