1 /// @file 2 /// Test Simultaneous Diagonalization 3 /// \test Simultaneous Diagonalization 4 #include <ceed.h> 5 #include <math.h> 6 7 int main(int argc, char **argv) { 8 Ceed ceed; 9 CeedInt P = 4, Q = 4; 10 CeedScalar M[P * P], K[P * P], X[P * P], lambda[P]; 11 CeedBasis basis; 12 13 CeedInit(argv[1], &ceed); 14 15 // Create mass, stiffness matrix 16 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, Q, CEED_GAUSS, &basis); 17 const CeedScalar *interp, *grad, *quad_weights; 18 CeedBasisGetInterp(basis, &interp); 19 CeedBasisGetGrad(basis, &grad); 20 CeedBasisGetQWeights(basis, &quad_weights); 21 for (int i = 0; i < P; i++) { 22 for (int j = 0; j < P; j++) { 23 CeedScalar sum_m = 0, sum_k = 0; 24 for (int k = 0; k < Q; k++) { 25 sum_m += interp[P * k + i] * quad_weights[k] * interp[P * k + j]; 26 sum_k += grad[P * k + i] * quad_weights[k] * grad[P * k + j]; 27 } 28 M[P * i + j] = sum_m; 29 K[P * i + j] = sum_k; 30 } 31 } 32 33 CeedSimultaneousDiagonalization(ceed, K, M, X, lambda, P); 34 35 // Check X^T M X = I 36 CeedScalar work[P * P]; 37 for (int i = 0; i < P; i++) { 38 for (int j = 0; j < P; j++) { 39 CeedScalar sum = 0; 40 for (int k = 0; k < P; k++) sum += M[P * i + k] * X[P * k + j]; 41 work[P * i + j] = sum; 42 } 43 } 44 for (int i = 0; i < P; i++) { 45 for (int j = 0; j < P; j++) { 46 CeedScalar sum = 0; 47 for (int k = 0; k < P; k++) sum += X[P * k + i] * work[P * k + j]; 48 M[P * i + j] = sum; 49 } 50 } 51 for (int i = 0; i < P; i++) { 52 for (int j = 0; j < P; j++) { 53 if (fabs(M[P * i + j] - (i == j ? 1.0 : 0.0)) > 100. * CEED_EPSILON) { 54 // LCOV_EXCL_START 55 printf("Error in diagonalization of M [%" CeedInt_FMT ", %" CeedInt_FMT "]: %f != %f\n", i, j, M[P * i + j], (i == j ? 1.0 : 0.0)); 56 // LCOV_EXCL_STOP 57 } 58 } 59 } 60 61 // Check X^T K X = Lamda 62 for (int i = 0; i < P; i++) { 63 for (int j = 0; j < P; j++) { 64 CeedScalar sum = 0; 65 for (int k = 0; k < P; k++) sum += K[P * i + k] * X[P * k + j]; 66 work[P * i + j] = sum; 67 } 68 } 69 for (int i = 0; i < P; i++) { 70 for (int j = 0; j < P; j++) { 71 CeedScalar sum = 0; 72 for (int k = 0; k < P; k++) sum += X[P * k + i] * work[P * k + j]; 73 K[P * i + j] = sum; 74 } 75 } 76 for (int i = 0; i < P; i++) { 77 for (int j = 0; j < P; j++) { 78 if (fabs(K[P * i + j] - (i == j ? lambda[i] : 0.0)) > 100. * CEED_EPSILON) { 79 // LCOV_EXCL_START 80 printf("Error in diagonalization of K [%" CeedInt_FMT ", %" CeedInt_FMT "]: %f != %f\n", i, j, K[P * i + j], (i == j ? lambda[i] : 0.0)); 81 // LCOV_EXCL_STOP 82 } 83 } 84 } 85 86 CeedBasisDestroy(&basis); 87 CeedDestroy(&ceed); 88 return 0; 89 } 90