13   CeedInt           p = 4, q = 4;  in main()  local
14   CeedScalar        M[p * p], K[p * p], X[p * p], lambda[p];  in main()
21   CeedBasisCreateTensorH1Lagrange(ceed, 1, 1, p, q, CEED_GAUSS, &basis);  in main()
25   for (int i = 0; i < p; i++) {  in main()
26     for (int j = 0; j < p; j++) {  in main()
29         sum_m += interpolation[p * k + i] * quadrature_weights[k] * interpolation[p * k + j];  in main()
30         sum_k += gradient[p * k + i] * quadrature_weights[k] * gradient[p * k + j];  in main()
32       M[p * i + j] = sum_m;  in main()
33       K[p * i + j] = sum_k;  in main()
37   CeedSimultaneousDiagonalization(ceed, K, M, X, lambda, p);  in main()
40   CeedScalar work_array[p * p];  in main()
41   for (int i = 0; i < p; i++) {  in main()
42     for (int j = 0; j < p; j++) {  in main()
44       for (int k = 0; k < p; k++) sum += M[p * i + k] * X[p * k + j];  in main()
45       work_array[p * i + j] = sum;  in main()
48   for (int i = 0; i < p; i++) {  in main()
49     for (int j = 0; j < p; j++) {  in main()
51       for (int k = 0; k < p; k++) sum += X[p * k + i] * work_array[p * k + j];  in main()
52       M[p * i + j] = sum;  in main()
55   for (int i = 0; i < p; i++) {  in main()
56     for (int j = 0; j < p; j++) {  in main()
57       if (fabs(M[p * i + j] - (i == j ? 1.0 : 0.0)) > 100. * CEED_EPSILON) {  in main()
59 …onalization of M [%" CeedInt_FMT ", %" CeedInt_FMT "]: %f != %f\n", i, j, M[p * i + j], (i == j ? …  in main()
66   for (int i = 0; i < p; i++) {  in main()
67     for (int j = 0; j < p; j++) {  in main()
69       for (int k = 0; k < p; k++) sum += K[p * i + k] * X[p * k + j];  in main()
70       work_array[p * i + j] = sum;  in main()
73   for (int i = 0; i < p; i++) {  in main()
74     for (int j = 0; j < p; j++) {  in main()
76       for (int k = 0; k < p; k++) sum += X[p * k + i] * work_array[p * k + j];  in main()
77       K[p * i + j] = sum;  in main()
80   for (int i = 0; i < p; i++) {  in main()
81     for (int j = 0; j < p; j++) {  in main()
82       if (fabs(K[p * i + j] - (i == j ? lambda[i] : 0.0)) > 1e3 * CEED_EPSILON) {  in main()
84 …onalization of K [%" CeedInt_FMT ", %" CeedInt_FMT "]: %f != %f\n", i, j, K[p * i + j], (i == j ? …  in main()