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