14411cf47Sjeremylt /// @file 252bfb9bbSJeremy L Thompson /// Test Symmetric Schur Decomposition 352bfb9bbSJeremy L Thompson /// \test Test Symmetric Schur Decomposition 457c64913Sjeremylt #include <ceed.h> 5673160d7Sjeremylt #include <ceed/backend.h> 6673160d7Sjeremylt #include <math.h> 757c64913Sjeremylt 857c64913Sjeremylt int main(int argc, char **argv) { 957c64913Sjeremylt Ceed ceed; 10*4fee36f0SJeremy L Thompson CeedInt p = 4; 11673160d7Sjeremylt CeedScalar M[16], Q[16], lambda[4], Q_lambda_Qt[16]; 12673160d7Sjeremylt CeedBasis basis; 13*4fee36f0SJeremy L Thompson const CeedScalar *interpolation, *quadrature_weights; 1457c64913Sjeremylt 1557c64913Sjeremylt CeedInit(argv[1], &ceed); 16288c0443SJeremy L Thompson 17673160d7Sjeremylt // Create mass matrix 18*4fee36f0SJeremy L Thompson CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, p, CEED_GAUSS, &basis); 19*4fee36f0SJeremy L Thompson CeedBasisGetInterp(basis, &interpolation); 20*4fee36f0SJeremy L Thompson CeedBasisGetQWeights(basis, &quadrature_weights); 21*4fee36f0SJeremy L Thompson for (int i = 0; i < p; i++) { 22*4fee36f0SJeremy L Thompson for (int j = 0; j < p; j++) { 23673160d7Sjeremylt CeedScalar sum = 0; 24*4fee36f0SJeremy L Thompson for (int k = 0; k < p; k++) sum += interpolation[p * k + i] * quadrature_weights[k] * interpolation[p * k + j]; 25*4fee36f0SJeremy L Thompson M[p * i + j] = sum; 26*4fee36f0SJeremy L Thompson Q[p * i + j] = sum; 2757c64913Sjeremylt } 282b730f8bSJeremy L Thompson } 29673160d7Sjeremylt 30*4fee36f0SJeremy L Thompson CeedSymmetricSchurDecomposition(ceed, Q, lambda, p); 31673160d7Sjeremylt 32673160d7Sjeremylt // Check diagonalization of M 33*4fee36f0SJeremy L Thompson for (int i = 0; i < p; i++) { 34*4fee36f0SJeremy L Thompson for (int j = 0; j < p; j++) { 35673160d7Sjeremylt CeedScalar sum = 0; 36*4fee36f0SJeremy L Thompson for (int k = 0; k < p; k++) sum += Q[p * i + k] * lambda[k] * Q[p * j + k]; 37*4fee36f0SJeremy L Thompson Q_lambda_Qt[p * i + j] = sum; 3852bfb9bbSJeremy L Thompson } 392b730f8bSJeremy L Thompson } 40*4fee36f0SJeremy L Thompson for (int i = 0; i < p; i++) { 41*4fee36f0SJeremy L Thompson for (int j = 0; j < p; j++) { 42*4fee36f0SJeremy L Thompson if (fabs(M[p * i + j] - Q_lambda_Qt[p * i + j]) > 100. * CEED_EPSILON) { 43673160d7Sjeremylt // LCOV_EXCL_START 44*4fee36f0SJeremy L Thompson printf("Error in diagonalization [%" CeedInt_FMT ", %" CeedInt_FMT "]: %f != %f\n", i, j, M[p * i + j], Q_lambda_Qt[p * i + j]); 45673160d7Sjeremylt // LCOV_EXCL_STOP 462b730f8bSJeremy L Thompson } 472b730f8bSJeremy L Thompson } 482b730f8bSJeremy L Thompson } 49673160d7Sjeremylt 50673160d7Sjeremylt CeedBasisDestroy(&basis); 5157c64913Sjeremylt CeedDestroy(&ceed); 5257c64913Sjeremylt return 0; 5357c64913Sjeremylt } 54