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 CeedSimultaneousDiagonalization(ceed, K, M, X, lambda, P); 33 34 // Check X^T M X = I 35 CeedScalar work[P*P]; 36 for (int i=0; i<P; i++) 37 for (int j=0; j<P; j++) { 38 CeedScalar sum = 0; 39 for (int k=0; k<P; k++) 40 sum += M[P*i+k]*X[P*k+j]; 41 work[P*i+j] = sum; 42 } 43 for (int i=0; i<P; i++) 44 for (int j=0; j<P; j++) { 45 CeedScalar sum = 0; 46 for (int k=0; k<P; k++) 47 sum += X[P*k+i]*work[P*k+j]; 48 M[P*i+j] = sum; 49 } 50 for (int i=0; i<P; i++) 51 for (int j=0; j<P; j++) 52 if (fabs(M[P*i+j] - (i == j ? 1.0 : 0.0)) > 100.*CEED_EPSILON) 53 // LCOV_EXCL_START 54 printf("Error in diagonalization of M [%" CeedInt_FMT 55 ", %" CeedInt_FMT "]: %f != %f\n", 56 i, j, M[P*i+j], (i == j ? 1.0 : 0.0)); 57 // LCOV_EXCL_STOP 58 59 // Check X^T K X = Lamda 60 for (int i=0; i<P; i++) 61 for (int j=0; j<P; j++) { 62 CeedScalar sum = 0; 63 for (int k=0; k<P; k++) 64 sum += K[P*i+k]*X[P*k+j]; 65 work[P*i+j] = sum; 66 } 67 for (int i=0; i<P; i++) 68 for (int j=0; j<P; j++) { 69 CeedScalar sum = 0; 70 for (int k=0; k<P; k++) 71 sum += X[P*k+i]*work[P*k+j]; 72 K[P*i+j] = sum; 73 } 74 for (int i=0; i<P; i++) 75 for (int j=0; j<P; j++) 76 if (fabs(K[P*i+j] - (i == j ? lambda[i] : 0.0)) > 100.*CEED_EPSILON) 77 // LCOV_EXCL_START 78 printf("Error in diagonalization of K [%" CeedInt_FMT 79 ", %" CeedInt_FMT "]: %f != %f\n", 80 i, j, K[P*i+j], (i == j ? lambda[i] : 0.0)); 81 // LCOV_EXCL_STOP 82 83 CeedBasisDestroy(&basis); 84 CeedDestroy(&ceed); 85 return 0; 86 } 87