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 const CeedScalar *interpolation, *gradient, *quadrature_weights; 13 14 CeedInit(argv[1], &ceed); 15 16 // Create mass, stiffness matrix 17 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis); 18 CeedBasisGetInterp(basis, &interpolation); 19 CeedBasisGetGrad(basis, &gradient); 20 CeedBasisGetQWeights(basis, &quadrature_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 += interpolation[p * k + i] * quadrature_weights[k] * interpolation[p * k + j]; 26 sum_k += gradient[p * k + i] * quadrature_weights[k] * gradient[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_array[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_array[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_array[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 = Lambda 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_array[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_array[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