1 //
2 // matrix_eig_test.c
3 //
4 // Test eigenvalue computation using Q/R decomposition. Eigenvalues
5 // will be on diagonal if they are all real.
6 //
7
8 #include <stdio.h>
9 #include <string.h>
10 #include "liquid.h"
11
main()12 int main() {
13 unsigned int num_iterations=8;
14
15 #if 0
16 unsigned int n=4;
17 float A[16] = {
18 1,2,3,4,
19 5,5,7,8,
20 6,4,8,7,
21 1,0,3,1};
22 #else
23 unsigned int n=3;
24 float A[9] = {
25 1.0f, 1.0f, 1.0f,
26 1.0f, 2.0f, 1.0f,
27 1.0f, 1.0f, 2.0f};
28 #endif
29
30 float eig[n];
31 float A0[n*n];
32 float Q[n*n], R[n*n];
33
34 memmove(A0, A, n*n*sizeof(float));
35
36 printf("\n");
37 printf("testing Q/R decomposition [Gram-Schmidt]\n");
38 matrixf_qrdecomp_gramschmidt(A,n,n,Q,R);
39 //matrixf_print(Q,n,n);
40 //matrixf_print(R,n,n);
41
42 unsigned int k;
43 for (k=0; k<num_iterations; k++) {
44 // compute Q/R decomposition
45 matrixf_qrdecomp_gramschmidt(A0,n,n,Q,R);
46
47 // compute A1 = R*Q
48 matrixf_mul(R,n,n, Q,n,n, A0,n,n);
49
50 matrixf_print(A0,n,n);
51 }
52
53 // extract eigen values
54 unsigned int i;
55 for (i=0; i<n; i++)
56 eig[i] = matrix_access(A0,n,n,i,i);
57
58 for (i=0; i<n; i++)
59 printf("eig[%3u] = %12.8f\n", i, eig[i]);
60
61 printf("done.\n");
62 return 0;
63 }
64
65