1 /// @file 2 /// Test Simultaneous Diagonalization 3 /// \test Simultaneous Diagonalization 4 #include <ceed.h> 5 6 int main(int argc, char **argv) { 7 Ceed ceed; 8 CeedScalar M[16] = {0.2, 0.0745355993, -0.0745355993, 0.0333333333, 9 0.0745355993, 1., 0.1666666667, -0.0745355993, 10 -0.0745355993, 0.1666666667, 1., 0.0745355993, 11 0.0333333333, -0.0745355993, 0.0745355993, 0.2 12 }; 13 CeedScalar K[16] = {3.0333333333, -3.4148928136, 0.4982261470, -0.1166666667, 14 -3.4148928136, 5.8333333333, -2.9166666667, 0.4982261470, 15 0.4982261470, -2.9166666667, 5.8333333333, -3.4148928136, 16 -0.1166666667, 0.4982261470, -3.4148928136, 3.0333333333 17 }; 18 CeedScalar x[16], lambda[4], xxt[16]; 19 20 CeedInit(argv[1], &ceed); 21 22 CeedSimultaneousDiagonalization(ceed, K, M, x, lambda, 4); 23 24 for (int i=0; i<4; i++) 25 for (int j=0; j<4; j++) { 26 xxt[j+4*i] = 0; 27 for (int k=0; k<4; k++) 28 xxt[j+4*i] += x[k+4*i]*x[k+4*j]; 29 } 30 31 fprintf(stdout, "x x^T:\n"); 32 for (int i=0; i<4; i++) { 33 for (int j=0; j<4; j++) { 34 if (xxt[j+4*i] <= 1E-14 && xxt[j+4*i] >= -1E-14) xxt[j+4*i] = 0; 35 fprintf(stdout, "%12.8f\t", xxt[j+4*i]); 36 } 37 fprintf(stdout, "\n"); 38 } 39 fprintf(stdout, "lambda:\n"); 40 for (int i=0; i<4; i++) { 41 if (lambda[i] <= 1E-14 && lambda[i] >= -1E-14) lambda[i] = 0; 42 fprintf(stdout, "%12.8f\n", lambda[i]); 43 } 44 CeedDestroy(&ceed); 45 return 0; 46 } 47