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