14411cf47Sjeremylt /// @file 252bfb9bbSJeremy L Thompson /// Test Symmetric Schur Decomposition 352bfb9bbSJeremy L Thompson /// \test Test Symmetric Schur Decomposition 457c64913Sjeremylt #include <ceed.h> 5*673160d7Sjeremylt #include <ceed/backend.h> 6*673160d7Sjeremylt #include <math.h> 757c64913Sjeremylt 857c64913Sjeremylt int main(int argc, char **argv) { 957c64913Sjeremylt Ceed ceed; 10*673160d7Sjeremylt CeedInt P = 4; 11*673160d7Sjeremylt CeedScalar M[16], Q[16], lambda[4], Q_lambda_Qt[16]; 12*673160d7Sjeremylt CeedBasis basis; 1357c64913Sjeremylt 1457c64913Sjeremylt CeedInit(argv[1], &ceed); 15288c0443SJeremy L Thompson 16*673160d7Sjeremylt // Create mass matrix 17*673160d7Sjeremylt CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P, CEED_GAUSS, &basis); 18*673160d7Sjeremylt const CeedScalar *interp, *quad_weights; 19*673160d7Sjeremylt CeedBasisGetInterp(basis, &interp); 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 = 0; 24*673160d7Sjeremylt for (int k=0; k<P; k++) 25*673160d7Sjeremylt sum += interp[P*k+i]*quad_weights[k]*interp[P*k+j]; 26*673160d7Sjeremylt M[P*i+j] = sum; 27*673160d7Sjeremylt Q[P*i+j] = sum; 2857c64913Sjeremylt } 29*673160d7Sjeremylt 30*673160d7Sjeremylt CeedSymmetricSchurDecomposition(ceed, Q, lambda, P); 31*673160d7Sjeremylt 32*673160d7Sjeremylt // Check diagonalization of M 33*673160d7Sjeremylt for (int i=0; i<P; i++) 34*673160d7Sjeremylt for (int j=0; j<P; j++) { 35*673160d7Sjeremylt CeedScalar sum = 0; 36*673160d7Sjeremylt for (int k=0; k<P; k++) 37*673160d7Sjeremylt sum += Q[P*i+k]*lambda[k]*Q[P*j+k]; 38*673160d7Sjeremylt Q_lambda_Qt[P*i+j] = sum; 3952bfb9bbSJeremy L Thompson } 40*673160d7Sjeremylt for (int i=0; i<P; i++) 41*673160d7Sjeremylt for (int j=0; j<P; j++) 42*673160d7Sjeremylt if (fabs(M[P*i+j] - Q_lambda_Qt[P*i+j]) > 1E-14) 43*673160d7Sjeremylt // LCOV_EXCL_START 44*673160d7Sjeremylt printf("Error in diagonalization [%d, %d]: %f != %f\n", 45*673160d7Sjeremylt i, j, M[P*i+j], Q_lambda_Qt[P*i+j]); 46*673160d7Sjeremylt // LCOV_EXCL_STOP 47*673160d7Sjeremylt 48*673160d7Sjeremylt CeedBasisDestroy(&basis); 4957c64913Sjeremylt CeedDestroy(&ceed); 5057c64913Sjeremylt return 0; 5157c64913Sjeremylt } 52