1 // Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at 2 // the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights 3 // reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed/ceed.h> 18 #include <ceed/backend.h> 19 #include <ceed-impl.h> 20 #include <math.h> 21 #include <stdbool.h> 22 #include <stdio.h> 23 #include <string.h> 24 25 /// @file 26 /// Implementation of CeedBasis interfaces 27 28 /// @cond DOXYGEN_SKIP 29 static struct CeedBasis_private ceed_basis_collocated; 30 /// @endcond 31 32 /// @addtogroup CeedBasisUser 33 /// @{ 34 35 /// Indicate that the quadrature points are collocated with the nodes 36 const CeedBasis CEED_BASIS_COLLOCATED = &ceed_basis_collocated; 37 38 /// @} 39 40 /// ---------------------------------------------------------------------------- 41 /// CeedBasis Library Internal Functions 42 /// ---------------------------------------------------------------------------- 43 /// @addtogroup CeedBasisDeveloper 44 /// @{ 45 46 /** 47 @brief Compute Householder reflection 48 49 Computes A = (I - b v v^T) A 50 where A is an mxn matrix indexed as A[i*row + j*col] 51 52 @param[in,out] A Matrix to apply Householder reflection to, in place 53 @param v Householder vector 54 @param b Scaling factor 55 @param m Number of rows in A 56 @param n Number of columns in A 57 @param row Row stride 58 @param col Col stride 59 60 @return An error code: 0 - success, otherwise - failure 61 62 @ref Developer 63 **/ 64 static int CeedHouseholderReflect(CeedScalar *A, const CeedScalar *v, 65 CeedScalar b, CeedInt m, CeedInt n, 66 CeedInt row, CeedInt col) { 67 for (CeedInt j=0; j<n; j++) { 68 CeedScalar w = A[0*row + j*col]; 69 for (CeedInt i=1; i<m; i++) 70 w += v[i] * A[i*row + j*col]; 71 A[0*row + j*col] -= b * w; 72 for (CeedInt i=1; i<m; i++) 73 A[i*row + j*col] -= b * w * v[i]; 74 } 75 return CEED_ERROR_SUCCESS; 76 } 77 78 /** 79 @brief Apply Householder Q matrix 80 81 Compute A = Q A where Q is mxm and A is mxn. 82 83 @param[in,out] A Matrix to apply Householder Q to, in place 84 @param Q Householder Q matrix 85 @param tau Householder scaling factors 86 @param t_mode Transpose mode for application 87 @param m Number of rows in A 88 @param n Number of columns in A 89 @param k Number of elementary reflectors in Q, k<m 90 @param row Row stride in A 91 @param col Col stride in A 92 93 @return An error code: 0 - success, otherwise - failure 94 95 @ref Developer 96 **/ 97 int CeedHouseholderApplyQ(CeedScalar *A, const CeedScalar *Q, 98 const CeedScalar *tau, CeedTransposeMode t_mode, 99 CeedInt m, CeedInt n, CeedInt k, 100 CeedInt row, CeedInt col) { 101 int ierr; 102 CeedScalar *v; 103 ierr = CeedMalloc(m, &v); CeedChk(ierr); 104 for (CeedInt ii=0; ii<k; ii++) { 105 CeedInt i = t_mode == CEED_TRANSPOSE ? ii : k-1-ii; 106 for (CeedInt j=i+1; j<m; j++) 107 v[j] = Q[j*k+i]; 108 // Apply Householder reflector (I - tau v v^T) collo_grad_1d^T 109 ierr = CeedHouseholderReflect(&A[i*row], &v[i], tau[i], m-i, n, row, col); 110 CeedChk(ierr); 111 } 112 ierr = CeedFree(&v); CeedChk(ierr); 113 return CEED_ERROR_SUCCESS; 114 } 115 116 /** 117 @brief Compute Givens rotation 118 119 Computes A = G A (or G^T A in transpose mode) 120 where A is an mxn matrix indexed as A[i*n + j*m] 121 122 @param[in,out] A Row major matrix to apply Givens rotation to, in place 123 @param c Cosine factor 124 @param s Sine factor 125 @param t_mode @ref CEED_NOTRANSPOSE to rotate the basis counter-clockwise, 126 which has the effect of rotating columns of A clockwise; 127 @ref CEED_TRANSPOSE for the opposite rotation 128 @param i First row/column to apply rotation 129 @param k Second row/column to apply rotation 130 @param m Number of rows in A 131 @param n Number of columns in A 132 133 @return An error code: 0 - success, otherwise - failure 134 135 @ref Developer 136 **/ 137 static int CeedGivensRotation(CeedScalar *A, CeedScalar c, CeedScalar s, 138 CeedTransposeMode t_mode, CeedInt i, CeedInt k, 139 CeedInt m, CeedInt n) { 140 CeedInt stride_j = 1, stride_ik = m, num_its = n; 141 if (t_mode == CEED_NOTRANSPOSE) { 142 stride_j = n; stride_ik = 1; num_its = m; 143 } 144 145 // Apply rotation 146 for (CeedInt j=0; j<num_its; j++) { 147 CeedScalar tau1 = A[i*stride_ik+j*stride_j], tau2 = A[k*stride_ik+j*stride_j]; 148 A[i*stride_ik+j*stride_j] = c*tau1 - s*tau2; 149 A[k*stride_ik+j*stride_j] = s*tau1 + c*tau2; 150 } 151 return CEED_ERROR_SUCCESS; 152 } 153 154 /** 155 @brief View an array stored in a CeedBasis 156 157 @param[in] name Name of array 158 @param[in] fp_fmt Printing format 159 @param[in] m Number of rows in array 160 @param[in] n Number of columns in array 161 @param[in] a Array to be viewed 162 @param[in] stream Stream to view to, e.g., stdout 163 164 @return An error code: 0 - success, otherwise - failure 165 166 @ref Developer 167 **/ 168 static int CeedScalarView(const char *name, const char *fp_fmt, CeedInt m, 169 CeedInt n, const CeedScalar *a, FILE *stream) { 170 for (int i=0; i<m; i++) { 171 if (m > 1) 172 fprintf(stream, "%12s[%d]:", name, i); 173 else 174 fprintf(stream, "%12s:", name); 175 for (int j=0; j<n; j++) 176 fprintf(stream, fp_fmt, fabs(a[i*n+j]) > 1E-14 ? a[i*n+j] : 0); 177 fputs("\n", stream); 178 } 179 return CEED_ERROR_SUCCESS; 180 } 181 182 /// @} 183 184 /// ---------------------------------------------------------------------------- 185 /// Ceed Backend API 186 /// ---------------------------------------------------------------------------- 187 /// @addtogroup CeedBasisBackend 188 /// @{ 189 190 /** 191 @brief Return collocated grad matrix 192 193 @param basis CeedBasis 194 @param[out] collo_grad_1d Row-major (Q_1d * Q_1d) matrix expressing derivatives of 195 basis functions at quadrature points 196 197 @return An error code: 0 - success, otherwise - failure 198 199 @ref Backend 200 **/ 201 int CeedBasisGetCollocatedGrad(CeedBasis basis, CeedScalar *collo_grad_1d) { 202 int i, j, k; 203 Ceed ceed; 204 CeedInt ierr, P_1d=(basis)->P_1d, Q_1d=(basis)->Q_1d; 205 CeedScalar *interp_1d, *grad_1d, *tau; 206 207 ierr = CeedMalloc(Q_1d*P_1d, &interp_1d); CeedChk(ierr); 208 ierr = CeedMalloc(Q_1d*P_1d, &grad_1d); CeedChk(ierr); 209 ierr = CeedMalloc(Q_1d, &tau); CeedChk(ierr); 210 memcpy(interp_1d, (basis)->interp_1d, Q_1d*P_1d*sizeof(basis)->interp_1d[0]); 211 memcpy(grad_1d, (basis)->grad_1d, Q_1d*P_1d*sizeof(basis)->interp_1d[0]); 212 213 // QR Factorization, interp_1d = Q R 214 ierr = CeedBasisGetCeed(basis, &ceed); CeedChk(ierr); 215 ierr = CeedQRFactorization(ceed, interp_1d, tau, Q_1d, P_1d); CeedChk(ierr); 216 // Note: This function is for backend use, so all errors are terminal 217 // and we do not need to clean up memory on failure. 218 219 // Apply Rinv, collo_grad_1d = grad_1d Rinv 220 for (i=0; i<Q_1d; i++) { // Row i 221 collo_grad_1d[Q_1d*i] = grad_1d[P_1d*i]/interp_1d[0]; 222 for (j=1; j<P_1d; j++) { // Column j 223 collo_grad_1d[j+Q_1d*i] = grad_1d[j+P_1d*i]; 224 for (k=0; k<j; k++) 225 collo_grad_1d[j+Q_1d*i] -= interp_1d[j+P_1d*k]*collo_grad_1d[k+Q_1d*i]; 226 collo_grad_1d[j+Q_1d*i] /= interp_1d[j+P_1d*j]; 227 } 228 for (j=P_1d; j<Q_1d; j++) 229 collo_grad_1d[j+Q_1d*i] = 0; 230 } 231 232 // Apply Qtranspose, collo_grad = collo_grad Q_transpose 233 ierr = CeedHouseholderApplyQ(collo_grad_1d, interp_1d, tau, CEED_NOTRANSPOSE, 234 Q_1d, Q_1d, P_1d, 1, Q_1d); CeedChk(ierr); 235 236 ierr = CeedFree(&interp_1d); CeedChk(ierr); 237 ierr = CeedFree(&grad_1d); CeedChk(ierr); 238 ierr = CeedFree(&tau); CeedChk(ierr); 239 return CEED_ERROR_SUCCESS; 240 } 241 242 /** 243 @brief Get tensor status for given CeedBasis 244 245 @param basis CeedBasis 246 @param[out] is_tensor Variable to store tensor status 247 248 @return An error code: 0 - success, otherwise - failure 249 250 @ref Backend 251 **/ 252 int CeedBasisIsTensor(CeedBasis basis, bool *is_tensor) { 253 *is_tensor = basis->tensor_basis; 254 return CEED_ERROR_SUCCESS; 255 } 256 257 /** 258 @brief Get backend data of a CeedBasis 259 260 @param basis CeedBasis 261 @param[out] data Variable to store data 262 263 @return An error code: 0 - success, otherwise - failure 264 265 @ref Backend 266 **/ 267 int CeedBasisGetData(CeedBasis basis, void *data) { 268 *(void **)data = basis->data; 269 return CEED_ERROR_SUCCESS; 270 } 271 272 /** 273 @brief Set backend data of a CeedBasis 274 275 @param[out] basis CeedBasis 276 @param data Data to set 277 278 @return An error code: 0 - success, otherwise - failure 279 280 @ref Backend 281 **/ 282 int CeedBasisSetData(CeedBasis basis, void *data) { 283 basis->data = data; 284 return CEED_ERROR_SUCCESS; 285 } 286 287 /** 288 @brief Increment the reference counter for a CeedBasis 289 290 @param basis Basis to increment the reference counter 291 292 @return An error code: 0 - success, otherwise - failure 293 294 @ref Backend 295 **/ 296 int CeedBasisReference(CeedBasis basis) { 297 basis->ref_count++; 298 return CEED_ERROR_SUCCESS; 299 } 300 301 /** 302 @brief Get dimension for given CeedElemTopology 303 304 @param topo CeedElemTopology 305 @param[out] dim Variable to store dimension of topology 306 307 @return An error code: 0 - success, otherwise - failure 308 309 @ref Backend 310 **/ 311 int CeedBasisGetTopologyDimension(CeedElemTopology topo, CeedInt *dim) { 312 *dim = (CeedInt) topo >> 16; 313 return CEED_ERROR_SUCCESS; 314 } 315 316 /** 317 @brief Get CeedTensorContract of a CeedBasis 318 319 @param basis CeedBasis 320 @param[out] contract Variable to store CeedTensorContract 321 322 @return An error code: 0 - success, otherwise - failure 323 324 @ref Backend 325 **/ 326 int CeedBasisGetTensorContract(CeedBasis basis, CeedTensorContract *contract) { 327 *contract = basis->contract; 328 return CEED_ERROR_SUCCESS; 329 } 330 331 /** 332 @brief Set CeedTensorContract of a CeedBasis 333 334 @param[out] basis CeedBasis 335 @param contract CeedTensorContract to set 336 337 @return An error code: 0 - success, otherwise - failure 338 339 @ref Backend 340 **/ 341 int CeedBasisSetTensorContract(CeedBasis basis, CeedTensorContract contract) { 342 int ierr; 343 basis->contract = contract; 344 ierr = CeedTensorContractReference(contract); CeedChk(ierr); 345 return CEED_ERROR_SUCCESS; 346 } 347 348 /** 349 @brief Return a reference implementation of matrix multiplication C = A B. 350 Note, this is a reference implementation for CPU CeedScalar pointers 351 that is not intended for high performance. 352 353 @param ceed A Ceed context for error handling 354 @param[in] mat_A Row-major matrix A 355 @param[in] mat_B Row-major matrix B 356 @param[out] mat_C Row-major output matrix C 357 @param m Number of rows of C 358 @param n Number of columns of C 359 @param kk Number of columns of A/rows of B 360 361 @return An error code: 0 - success, otherwise - failure 362 363 @ref Utility 364 **/ 365 int CeedMatrixMultiply(Ceed ceed, const CeedScalar *mat_A, 366 const CeedScalar *mat_B, CeedScalar *mat_C, CeedInt m, 367 CeedInt n, CeedInt kk) { 368 for (CeedInt i=0; i<m; i++) 369 for (CeedInt j=0; j<n; j++) { 370 CeedScalar sum = 0; 371 for (CeedInt k=0; k<kk; k++) 372 sum += mat_A[k+i*kk]*mat_B[j+k*n]; 373 mat_C[j+i*n] = sum; 374 } 375 return CEED_ERROR_SUCCESS; 376 } 377 378 /// @} 379 380 /// ---------------------------------------------------------------------------- 381 /// CeedBasis Public API 382 /// ---------------------------------------------------------------------------- 383 /// @addtogroup CeedBasisUser 384 /// @{ 385 386 /** 387 @brief Create a tensor-product basis for H^1 discretizations 388 389 @param ceed A Ceed object where the CeedBasis will be created 390 @param dim Topological dimension 391 @param num_comp Number of field components (1 for scalar fields) 392 @param P_1d Number of nodes in one dimension 393 @param Q_1d Number of quadrature points in one dimension 394 @param interp_1d Row-major (Q_1d * P_1d) matrix expressing the values of nodal 395 basis functions at quadrature points 396 @param grad_1d Row-major (Q_1d * P_1d) matrix expressing derivatives of nodal 397 basis functions at quadrature points 398 @param q_ref_1d Array of length Q_1d holding the locations of quadrature points 399 on the 1D reference element [-1, 1] 400 @param q_weight_1d Array of length Q_1d holding the quadrature weights on the 401 reference element 402 @param[out] basis Address of the variable where the newly created 403 CeedBasis will be stored. 404 405 @return An error code: 0 - success, otherwise - failure 406 407 @ref User 408 **/ 409 int CeedBasisCreateTensorH1(Ceed ceed, CeedInt dim, CeedInt num_comp, 410 CeedInt P_1d, CeedInt Q_1d, 411 const CeedScalar *interp_1d, 412 const CeedScalar *grad_1d, const CeedScalar *q_ref_1d, 413 const CeedScalar *q_weight_1d, CeedBasis *basis) { 414 int ierr; 415 416 if (!ceed->BasisCreateTensorH1) { 417 Ceed delegate; 418 ierr = CeedGetObjectDelegate(ceed, &delegate, "Basis"); CeedChk(ierr); 419 420 if (!delegate) 421 // LCOV_EXCL_START 422 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 423 "Backend does not support BasisCreateTensorH1"); 424 // LCOV_EXCL_STOP 425 426 ierr = CeedBasisCreateTensorH1(delegate, dim, num_comp, P_1d, 427 Q_1d, interp_1d, grad_1d, q_ref_1d, 428 q_weight_1d, basis); CeedChk(ierr); 429 return CEED_ERROR_SUCCESS; 430 } 431 432 if (dim<1) 433 // LCOV_EXCL_START 434 return CeedError(ceed, CEED_ERROR_DIMENSION, 435 "Basis dimension must be a positive value"); 436 // LCOV_EXCL_STOP 437 CeedElemTopology topo = dim == 1 ? CEED_LINE 438 : dim == 2 ? CEED_QUAD 439 : CEED_HEX; 440 441 ierr = CeedCalloc(1, basis); CeedChk(ierr); 442 (*basis)->ceed = ceed; 443 ierr = CeedReference(ceed); CeedChk(ierr); 444 (*basis)->ref_count = 1; 445 (*basis)->tensor_basis = 1; 446 (*basis)->dim = dim; 447 (*basis)->topo = topo; 448 (*basis)->num_comp = num_comp; 449 (*basis)->P_1d = P_1d; 450 (*basis)->Q_1d = Q_1d; 451 (*basis)->P = CeedIntPow(P_1d, dim); 452 (*basis)->Q = CeedIntPow(Q_1d, dim); 453 ierr = CeedMalloc(Q_1d,&(*basis)->q_ref_1d); CeedChk(ierr); 454 ierr = CeedMalloc(Q_1d,&(*basis)->q_weight_1d); CeedChk(ierr); 455 memcpy((*basis)->q_ref_1d, q_ref_1d, Q_1d*sizeof(q_ref_1d[0])); 456 memcpy((*basis)->q_weight_1d, q_weight_1d, Q_1d*sizeof(q_weight_1d[0])); 457 ierr = CeedMalloc(Q_1d*P_1d,&(*basis)->interp_1d); CeedChk(ierr); 458 ierr = CeedMalloc(Q_1d*P_1d,&(*basis)->grad_1d); CeedChk(ierr); 459 memcpy((*basis)->interp_1d, interp_1d, Q_1d*P_1d*sizeof(interp_1d[0])); 460 memcpy((*basis)->grad_1d, grad_1d, Q_1d*P_1d*sizeof(grad_1d[0])); 461 ierr = ceed->BasisCreateTensorH1(dim, P_1d, Q_1d, interp_1d, grad_1d, q_ref_1d, 462 q_weight_1d, *basis); CeedChk(ierr); 463 return CEED_ERROR_SUCCESS; 464 } 465 466 /** 467 @brief Create a tensor-product Lagrange basis 468 469 @param ceed A Ceed object where the CeedBasis will be created 470 @param dim Topological dimension of element 471 @param num_comp Number of field components (1 for scalar fields) 472 @param P Number of Gauss-Lobatto nodes in one dimension. The 473 polynomial degree of the resulting Q_k element is k=P-1. 474 @param Q Number of quadrature points in one dimension. 475 @param quad_mode Distribution of the Q quadrature points (affects order of 476 accuracy for the quadrature) 477 @param[out] basis Address of the variable where the newly created 478 CeedBasis will be stored. 479 480 @return An error code: 0 - success, otherwise - failure 481 482 @ref User 483 **/ 484 int CeedBasisCreateTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt num_comp, 485 CeedInt P, CeedInt Q, CeedQuadMode quad_mode, 486 CeedBasis *basis) { 487 // Allocate 488 int ierr, ierr2, i, j, k; 489 CeedScalar c1, c2, c3, c4, dx, *nodes, *interp_1d, *grad_1d, *q_ref_1d, 490 *q_weight_1d; 491 492 if (dim<1) 493 // LCOV_EXCL_START 494 return CeedError(ceed, CEED_ERROR_DIMENSION, 495 "Basis dimension must be a positive value"); 496 // LCOV_EXCL_STOP 497 498 // Get Nodes and Weights 499 ierr = CeedCalloc(P*Q, &interp_1d); CeedChk(ierr); 500 ierr = CeedCalloc(P*Q, &grad_1d); CeedChk(ierr); 501 ierr = CeedCalloc(P, &nodes); CeedChk(ierr); 502 ierr = CeedCalloc(Q, &q_ref_1d); CeedChk(ierr); 503 ierr = CeedCalloc(Q, &q_weight_1d); CeedChk(ierr); 504 ierr = CeedLobattoQuadrature(P, nodes, NULL); 505 if (ierr) { goto cleanup; } CeedChk(ierr); 506 switch (quad_mode) { 507 case CEED_GAUSS: 508 ierr = CeedGaussQuadrature(Q, q_ref_1d, q_weight_1d); 509 break; 510 case CEED_GAUSS_LOBATTO: 511 ierr = CeedLobattoQuadrature(Q, q_ref_1d, q_weight_1d); 512 break; 513 } 514 if (ierr) { goto cleanup; } CeedChk(ierr); 515 516 // Build B, D matrix 517 // Fornberg, 1998 518 for (i = 0; i < Q; i++) { 519 c1 = 1.0; 520 c3 = nodes[0] - q_ref_1d[i]; 521 interp_1d[i*P+0] = 1.0; 522 for (j = 1; j < P; j++) { 523 c2 = 1.0; 524 c4 = c3; 525 c3 = nodes[j] - q_ref_1d[i]; 526 for (k = 0; k < j; k++) { 527 dx = nodes[j] - nodes[k]; 528 c2 *= dx; 529 if (k == j - 1) { 530 grad_1d[i*P + j] = c1*(interp_1d[i*P + k] - c4*grad_1d[i*P + k]) / c2; 531 interp_1d[i*P + j] = - c1*c4*interp_1d[i*P + k] / c2; 532 } 533 grad_1d[i*P + k] = (c3*grad_1d[i*P + k] - interp_1d[i*P + k]) / dx; 534 interp_1d[i*P + k] = c3*interp_1d[i*P + k] / dx; 535 } 536 c1 = c2; 537 } 538 } 539 // Pass to CeedBasisCreateTensorH1 540 ierr = CeedBasisCreateTensorH1(ceed, dim, num_comp, P, Q, interp_1d, grad_1d, 541 q_ref_1d, q_weight_1d, basis); CeedChk(ierr); 542 cleanup: 543 ierr2 = CeedFree(&interp_1d); CeedChk(ierr2); 544 ierr2 = CeedFree(&grad_1d); CeedChk(ierr2); 545 ierr2 = CeedFree(&nodes); CeedChk(ierr2); 546 ierr2 = CeedFree(&q_ref_1d); CeedChk(ierr2); 547 ierr2 = CeedFree(&q_weight_1d); CeedChk(ierr2); 548 CeedChk(ierr); 549 return CEED_ERROR_SUCCESS; 550 } 551 552 /** 553 @brief Create a non tensor-product basis for H^1 discretizations 554 555 @param ceed A Ceed object where the CeedBasis will be created 556 @param topo Topology of element, e.g. hypercube, simplex, ect 557 @param num_comp Number of field components (1 for scalar fields) 558 @param num_nodes Total number of nodes 559 @param num_qpts Total number of quadrature points 560 @param interp Row-major (num_qpts * num_nodes) matrix expressing the values of 561 nodal basis functions at quadrature points 562 @param grad Row-major (num_qpts * dim * num_nodes) matrix expressing 563 derivatives of nodal basis functions at quadrature points 564 @param q_ref Array of length num_qpts holding the locations of quadrature 565 points on the reference element [-1, 1] 566 @param q_weight Array of length num_qpts holding the quadrature weights on the 567 reference element 568 @param[out] basis Address of the variable where the newly created 569 CeedBasis will be stored. 570 571 @return An error code: 0 - success, otherwise - failure 572 573 @ref User 574 **/ 575 int CeedBasisCreateH1(Ceed ceed, CeedElemTopology topo, CeedInt num_comp, 576 CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp, 577 const CeedScalar *grad, const CeedScalar *q_ref, 578 const CeedScalar *q_weight, CeedBasis *basis) { 579 int ierr; 580 CeedInt P = num_nodes, Q = num_qpts, dim = 0; 581 582 if (!ceed->BasisCreateH1) { 583 Ceed delegate; 584 ierr = CeedGetObjectDelegate(ceed, &delegate, "Basis"); CeedChk(ierr); 585 586 if (!delegate) 587 // LCOV_EXCL_START 588 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 589 "Backend does not support BasisCreateH1"); 590 // LCOV_EXCL_STOP 591 592 ierr = CeedBasisCreateH1(delegate, topo, num_comp, num_nodes, 593 num_qpts, interp, grad, q_ref, 594 q_weight, basis); CeedChk(ierr); 595 return CEED_ERROR_SUCCESS; 596 } 597 598 ierr = CeedCalloc(1,basis); CeedChk(ierr); 599 600 ierr = CeedBasisGetTopologyDimension(topo, &dim); CeedChk(ierr); 601 602 (*basis)->ceed = ceed; 603 ierr = CeedReference(ceed); CeedChk(ierr); 604 (*basis)->ref_count = 1; 605 (*basis)->tensor_basis = 0; 606 (*basis)->dim = dim; 607 (*basis)->topo = topo; 608 (*basis)->num_comp = num_comp; 609 (*basis)->P = P; 610 (*basis)->Q = Q; 611 ierr = CeedMalloc(Q*dim,&(*basis)->q_ref_1d); CeedChk(ierr); 612 ierr = CeedMalloc(Q,&(*basis)->q_weight_1d); CeedChk(ierr); 613 memcpy((*basis)->q_ref_1d, q_ref, Q*dim*sizeof(q_ref[0])); 614 memcpy((*basis)->q_weight_1d, q_weight, Q*sizeof(q_weight[0])); 615 ierr = CeedMalloc(Q*P, &(*basis)->interp); CeedChk(ierr); 616 ierr = CeedMalloc(dim*Q*P, &(*basis)->grad); CeedChk(ierr); 617 memcpy((*basis)->interp, interp, Q*P*sizeof(interp[0])); 618 memcpy((*basis)->grad, grad, dim*Q*P*sizeof(grad[0])); 619 ierr = ceed->BasisCreateH1(topo, dim, P, Q, interp, grad, q_ref, 620 q_weight, *basis); CeedChk(ierr); 621 return CEED_ERROR_SUCCESS; 622 } 623 624 /** 625 @brief Copy the pointer to a CeedBasis. Both pointers should 626 be destroyed with `CeedBasisDestroy()`; 627 Note: If `*basis_copy` is non-NULL, then it is assumed that 628 `*basis_copy` is a pointer to a CeedBasis. This CeedBasis 629 will be destroyed if `*basis_copy` is the only 630 reference to this CeedBasis. 631 632 @param basis CeedBasis to copy reference to 633 @param[out] basis_copy Variable to store copied reference 634 635 @return An error code: 0 - success, otherwise - failure 636 637 @ref User 638 **/ 639 int CeedBasisReferenceCopy(CeedBasis basis, CeedBasis *basis_copy) { 640 int ierr; 641 642 ierr = CeedBasisReference(basis); CeedChk(ierr); 643 ierr = CeedBasisDestroy(basis_copy); CeedChk(ierr); 644 *basis_copy = basis; 645 return CEED_ERROR_SUCCESS; 646 } 647 648 /** 649 @brief View a CeedBasis 650 651 @param basis CeedBasis to view 652 @param stream Stream to view to, e.g., stdout 653 654 @return An error code: 0 - success, otherwise - failure 655 656 @ref User 657 **/ 658 int CeedBasisView(CeedBasis basis, FILE *stream) { 659 int ierr; 660 661 if (basis->tensor_basis) { 662 fprintf(stream, "CeedBasis: dim=%d P=%d Q=%d\n", basis->dim, basis->P_1d, 663 basis->Q_1d); 664 ierr = CeedScalarView("qref1d", "\t% 12.8f", 1, basis->Q_1d, basis->q_ref_1d, 665 stream); CeedChk(ierr); 666 ierr = CeedScalarView("qweight1d", "\t% 12.8f", 1, basis->Q_1d, 667 basis->q_weight_1d, stream); CeedChk(ierr); 668 ierr = CeedScalarView("interp1d", "\t% 12.8f", basis->Q_1d, basis->P_1d, 669 basis->interp_1d, stream); CeedChk(ierr); 670 ierr = CeedScalarView("grad1d", "\t% 12.8f", basis->Q_1d, basis->P_1d, 671 basis->grad_1d, stream); CeedChk(ierr); 672 } else { 673 fprintf(stream, "CeedBasis: dim=%d P=%d Q=%d\n", basis->dim, basis->P, 674 basis->Q); 675 ierr = CeedScalarView("qref", "\t% 12.8f", 1, basis->Q*basis->dim, 676 basis->q_ref_1d, 677 stream); CeedChk(ierr); 678 ierr = CeedScalarView("qweight", "\t% 12.8f", 1, basis->Q, basis->q_weight_1d, 679 stream); CeedChk(ierr); 680 ierr = CeedScalarView("interp", "\t% 12.8f", basis->Q, basis->P, 681 basis->interp, stream); CeedChk(ierr); 682 ierr = CeedScalarView("grad", "\t% 12.8f", basis->dim*basis->Q, basis->P, 683 basis->grad, stream); CeedChk(ierr); 684 } 685 return CEED_ERROR_SUCCESS; 686 } 687 688 /** 689 @brief Apply basis evaluation from nodes to quadrature points or vice versa 690 691 @param basis CeedBasis to evaluate 692 @param num_elem The number of elements to apply the basis evaluation to; 693 the backend will specify the ordering in 694 CeedElemRestrictionCreateBlocked() 695 @param t_mode \ref CEED_NOTRANSPOSE to evaluate from nodes to quadrature 696 points, \ref CEED_TRANSPOSE to apply the transpose, mapping 697 from quadrature points to nodes 698 @param eval_mode \ref CEED_EVAL_NONE to use values directly, 699 \ref CEED_EVAL_INTERP to use interpolated values, 700 \ref CEED_EVAL_GRAD to use gradients, 701 \ref CEED_EVAL_WEIGHT to use quadrature weights. 702 @param[in] u Input CeedVector 703 @param[out] v Output CeedVector 704 705 @return An error code: 0 - success, otherwise - failure 706 707 @ref User 708 **/ 709 int CeedBasisApply(CeedBasis basis, CeedInt num_elem, CeedTransposeMode t_mode, 710 CeedEvalMode eval_mode, CeedVector u, CeedVector v) { 711 int ierr; 712 CeedInt u_length = 0, v_length, dim, num_comp, num_nodes, num_qpts; 713 ierr = CeedBasisGetDimension(basis, &dim); CeedChk(ierr); 714 ierr = CeedBasisGetNumComponents(basis, &num_comp); CeedChk(ierr); 715 ierr = CeedBasisGetNumNodes(basis, &num_nodes); CeedChk(ierr); 716 ierr = CeedBasisGetNumQuadraturePoints(basis, &num_qpts); CeedChk(ierr); 717 ierr = CeedVectorGetLength(v, &v_length); CeedChk(ierr); 718 if (u) { 719 ierr = CeedVectorGetLength(u, &u_length); CeedChk(ierr); 720 } 721 722 if (!basis->Apply) 723 // LCOV_EXCL_START 724 return CeedError(basis->ceed, CEED_ERROR_UNSUPPORTED, 725 "Backend does not support BasisApply"); 726 // LCOV_EXCL_STOP 727 728 // Check compatibility of topological and geometrical dimensions 729 if ((t_mode == CEED_TRANSPOSE && (v_length%num_nodes != 0 || 730 u_length%num_qpts != 0)) || 731 (t_mode == CEED_NOTRANSPOSE && (u_length%num_nodes != 0 || 732 v_length%num_qpts != 0))) 733 // LCOV_EXCL_START 734 return CeedError(basis->ceed, CEED_ERROR_DIMENSION, 735 "Length of input/output vectors " 736 "incompatible with basis dimensions"); 737 // LCOV_EXCL_STOP 738 739 // Check vector lengths to prevent out of bounds issues 740 bool bad_dims = false; 741 switch (eval_mode) { 742 case CEED_EVAL_NONE: 743 case CEED_EVAL_INTERP: bad_dims = 744 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 745 v_length < num_elem*num_comp*num_nodes)) || 746 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 747 u_length < num_elem*num_comp*num_nodes))); 748 break; 749 case CEED_EVAL_GRAD: bad_dims = 750 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts*dim || 751 v_length < num_elem*num_comp*num_nodes)) || 752 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp*dim || 753 u_length < num_elem*num_comp*num_nodes))); 754 break; 755 case CEED_EVAL_WEIGHT: 756 bad_dims = v_length < num_elem*num_qpts; 757 break; 758 // LCOV_EXCL_START 759 case CEED_EVAL_DIV: bad_dims = 760 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 761 v_length < num_elem*num_comp*num_nodes)) || 762 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 763 u_length < num_elem*num_comp*num_nodes))); 764 break; 765 case CEED_EVAL_CURL: bad_dims = 766 ((t_mode == CEED_TRANSPOSE && (u_length < num_elem*num_comp*num_qpts || 767 v_length < num_elem*num_comp*num_nodes)) || 768 (t_mode == CEED_NOTRANSPOSE && (v_length < num_elem*num_qpts*num_comp || 769 u_length < num_elem*num_comp*num_nodes))); 770 break; 771 // LCOV_EXCL_STOP 772 } 773 if (bad_dims) 774 // LCOV_EXCL_START 775 return CeedError(basis->ceed, CEED_ERROR_DIMENSION, 776 "Input/output vectors too short for basis and evaluation mode"); 777 // LCOV_EXCL_STOP 778 779 ierr = basis->Apply(basis, num_elem, t_mode, eval_mode, u, v); CeedChk(ierr); 780 return CEED_ERROR_SUCCESS; 781 } 782 783 /** 784 @brief Get Ceed associated with a CeedBasis 785 786 @param basis CeedBasis 787 @param[out] ceed Variable to store Ceed 788 789 @return An error code: 0 - success, otherwise - failure 790 791 @ref Advanced 792 **/ 793 int CeedBasisGetCeed(CeedBasis basis, Ceed *ceed) { 794 *ceed = basis->ceed; 795 return CEED_ERROR_SUCCESS; 796 } 797 798 /** 799 @brief Get dimension for given CeedBasis 800 801 @param basis CeedBasis 802 @param[out] dim Variable to store dimension of basis 803 804 @return An error code: 0 - success, otherwise - failure 805 806 @ref Advanced 807 **/ 808 int CeedBasisGetDimension(CeedBasis basis, CeedInt *dim) { 809 *dim = basis->dim; 810 return CEED_ERROR_SUCCESS; 811 } 812 813 /** 814 @brief Get topology for given CeedBasis 815 816 @param basis CeedBasis 817 @param[out] topo Variable to store topology of basis 818 819 @return An error code: 0 - success, otherwise - failure 820 821 @ref Advanced 822 **/ 823 int CeedBasisGetTopology(CeedBasis basis, CeedElemTopology *topo) { 824 *topo = basis->topo; 825 return CEED_ERROR_SUCCESS; 826 } 827 828 /** 829 @brief Get number of components for given CeedBasis 830 831 @param basis CeedBasis 832 @param[out] num_comp Variable to store number of components of basis 833 834 @return An error code: 0 - success, otherwise - failure 835 836 @ref Advanced 837 **/ 838 int CeedBasisGetNumComponents(CeedBasis basis, CeedInt *num_comp) { 839 *num_comp = basis->num_comp; 840 return CEED_ERROR_SUCCESS; 841 } 842 843 /** 844 @brief Get total number of nodes (in dim dimensions) of a CeedBasis 845 846 @param basis CeedBasis 847 @param[out] P Variable to store number of nodes 848 849 @return An error code: 0 - success, otherwise - failure 850 851 @ref Utility 852 **/ 853 int CeedBasisGetNumNodes(CeedBasis basis, CeedInt *P) { 854 *P = basis->P; 855 return CEED_ERROR_SUCCESS; 856 } 857 858 /** 859 @brief Get total number of nodes (in 1 dimension) of a CeedBasis 860 861 @param basis CeedBasis 862 @param[out] P_1d Variable to store number of nodes 863 864 @return An error code: 0 - success, otherwise - failure 865 866 @ref Advanced 867 **/ 868 int CeedBasisGetNumNodes1D(CeedBasis basis, CeedInt *P_1d) { 869 if (!basis->tensor_basis) 870 // LCOV_EXCL_START 871 return CeedError(basis->ceed, CEED_ERROR_MINOR, 872 "Cannot supply P_1d for non-tensor basis"); 873 // LCOV_EXCL_STOP 874 875 *P_1d = basis->P_1d; 876 return CEED_ERROR_SUCCESS; 877 } 878 879 /** 880 @brief Get total number of quadrature points (in dim dimensions) of a CeedBasis 881 882 @param basis CeedBasis 883 @param[out] Q Variable to store number of quadrature points 884 885 @return An error code: 0 - success, otherwise - failure 886 887 @ref Utility 888 **/ 889 int CeedBasisGetNumQuadraturePoints(CeedBasis basis, CeedInt *Q) { 890 *Q = basis->Q; 891 return CEED_ERROR_SUCCESS; 892 } 893 894 /** 895 @brief Get total number of quadrature points (in 1 dimension) of a CeedBasis 896 897 @param basis CeedBasis 898 @param[out] Q_1d Variable to store number of quadrature points 899 900 @return An error code: 0 - success, otherwise - failure 901 902 @ref Advanced 903 **/ 904 int CeedBasisGetNumQuadraturePoints1D(CeedBasis basis, CeedInt *Q_1d) { 905 if (!basis->tensor_basis) 906 // LCOV_EXCL_START 907 return CeedError(basis->ceed, CEED_ERROR_MINOR, 908 "Cannot supply Q_1d for non-tensor basis"); 909 // LCOV_EXCL_STOP 910 911 *Q_1d = basis->Q_1d; 912 return CEED_ERROR_SUCCESS; 913 } 914 915 /** 916 @brief Get reference coordinates of quadrature points (in dim dimensions) 917 of a CeedBasis 918 919 @param basis CeedBasis 920 @param[out] q_ref Variable to store reference coordinates of quadrature points 921 922 @return An error code: 0 - success, otherwise - failure 923 924 @ref Advanced 925 **/ 926 int CeedBasisGetQRef(CeedBasis basis, const CeedScalar **q_ref) { 927 *q_ref = basis->q_ref_1d; 928 return CEED_ERROR_SUCCESS; 929 } 930 931 /** 932 @brief Get quadrature weights of quadrature points (in dim dimensions) 933 of a CeedBasis 934 935 @param basis CeedBasis 936 @param[out] q_weight Variable to store quadrature weights 937 938 @return An error code: 0 - success, otherwise - failure 939 940 @ref Advanced 941 **/ 942 int CeedBasisGetQWeights(CeedBasis basis, const CeedScalar **q_weight) { 943 *q_weight = basis->q_weight_1d; 944 return CEED_ERROR_SUCCESS; 945 } 946 947 /** 948 @brief Get interpolation matrix of a CeedBasis 949 950 @param basis CeedBasis 951 @param[out] interp Variable to store interpolation matrix 952 953 @return An error code: 0 - success, otherwise - failure 954 955 @ref Advanced 956 **/ 957 int CeedBasisGetInterp(CeedBasis basis, const CeedScalar **interp) { 958 if (!basis->interp && basis->tensor_basis) { 959 // Allocate 960 int ierr; 961 ierr = CeedMalloc(basis->Q*basis->P, &basis->interp); CeedChk(ierr); 962 963 // Initialize 964 for (CeedInt i=0; i<basis->Q*basis->P; i++) 965 basis->interp[i] = 1.0; 966 967 // Calculate 968 for (CeedInt d=0; d<basis->dim; d++) 969 for (CeedInt qpt=0; qpt<basis->Q; qpt++) 970 for (CeedInt node=0; node<basis->P; node++) { 971 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 972 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 973 basis->interp[qpt*(basis->P)+node] *= basis->interp_1d[q*basis->P_1d+p]; 974 } 975 } 976 *interp = basis->interp; 977 return CEED_ERROR_SUCCESS; 978 } 979 980 /** 981 @brief Get 1D interpolation matrix of a tensor product CeedBasis 982 983 @param basis CeedBasis 984 @param[out] interp_1d Variable to store interpolation matrix 985 986 @return An error code: 0 - success, otherwise - failure 987 988 @ref Backend 989 **/ 990 int CeedBasisGetInterp1D(CeedBasis basis, const CeedScalar **interp_1d) { 991 if (!basis->tensor_basis) 992 // LCOV_EXCL_START 993 return CeedError(basis->ceed, CEED_ERROR_MINOR, 994 "CeedBasis is not a tensor product basis."); 995 // LCOV_EXCL_STOP 996 997 *interp_1d = basis->interp_1d; 998 return CEED_ERROR_SUCCESS; 999 } 1000 1001 /** 1002 @brief Get gradient matrix of a CeedBasis 1003 1004 @param basis CeedBasis 1005 @param[out] grad Variable to store gradient matrix 1006 1007 @return An error code: 0 - success, otherwise - failure 1008 1009 @ref Advanced 1010 **/ 1011 int CeedBasisGetGrad(CeedBasis basis, const CeedScalar **grad) { 1012 if (!basis->grad && basis->tensor_basis) { 1013 // Allocate 1014 int ierr; 1015 ierr = CeedMalloc(basis->dim*basis->Q*basis->P, &basis->grad); 1016 CeedChk(ierr); 1017 1018 // Initialize 1019 for (CeedInt i=0; i<basis->dim*basis->Q*basis->P; i++) 1020 basis->grad[i] = 1.0; 1021 1022 // Calculate 1023 for (CeedInt d=0; d<basis->dim; d++) 1024 for (CeedInt i=0; i<basis->dim; i++) 1025 for (CeedInt qpt=0; qpt<basis->Q; qpt++) 1026 for (CeedInt node=0; node<basis->P; node++) { 1027 CeedInt p = (node / CeedIntPow(basis->P_1d, d)) % basis->P_1d; 1028 CeedInt q = (qpt / CeedIntPow(basis->Q_1d, d)) % basis->Q_1d; 1029 if (i == d) 1030 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *= 1031 basis->grad_1d[q*basis->P_1d+p]; 1032 else 1033 basis->grad[(i*basis->Q+qpt)*(basis->P)+node] *= 1034 basis->interp_1d[q*basis->P_1d+p]; 1035 } 1036 } 1037 *grad = basis->grad; 1038 return CEED_ERROR_SUCCESS; 1039 } 1040 1041 /** 1042 @brief Get 1D gradient matrix of a tensor product CeedBasis 1043 1044 @param basis CeedBasis 1045 @param[out] grad_1d Variable to store gradient matrix 1046 1047 @return An error code: 0 - success, otherwise - failure 1048 1049 @ref Advanced 1050 **/ 1051 int CeedBasisGetGrad1D(CeedBasis basis, const CeedScalar **grad_1d) { 1052 if (!basis->tensor_basis) 1053 // LCOV_EXCL_START 1054 return CeedError(basis->ceed, CEED_ERROR_MINOR, 1055 "CeedBasis is not a tensor product basis."); 1056 // LCOV_EXCL_STOP 1057 1058 *grad_1d = basis->grad_1d; 1059 return CEED_ERROR_SUCCESS; 1060 } 1061 1062 /** 1063 @brief Destroy a CeedBasis 1064 1065 @param basis CeedBasis to destroy 1066 1067 @return An error code: 0 - success, otherwise - failure 1068 1069 @ref User 1070 **/ 1071 int CeedBasisDestroy(CeedBasis *basis) { 1072 int ierr; 1073 1074 if (!*basis || --(*basis)->ref_count > 0) return CEED_ERROR_SUCCESS; 1075 if ((*basis)->Destroy) { 1076 ierr = (*basis)->Destroy(*basis); CeedChk(ierr); 1077 } 1078 if ((*basis)->contract) { 1079 ierr = CeedTensorContractDestroy(&(*basis)->contract); CeedChk(ierr); 1080 } 1081 ierr = CeedFree(&(*basis)->interp); CeedChk(ierr); 1082 ierr = CeedFree(&(*basis)->interp_1d); CeedChk(ierr); 1083 ierr = CeedFree(&(*basis)->grad); CeedChk(ierr); 1084 ierr = CeedFree(&(*basis)->grad_1d); CeedChk(ierr); 1085 ierr = CeedFree(&(*basis)->q_ref_1d); CeedChk(ierr); 1086 ierr = CeedFree(&(*basis)->q_weight_1d); CeedChk(ierr); 1087 ierr = CeedDestroy(&(*basis)->ceed); CeedChk(ierr); 1088 ierr = CeedFree(basis); CeedChk(ierr); 1089 return CEED_ERROR_SUCCESS; 1090 } 1091 1092 /** 1093 @brief Construct a Gauss-Legendre quadrature 1094 1095 @param Q Number of quadrature points (integrates polynomials of 1096 degree 2*Q-1 exactly) 1097 @param[out] q_ref_1d Array of length Q to hold the abscissa on [-1, 1] 1098 @param[out] q_weight_1d Array of length Q to hold the weights 1099 1100 @return An error code: 0 - success, otherwise - failure 1101 1102 @ref Utility 1103 **/ 1104 int CeedGaussQuadrature(CeedInt Q, CeedScalar *q_ref_1d, 1105 CeedScalar *q_weight_1d) { 1106 // Allocate 1107 CeedScalar P0, P1, P2, dP2, xi, wi, PI = 4.0*atan(1.0); 1108 // Build q_ref_1d, q_weight_1d 1109 for (int i = 0; i <= Q/2; i++) { 1110 // Guess 1111 xi = cos(PI*(CeedScalar)(2*i+1)/((CeedScalar)(2*Q))); 1112 // Pn(xi) 1113 P0 = 1.0; 1114 P1 = xi; 1115 P2 = 0.0; 1116 for (int j = 2; j <= Q; j++) { 1117 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1118 P0 = P1; 1119 P1 = P2; 1120 } 1121 // First Newton Step 1122 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1123 xi = xi-P2/dP2; 1124 // Newton to convergence 1125 for (int k=0; k<100 && fabs(P2)>10*CEED_EPSILON; k++) { 1126 P0 = 1.0; 1127 P1 = xi; 1128 for (int j = 2; j <= Q; j++) { 1129 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1130 P0 = P1; 1131 P1 = P2; 1132 } 1133 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1134 xi = xi-P2/dP2; 1135 } 1136 // Save xi, wi 1137 wi = 2.0/((1.0-xi*xi)*dP2*dP2); 1138 q_weight_1d[i] = wi; 1139 q_weight_1d[Q-1-i] = wi; 1140 q_ref_1d[i] = -xi; 1141 q_ref_1d[Q-1-i]= xi; 1142 } 1143 return CEED_ERROR_SUCCESS; 1144 } 1145 1146 /** 1147 @brief Construct a Gauss-Legendre-Lobatto quadrature 1148 1149 @param Q Number of quadrature points (integrates polynomials of 1150 degree 2*Q-3 exactly) 1151 @param[out] q_ref_1d Array of length Q to hold the abscissa on [-1, 1] 1152 @param[out] q_weight_1d Array of length Q to hold the weights 1153 1154 @return An error code: 0 - success, otherwise - failure 1155 1156 @ref Utility 1157 **/ 1158 int CeedLobattoQuadrature(CeedInt Q, CeedScalar *q_ref_1d, 1159 CeedScalar *q_weight_1d) { 1160 // Allocate 1161 CeedScalar P0, P1, P2, dP2, d2P2, xi, wi, PI = 4.0*atan(1.0); 1162 // Build q_ref_1d, q_weight_1d 1163 // Set endpoints 1164 if (Q < 2) 1165 // LCOV_EXCL_START 1166 return CeedError(NULL, CEED_ERROR_DIMENSION, 1167 "Cannot create Lobatto quadrature with Q=%d < 2 points", Q); 1168 // LCOV_EXCL_STOP 1169 wi = 2.0/((CeedScalar)(Q*(Q-1))); 1170 if (q_weight_1d) { 1171 q_weight_1d[0] = wi; 1172 q_weight_1d[Q-1] = wi; 1173 } 1174 q_ref_1d[0] = -1.0; 1175 q_ref_1d[Q-1] = 1.0; 1176 // Interior 1177 for (int i = 1; i <= (Q-1)/2; i++) { 1178 // Guess 1179 xi = cos(PI*(CeedScalar)(i)/(CeedScalar)(Q-1)); 1180 // Pn(xi) 1181 P0 = 1.0; 1182 P1 = xi; 1183 P2 = 0.0; 1184 for (int j = 2; j < Q; j++) { 1185 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1186 P0 = P1; 1187 P1 = P2; 1188 } 1189 // First Newton step 1190 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1191 d2P2 = (2*xi*dP2 - (CeedScalar)(Q*(Q-1))*P2)/(1.0-xi*xi); 1192 xi = xi-dP2/d2P2; 1193 // Newton to convergence 1194 for (int k=0; k<100 && fabs(dP2)>10*CEED_EPSILON; k++) { 1195 P0 = 1.0; 1196 P1 = xi; 1197 for (int j = 2; j < Q; j++) { 1198 P2 = (((CeedScalar)(2*j-1))*xi*P1-((CeedScalar)(j-1))*P0)/((CeedScalar)(j)); 1199 P0 = P1; 1200 P1 = P2; 1201 } 1202 dP2 = (xi*P2 - P0)*(CeedScalar)Q/(xi*xi-1.0); 1203 d2P2 = (2*xi*dP2 - (CeedScalar)(Q*(Q-1))*P2)/(1.0-xi*xi); 1204 xi = xi-dP2/d2P2; 1205 } 1206 // Save xi, wi 1207 wi = 2.0/(((CeedScalar)(Q*(Q-1)))*P2*P2); 1208 if (q_weight_1d) { 1209 q_weight_1d[i] = wi; 1210 q_weight_1d[Q-1-i] = wi; 1211 } 1212 q_ref_1d[i] = -xi; 1213 q_ref_1d[Q-1-i]= xi; 1214 } 1215 return CEED_ERROR_SUCCESS; 1216 } 1217 1218 /** 1219 @brief Return QR Factorization of a matrix 1220 1221 @param ceed A Ceed context for error handling 1222 @param[in,out] mat Row-major matrix to be factorized in place 1223 @param[in,out] tau Vector of length m of scaling factors 1224 @param m Number of rows 1225 @param n Number of columns 1226 1227 @return An error code: 0 - success, otherwise - failure 1228 1229 @ref Utility 1230 **/ 1231 int CeedQRFactorization(Ceed ceed, CeedScalar *mat, CeedScalar *tau, 1232 CeedInt m, CeedInt n) { 1233 CeedScalar v[m]; 1234 1235 // Check m >= n 1236 if (n > m) 1237 // LCOV_EXCL_START 1238 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1239 "Cannot compute QR factorization with n > m"); 1240 // LCOV_EXCL_STOP 1241 1242 for (CeedInt i=0; i<n; i++) { 1243 if (i >= m-1) { // last row of matrix, no reflection needed 1244 tau[i] = 0.; 1245 break; 1246 } 1247 // Calculate Householder vector, magnitude 1248 CeedScalar sigma = 0.0; 1249 v[i] = mat[i+n*i]; 1250 for (CeedInt j=i+1; j<m; j++) { 1251 v[j] = mat[i+n*j]; 1252 sigma += v[j] * v[j]; 1253 } 1254 CeedScalar norm = sqrt(v[i]*v[i] + sigma); // norm of v[i:m] 1255 CeedScalar R_ii = -copysign(norm, v[i]); 1256 v[i] -= R_ii; 1257 // norm of v[i:m] after modification above and scaling below 1258 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1259 // tau = 2 / (norm*norm) 1260 tau[i] = 2 * v[i]*v[i] / (v[i]*v[i] + sigma); 1261 for (CeedInt j=i+1; j<m; j++) 1262 v[j] /= v[i]; 1263 1264 // Apply Householder reflector to lower right panel 1265 CeedHouseholderReflect(&mat[i*n+i+1], &v[i], tau[i], m-i, n-i-1, n, 1); 1266 // Save v 1267 mat[i+n*i] = R_ii; 1268 for (CeedInt j=i+1; j<m; j++) 1269 mat[i+n*j] = v[j]; 1270 } 1271 return CEED_ERROR_SUCCESS; 1272 } 1273 1274 /** 1275 @brief Return symmetric Schur decomposition of the symmetric matrix mat via 1276 symmetric QR factorization 1277 1278 @param ceed A Ceed context for error handling 1279 @param[in,out] mat Row-major matrix to be factorized in place 1280 @param[out] lambda Vector of length n of eigenvalues 1281 @param n Number of rows/columns 1282 1283 @return An error code: 0 - success, otherwise - failure 1284 1285 @ref Utility 1286 **/ 1287 CeedPragmaOptimizeOff 1288 int CeedSymmetricSchurDecomposition(Ceed ceed, CeedScalar *mat, 1289 CeedScalar *lambda, CeedInt n) { 1290 // Check bounds for clang-tidy 1291 if (n<2) 1292 // LCOV_EXCL_START 1293 return CeedError(ceed, CEED_ERROR_UNSUPPORTED, 1294 "Cannot compute symmetric Schur decomposition of scalars"); 1295 // LCOV_EXCL_STOP 1296 1297 CeedScalar v[n-1], tau[n-1], mat_T[n*n]; 1298 1299 // Copy mat to mat_T and set mat to I 1300 memcpy(mat_T, mat, n*n*sizeof(mat[0])); 1301 for (CeedInt i=0; i<n; i++) 1302 for (CeedInt j=0; j<n; j++) 1303 mat[j+n*i] = (i==j) ? 1 : 0; 1304 1305 // Reduce to tridiagonal 1306 for (CeedInt i=0; i<n-1; i++) { 1307 // Calculate Householder vector, magnitude 1308 CeedScalar sigma = 0.0; 1309 v[i] = mat_T[i+n*(i+1)]; 1310 for (CeedInt j=i+1; j<n-1; j++) { 1311 v[j] = mat_T[i+n*(j+1)]; 1312 sigma += v[j] * v[j]; 1313 } 1314 CeedScalar norm = sqrt(v[i]*v[i] + sigma); // norm of v[i:n-1] 1315 CeedScalar R_ii = -copysign(norm, v[i]); 1316 v[i] -= R_ii; 1317 // norm of v[i:m] after modification above and scaling below 1318 // norm = sqrt(v[i]*v[i] + sigma) / v[i]; 1319 // tau = 2 / (norm*norm) 1320 tau[i] = i == n - 2 ? 2 : 2 * v[i]*v[i] / (v[i]*v[i] + sigma); 1321 for (CeedInt j=i+1; j<n-1; j++) 1322 v[j] /= v[i]; 1323 1324 // Update sub and super diagonal 1325 for (CeedInt j=i+2; j<n; j++) { 1326 mat_T[i+n*j] = 0; mat_T[j+n*i] = 0; 1327 } 1328 // Apply symmetric Householder reflector to lower right panel 1329 CeedHouseholderReflect(&mat_T[(i+1)+n*(i+1)], &v[i], tau[i], 1330 n-(i+1), n-(i+1), n, 1); 1331 CeedHouseholderReflect(&mat_T[(i+1)+n*(i+1)], &v[i], tau[i], 1332 n-(i+1), n-(i+1), 1, n); 1333 1334 // Save v 1335 mat_T[i+n*(i+1)] = R_ii; 1336 mat_T[(i+1)+n*i] = R_ii; 1337 for (CeedInt j=i+1; j<n-1; j++) { 1338 mat_T[i+n*(j+1)] = v[j]; 1339 } 1340 } 1341 // Backwards accumulation of Q 1342 for (CeedInt i=n-2; i>=0; i--) { 1343 if (tau[i] > 0.0) { 1344 v[i] = 1; 1345 for (CeedInt j=i+1; j<n-1; j++) { 1346 v[j] = mat_T[i+n*(j+1)]; 1347 mat_T[i+n*(j+1)] = 0; 1348 } 1349 CeedHouseholderReflect(&mat[(i+1)+n*(i+1)], &v[i], tau[i], 1350 n-(i+1), n-(i+1), n, 1); 1351 } 1352 } 1353 1354 // Reduce sub and super diagonal 1355 CeedInt p = 0, q = 0, itr = 0, max_itr = n*n*n*n; 1356 CeedScalar tol = CEED_EPSILON; 1357 1358 while (itr < max_itr) { 1359 // Update p, q, size of reduced portions of diagonal 1360 p = 0; q = 0; 1361 for (CeedInt i=n-2; i>=0; i--) { 1362 if (fabs(mat_T[i+n*(i+1)]) < tol) 1363 q += 1; 1364 else 1365 break; 1366 } 1367 for (CeedInt i=0; i<n-q-1; i++) { 1368 if (fabs(mat_T[i+n*(i+1)]) < tol) 1369 p += 1; 1370 else 1371 break; 1372 } 1373 if (q == n-1) break; // Finished reducing 1374 1375 // Reduce tridiagonal portion 1376 CeedScalar t_nn = mat_T[(n-1-q)+n*(n-1-q)], 1377 t_nnm1 = mat_T[(n-2-q)+n*(n-1-q)]; 1378 CeedScalar d = (mat_T[(n-2-q)+n*(n-2-q)] - t_nn)/2; 1379 CeedScalar mu = t_nn - t_nnm1*t_nnm1 / 1380 (d + copysign(sqrt(d*d + t_nnm1*t_nnm1), d)); 1381 CeedScalar x = mat_T[p+n*p] - mu; 1382 CeedScalar z = mat_T[p+n*(p+1)]; 1383 for (CeedInt k=p; k<n-q-1; k++) { 1384 // Compute Givens rotation 1385 CeedScalar c = 1, s = 0; 1386 if (fabs(z) > tol) { 1387 if (fabs(z) > fabs(x)) { 1388 CeedScalar tau = -x/z; 1389 s = 1/sqrt(1+tau*tau), c = s*tau; 1390 } else { 1391 CeedScalar tau = -z/x; 1392 c = 1/sqrt(1+tau*tau), s = c*tau; 1393 } 1394 } 1395 1396 // Apply Givens rotation to T 1397 CeedGivensRotation(mat_T, c, s, CEED_NOTRANSPOSE, k, k+1, n, n); 1398 CeedGivensRotation(mat_T, c, s, CEED_TRANSPOSE, k, k+1, n, n); 1399 1400 // Apply Givens rotation to Q 1401 CeedGivensRotation(mat, c, s, CEED_NOTRANSPOSE, k, k+1, n, n); 1402 1403 // Update x, z 1404 if (k < n-q-2) { 1405 x = mat_T[k+n*(k+1)]; 1406 z = mat_T[k+n*(k+2)]; 1407 } 1408 } 1409 itr++; 1410 } 1411 1412 // Save eigenvalues 1413 for (CeedInt i=0; i<n; i++) 1414 lambda[i] = mat_T[i+n*i]; 1415 1416 // Check convergence 1417 if (itr == max_itr && q < n-1) 1418 // LCOV_EXCL_START 1419 return CeedError(ceed, CEED_ERROR_MINOR, 1420 "Symmetric QR failed to converge"); 1421 // LCOV_EXCL_STOP 1422 return CEED_ERROR_SUCCESS; 1423 } 1424 CeedPragmaOptimizeOn 1425 1426 /** 1427 @brief Return Simultaneous Diagonalization of two matrices. This solves the 1428 generalized eigenvalue problem A x = lambda B x, where A and B 1429 are symmetric and B is positive definite. We generate the matrix X 1430 and vector Lambda such that X^T A X = Lambda and X^T B X = I. This 1431 is equivalent to the LAPACK routine 'sygv' with TYPE = 1. 1432 1433 @param ceed A Ceed context for error handling 1434 @param[in] mat_A Row-major matrix to be factorized with eigenvalues 1435 @param[in] mat_B Row-major matrix to be factorized to identity 1436 @param[out] mat_X Row-major orthogonal matrix 1437 @param[out] lambda Vector of length n of generalized eigenvalues 1438 @param n Number of rows/columns 1439 1440 @return An error code: 0 - success, otherwise - failure 1441 1442 @ref Utility 1443 **/ 1444 CeedPragmaOptimizeOff 1445 int CeedSimultaneousDiagonalization(Ceed ceed, CeedScalar *mat_A, 1446 CeedScalar *mat_B, CeedScalar *mat_X, 1447 CeedScalar *lambda, CeedInt n) { 1448 int ierr; 1449 CeedScalar *mat_C, *mat_G, *vec_D; 1450 ierr = CeedCalloc(n*n, &mat_C); CeedChk(ierr); 1451 ierr = CeedCalloc(n*n, &mat_G); CeedChk(ierr); 1452 ierr = CeedCalloc(n, &vec_D); CeedChk(ierr); 1453 1454 // Compute B = G D G^T 1455 memcpy(mat_G, mat_B, n*n*sizeof(mat_B[0])); 1456 ierr = CeedSymmetricSchurDecomposition(ceed, mat_G, vec_D, n); CeedChk(ierr); 1457 1458 // Sort eigenvalues 1459 for (CeedInt i=n-1; i>=0; i--) 1460 for (CeedInt j=0; j<i; j++) { 1461 if (fabs(vec_D[j]) > fabs(vec_D[j+1])) { 1462 CeedScalar temp; 1463 temp = vec_D[j]; vec_D[j] = vec_D[j+1]; vec_D[j+1] = temp; 1464 for (CeedInt k=0; k<n; k++) { 1465 temp = mat_G[k*n+j]; mat_G[k*n+j] = mat_G[k*n+j+1]; mat_G[k*n+j+1] = temp; 1466 } 1467 } 1468 } 1469 1470 // Compute C = (G D^1/2)^-1 A (G D^1/2)^-T 1471 // = D^-1/2 G^T A G D^-1/2 1472 // -- D = D^-1/2 1473 for (CeedInt i=0; i<n; i++) 1474 vec_D[i] = 1./sqrt(vec_D[i]); 1475 // -- G = G D^-1/2 1476 // -- C = D^-1/2 G^T 1477 for (CeedInt i=0; i<n; i++) 1478 for (CeedInt j=0; j<n; j++) { 1479 mat_G[i*n+j] *= vec_D[j]; 1480 mat_C[j*n+i] = mat_G[i*n+j]; 1481 } 1482 // -- X = (D^-1/2 G^T) A 1483 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_C, 1484 (const CeedScalar *)mat_A, mat_X, n, n, n); 1485 CeedChk(ierr); 1486 // -- C = (D^-1/2 G^T A) (G D^-1/2) 1487 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_X, 1488 (const CeedScalar *)mat_G, mat_C, n, n, n); 1489 CeedChk(ierr); 1490 1491 // Compute Q^T C Q = lambda 1492 ierr = CeedSymmetricSchurDecomposition(ceed, mat_C, lambda, n); CeedChk(ierr); 1493 1494 // Sort eigenvalues 1495 for (CeedInt i=n-1; i>=0; i--) 1496 for (CeedInt j=0; j<i; j++) { 1497 if (fabs(lambda[j]) > fabs(lambda[j+1])) { 1498 CeedScalar temp; 1499 temp = lambda[j]; lambda[j] = lambda[j+1]; lambda[j+1] = temp; 1500 for (CeedInt k=0; k<n; k++) { 1501 temp = mat_C[k*n+j]; mat_C[k*n+j] = mat_C[k*n+j+1]; mat_C[k*n+j+1] = temp; 1502 } 1503 } 1504 } 1505 1506 // Set X = (G D^1/2)^-T Q 1507 // = G D^-1/2 Q 1508 ierr = CeedMatrixMultiply(ceed, (const CeedScalar *)mat_G, 1509 (const CeedScalar *)mat_C, mat_X, n, n, n); 1510 CeedChk(ierr); 1511 1512 // Cleanup 1513 ierr = CeedFree(&mat_C); CeedChk(ierr); 1514 ierr = CeedFree(&mat_G); CeedChk(ierr); 1515 ierr = CeedFree(&vec_D); CeedChk(ierr); 1516 return CEED_ERROR_SUCCESS; 1517 } 1518 CeedPragmaOptimizeOn 1519 1520 /// @} 1521