1 /// @file 2 /// Test Symmetric Schur Decomposition 3 /// \test Test Symmetric Schur Decomposition 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; 14 CeedScalar M[16], Q[16], lambda[4], Q_lambda_Qt[16]; 15 CeedBasis basis; 16 const CeedScalar *interpolation, *quadrature_weights; 17 18 CeedInit(argv[1], &ceed); 19 20 // Create mass matrix 21 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, p, CEED_GAUSS, &basis); 22 CeedBasisGetInterp(basis, &interpolation); 23 CeedBasisGetQWeights(basis, &quadrature_weights); 24 for (int i = 0; i < p; i++) { 25 for (int j = 0; j < p; j++) { 26 CeedScalar sum = 0; 27 for (int k = 0; k < p; k++) sum += interpolation[p * k + i] * quadrature_weights[k] * interpolation[p * k + j]; 28 M[p * i + j] = sum; 29 Q[p * i + j] = sum; 30 } 31 } 32 33 CeedSymmetricSchurDecomposition(ceed, Q, lambda, p); 34 35 // Check diagonalization of M 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++) sum += Q[p * i + k] * lambda[k] * Q[p * j + k]; 40 Q_lambda_Qt[p * i + j] = sum; 41 } 42 } 43 for (int i = 0; i < p; i++) { 44 for (int j = 0; j < p; j++) { 45 if (fabs(M[p * i + j] - Q_lambda_Qt[p * i + j]) > 100. * CEED_EPSILON) { 46 // LCOV_EXCL_START 47 printf("Error in diagonalization [%" CeedInt_FMT ", %" CeedInt_FMT "]: %f != %f\n", i, j, M[p * i + j], Q_lambda_Qt[p * i + j]); 48 // LCOV_EXCL_STOP 49 } 50 } 51 } 52 53 CeedBasisDestroy(&basis); 54 CeedDestroy(&ceed); 55 return 0; 56 } 57