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
main(int argc,char ** argv)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