1 /// @file 2 /// Test Symmetric Schur Decomposition 3 /// \test Test Symmetric Schur Decomposition 4 #include <ceed.h> 5 #include <ceed/backend.h> 6 #include <math.h> 7 8 int main(int argc, char **argv) { 9 Ceed ceed; 10 CeedInt P = 4; 11 CeedScalar M[16], Q[16], lambda[4], Q_lambda_Qt[16]; 12 CeedBasis basis; 13 14 CeedInit(argv[1], &ceed); 15 16 // Create mass matrix 17 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, P, P, CEED_GAUSS, &basis); 18 const CeedScalar *interp, *quad_weights; 19 CeedBasisGetInterp(basis, &interp); 20 CeedBasisGetQWeights(basis, &quad_weights); 21 for (int i=0; i<P; i++) 22 for (int j=0; j<P; j++) { 23 CeedScalar sum = 0; 24 for (int k=0; k<P; k++) 25 sum += interp[P*k+i]*quad_weights[k]*interp[P*k+j]; 26 M[P*i+j] = sum; 27 Q[P*i+j] = sum; 28 } 29 30 CeedSymmetricSchurDecomposition(ceed, Q, lambda, P); 31 32 // Check diagonalization of M 33 for (int i=0; i<P; i++) 34 for (int j=0; j<P; j++) { 35 CeedScalar sum = 0; 36 for (int k=0; k<P; k++) 37 sum += Q[P*i+k]*lambda[k]*Q[P*j+k]; 38 Q_lambda_Qt[P*i+j] = sum; 39 } 40 for (int i=0; i<P; i++) 41 for (int j=0; j<P; j++) 42 if (fabs(M[P*i+j] - Q_lambda_Qt[P*i+j]) > 1E-14) 43 // LCOV_EXCL_START 44 printf("Error in diagonalization [%d, %d]: %f != %f\n", 45 i, j, M[P*i+j], Q_lambda_Qt[P*i+j]); 46 // LCOV_EXCL_STOP 47 48 CeedBasisDestroy(&basis); 49 CeedDestroy(&ceed); 50 return 0; 51 } 52