xref: /libCEED/tests/t305-basis.c (revision fb5510374522664021f5faa53a016cd7b340b623)
14411cf47Sjeremylt /// @file
252bfb9bbSJeremy L Thompson /// Test Simultaneous Diagonalization
352bfb9bbSJeremy L Thompson /// \test Simultaneous Diagonalization
457c64913Sjeremylt #include <ceed.h>
557c64913Sjeremylt 
657c64913Sjeremylt int main(int argc, char **argv) {
757c64913Sjeremylt   Ceed ceed;
8*fb551037Sjeremylt   CeedScalar M[16] = {0.2, 0.0745355993, -0.0745355993, 0.0333333333,
9*fb551037Sjeremylt                       0.0745355993, 1., 0.1666666667, -0.0745355993,
10*fb551037Sjeremylt                       -0.0745355993, 0.1666666667, 1., 0.0745355993,
11*fb551037Sjeremylt                       0.0333333333, -0.0745355993, 0.0745355993, 0.2
127f823360Sjeremylt                      };
13*fb551037Sjeremylt   CeedScalar K[16] = {3.0333333333, -3.4148928136, 0.4982261470, -0.1166666667,
14*fb551037Sjeremylt                       -3.4148928136, 5.8333333333, -2.9166666667, 0.4982261470,
15*fb551037Sjeremylt                       0.4982261470, -2.9166666667, 5.8333333333, -3.4148928136,
16*fb551037Sjeremylt                       -0.1166666667, 0.4982261470, -3.4148928136, 3.0333333333
177f823360Sjeremylt                      };
18*fb551037Sjeremylt   CeedScalar x[16], lambda[4], xxt[16];
1957c64913Sjeremylt 
2057c64913Sjeremylt   CeedInit(argv[1], &ceed);
21aedaa0e5Sjeremylt 
2252bfb9bbSJeremy L Thompson   CeedSimultaneousDiagonalization(ceed, K, M, x, lambda, 4);
23*fb551037Sjeremylt 
24*fb551037Sjeremylt   for (int i=0; i<4; i++)
25*fb551037Sjeremylt     for (int j=0; j<4; j++) {
26*fb551037Sjeremylt       xxt[j+4*i] = 0;
27*fb551037Sjeremylt       for (int k=0; k<4; k++)
28*fb551037Sjeremylt         xxt[j+4*i] += x[k+4*i]*x[k+4*j];
29*fb551037Sjeremylt     }
30*fb551037Sjeremylt 
31*fb551037Sjeremylt   fprintf(stdout, "x x^T:\n");
3252bfb9bbSJeremy L Thompson   for (int i=0; i<4; i++) {
3352bfb9bbSJeremy L Thompson     for (int j=0; j<4; j++) {
34*fb551037Sjeremylt       if (xxt[j+4*i] <= 1E-14 && xxt[j+4*i] >= -1E-14) xxt[j+4*i] = 0;
35*fb551037Sjeremylt       fprintf(stdout, "%12.8f\t", xxt[j+4*i]);
3652bfb9bbSJeremy L Thompson     }
3752bfb9bbSJeremy L Thompson     fprintf(stdout, "\n");
3852bfb9bbSJeremy L Thompson   }
3952bfb9bbSJeremy L Thompson   fprintf(stdout, "lambda:\n");
4052bfb9bbSJeremy L Thompson   for (int i=0; i<4; i++) {
4152bfb9bbSJeremy L Thompson     if (lambda[i] <= 1E-14 && lambda[i] >= -1E-14) lambda[i] = 0;
4252bfb9bbSJeremy L Thompson     fprintf(stdout, "%12.8f\n", lambda[i]);
4352bfb9bbSJeremy L Thompson   }
4457c64913Sjeremylt   CeedDestroy(&ceed);
4557c64913Sjeremylt   return 0;
4657c64913Sjeremylt }
47