Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 6 commits into from
Mar 17, 2023
Merged
Show file tree
Hide file tree
Changes from 5 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
105 changes: 0 additions & 105 deletions examples/python/tutorial-3-basis.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -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": {
Expand Down
7 changes: 5 additions & 2 deletions include/ceed/backend.h
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down Expand Up @@ -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
3 changes: 0 additions & 3 deletions include/ceed/ceed.h
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
Loading