1 /// @file
2 /// Test Simultaneous Diagonalization
3 /// \test Simultaneous Diagonalization
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, q = 4;
14 CeedScalar M[p * p], K[p * p], X[p * p], lambda[p];
15 CeedBasis basis;
16 const CeedScalar *interpolation, *gradient, *quadrature_weights;
17
18 CeedInit(argv[1], &ceed);
19
20 // Create mass, stiffness matrix
21 CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis);
22 CeedBasisGetInterp(basis, &interpolation);
23 CeedBasisGetGrad(basis, &gradient);
24 CeedBasisGetQWeights(basis, &quadrature_weights);
25 for (int i = 0; i < p; i++) {
26 for (int j = 0; j < p; j++) {
27 CeedScalar sum_m = 0, sum_k = 0;
28 for (int k = 0; k < q; k++) {
29 sum_m += interpolation[p * k + i] * quadrature_weights[k] * interpolation[p * k + j];
30 sum_k += gradient[p * k + i] * quadrature_weights[k] * gradient[p * k + j];
31 }
32 M[p * i + j] = sum_m;
33 K[p * i + j] = sum_k;
34 }
35 }
36
37 CeedSimultaneousDiagonalization(ceed, K, M, X, lambda, p);
38
39 // Check X^T M X = I
40 CeedScalar work_array[p * p];
41 for (int i = 0; i < p; i++) {
42 for (int j = 0; j < p; j++) {
43 CeedScalar sum = 0;
44 for (int k = 0; k < p; k++) sum += M[p * i + k] * X[p * k + j];
45 work_array[p * i + j] = sum;
46 }
47 }
48 for (int i = 0; i < p; i++) {
49 for (int j = 0; j < p; j++) {
50 CeedScalar sum = 0;
51 for (int k = 0; k < p; k++) sum += X[p * k + i] * work_array[p * k + j];
52 M[p * i + j] = sum;
53 }
54 }
55 for (int i = 0; i < p; i++) {
56 for (int j = 0; j < p; j++) {
57 if (fabs(M[p * i + j] - (i == j ? 1.0 : 0.0)) > 100. * CEED_EPSILON) {
58 // LCOV_EXCL_START
59 printf("Error in diagonalization of M [%" CeedInt_FMT ", %" CeedInt_FMT "]: %f != %f\n", i, j, M[p * i + j], (i == j ? 1.0 : 0.0));
60 // LCOV_EXCL_STOP
61 }
62 }
63 }
64
65 // Check X^T K X = Lambda
66 for (int i = 0; i < p; i++) {
67 for (int j = 0; j < p; j++) {
68 CeedScalar sum = 0;
69 for (int k = 0; k < p; k++) sum += K[p * i + k] * X[p * k + j];
70 work_array[p * i + j] = sum;
71 }
72 }
73 for (int i = 0; i < p; i++) {
74 for (int j = 0; j < p; j++) {
75 CeedScalar sum = 0;
76 for (int k = 0; k < p; k++) sum += X[p * k + i] * work_array[p * k + j];
77 K[p * i + j] = sum;
78 }
79 }
80 for (int i = 0; i < p; i++) {
81 for (int j = 0; j < p; j++) {
82 if (fabs(K[p * i + j] - (i == j ? lambda[i] : 0.0)) > 1e3 * CEED_EPSILON) {
83 // LCOV_EXCL_START
84 printf("Error in diagonalization of K [%" CeedInt_FMT ", %" CeedInt_FMT "]: %f != %f\n", i, j, K[p * i + j], (i == j ? lambda[i] : 0.0));
85 // LCOV_EXCL_STOP
86 }
87 }
88 }
89
90 CeedBasisDestroy(&basis);
91 CeedDestroy(&ceed);
92 return 0;
93 }
94