14411cf47Sjeremylt /// @file 252bfb9bbSJeremy L Thompson /// Test Simultaneous Diagonalization 352bfb9bbSJeremy L Thompson /// \test Simultaneous Diagonalization 457c64913Sjeremylt #include <ceed.h> 5*673160d7Sjeremylt #include <math.h> 657c64913Sjeremylt 757c64913Sjeremylt int main(int argc, char **argv) { 857c64913Sjeremylt Ceed ceed; 9*673160d7Sjeremylt CeedInt P = 4, Q = 4; 10*673160d7Sjeremylt CeedScalar M[P*P], K[P*P], X[P*P], lambda[P]; 11*673160d7Sjeremylt CeedBasis basis; 1257c64913Sjeremylt 1357c64913Sjeremylt CeedInit(argv[1], &ceed); 14aedaa0e5Sjeremylt 15*673160d7Sjeremylt // Create mass, stiffness matrix 16*673160d7Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, Q, CEED_GAUSS, &basis); 17*673160d7Sjeremylt const CeedScalar *interp, *grad, *quad_weights; 18*673160d7Sjeremylt CeedBasisGetInterp(basis, &interp); 19*673160d7Sjeremylt CeedBasisGetGrad(basis, &grad); 20*673160d7Sjeremylt CeedBasisGetQWeights(basis, &quad_weights); 21*673160d7Sjeremylt for (int i=0; i<P; i++) 22*673160d7Sjeremylt for (int j=0; j<P; j++) { 23*673160d7Sjeremylt CeedScalar sum_m = 0, sum_k = 0; 24*673160d7Sjeremylt for (int k=0; k<Q; k++) { 25*673160d7Sjeremylt sum_m += interp[P*k+i]*quad_weights[k]*interp[P*k+j]; 26*673160d7Sjeremylt sum_k += grad[P*k+i]*quad_weights[k]*grad[P*k+j]; 27*673160d7Sjeremylt } 28*673160d7Sjeremylt M[P*i+j] = sum_m; 29*673160d7Sjeremylt K[P*i+j] = sum_k; 30fb551037Sjeremylt } 31fb551037Sjeremylt 32*673160d7Sjeremylt CeedSimultaneousDiagonalization(ceed, K, M, X, lambda, P); 33*673160d7Sjeremylt 34*673160d7Sjeremylt // Check X^T M X = I 35*673160d7Sjeremylt CeedScalar work[P*P]; 36*673160d7Sjeremylt for (int i=0; i<P; i++) 37*673160d7Sjeremylt for (int j=0; j<P; j++) { 38*673160d7Sjeremylt CeedScalar sum = 0; 39*673160d7Sjeremylt for (int k=0; k<P; k++) 40*673160d7Sjeremylt sum += M[P*i+k]*X[P*k+j]; 41*673160d7Sjeremylt work[P*i+j] = sum; 4252bfb9bbSJeremy L Thompson } 43*673160d7Sjeremylt for (int i=0; i<P; i++) 44*673160d7Sjeremylt for (int j=0; j<P; j++) { 45*673160d7Sjeremylt CeedScalar sum = 0; 46*673160d7Sjeremylt for (int k=0; k<P; k++) 47*673160d7Sjeremylt sum += X[P*k+i]*work[P*k+j]; 48*673160d7Sjeremylt M[P*i+j] = sum; 4952bfb9bbSJeremy L Thompson } 50*673160d7Sjeremylt for (int i=0; i<P; i++) 51*673160d7Sjeremylt for (int j=0; j<P; j++) 52*673160d7Sjeremylt if (fabs(M[P*i+j] - (i == j ? 1.0 : 0.0)) > 1E-14) 53*673160d7Sjeremylt // LCOV_EXCL_START 54*673160d7Sjeremylt printf("Error in diagonalization of M [%d, %d]: %f != %f\n", 55*673160d7Sjeremylt i, j, M[P*i+j], (i == j ? 1.0 : 0.0)); 56*673160d7Sjeremylt // LCOV_EXCL_STOP 57*673160d7Sjeremylt 58*673160d7Sjeremylt // Check X^T K X = Lamda 59*673160d7Sjeremylt for (int i=0; i<P; i++) 60*673160d7Sjeremylt for (int j=0; j<P; j++) { 61*673160d7Sjeremylt CeedScalar sum = 0; 62*673160d7Sjeremylt for (int k=0; k<P; k++) 63*673160d7Sjeremylt sum += K[P*i+k]*X[P*k+j]; 64*673160d7Sjeremylt work[P*i+j] = sum; 6552bfb9bbSJeremy L Thompson } 66*673160d7Sjeremylt for (int i=0; i<P; i++) 67*673160d7Sjeremylt for (int j=0; j<P; j++) { 68*673160d7Sjeremylt CeedScalar sum = 0; 69*673160d7Sjeremylt for (int k=0; k<P; k++) 70*673160d7Sjeremylt sum += X[P*k+i]*work[P*k+j]; 71*673160d7Sjeremylt K[P*i+j] = sum; 72*673160d7Sjeremylt } 73*673160d7Sjeremylt for (int i=0; i<P; i++) 74*673160d7Sjeremylt for (int j=0; j<P; j++) 75*673160d7Sjeremylt if (fabs(K[P*i+j] - (i == j ? lambda[i] : 0.0)) > 1E-14) 76*673160d7Sjeremylt // LCOV_EXCL_START 77*673160d7Sjeremylt printf("Error in diagonalization of K [%d, %d]: %f != %f\n", 78*673160d7Sjeremylt i, j, K[P*i+j], (i == j ? lambda[i] : 0.0)); 79*673160d7Sjeremylt // LCOV_EXCL_STOP 80*673160d7Sjeremylt 81*673160d7Sjeremylt CeedBasisDestroy(&basis); 8257c64913Sjeremylt CeedDestroy(&ceed); 8357c64913Sjeremylt return 0; 8457c64913Sjeremylt } 85