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