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