1 /***********************************************************************/
2 /* Open Visualization Data Explorer */
3 /* (C) Copyright IBM Corp. 1989,1999 */
4 /* ALL RIGHTS RESERVED */
5 /* This code licensed under the */
6 /* "IBM PUBLIC LICENSE - Open Visualization Data Explorer" */
7 /***********************************************************************/
8 /*
9 * eigen.c - routines to calculate eigenvectors and eigenvalues of a
10 * symmetric matrix
11 */
12 #include <string.h>
13 #include <dx/dx.h>
14 #include <math.h>
15 #include "eigen.h"
16
17 #include <dxconfig.h>
18
19
20 #define CYCLE(a,i,j,k,l) g=a[i][j]; h=a[k][l]; a[i][j]=g-s*(h+g*tau);\
21 a[k][l]=h+s*(g-h*tau);
22 Error
_dxfEigen(float ** a,int n,float d[],float ** v)23 _dxfEigen(float **a, int n, float d[], float **v)
24 {
25 float threshold, angle, tau, t, sum, s, h, g, c, *b, *z;
26 int j, Q, P, i, row, col;
27
28 /* check if matrix symmetric */
29 for (row = 1; row <= n; row++) {
30 for (col = 1; col <= n; col++) {
31 if (a[row][col] != a[col][row]) {
32 DXSetError (ERROR_DATA_INVALID, "A tensor must be a 3X3 symmetric matrix.");
33 return ERROR;
34 }
35 }
36 }
37
38 b = _dxfEigenVector(1, n);
39 z = _dxfEigenVector(1, n);
40
41 for (P = 1; P <= n; P++) {
42 for (Q = 1; Q <= n; Q++) v[P][Q]=0.0;
43 v[P][P] = 1.0;
44 }
45
46 for (P = 1; P <= n; P++) {
47 b[P] = d[P]= a[P][P];
48 z[P] = 0.0;
49 }
50
51 for (i = 1; i <= 50; i++) {
52 sum = 0.0;
53 for (P = 1; P <= n-1; P++) {
54 for (Q = P+1; Q<=n; Q++)
55 sum += fabs(a[P][Q]);
56 }
57 if (sum == 0.0) {
58 _dxfEigenFreeVector(z,1);
59 _dxfEigenFreeVector(b,1);
60 return OK;
61 }
62 if (i < 4)
63 threshold = 0.2 * sum / (n*n);
64 else
65 threshold = 0.0;
66 for (P = 1; P <= n-1; P++) {
67 for (Q = P+1; Q <= n; Q++) {
68 g=100.0*fabs(a[P][Q]);
69 if (i > 4 && (float)(fabs(d[P]) + g) == (float)fabs(d[P])
70 && (float)(fabs(d[Q]) + g) == (float)fabs(d[Q]))
71 a[P][Q]=0.0;
72 else if (fabs(a[P][Q]) > threshold) {
73 h = d[Q] - d[P];
74 if ((float)(fabs(h) + g) == (float)fabs(h))
75 t = (a[P][Q]) / h;
76 else {
77 angle = 0.5 * h / (a[P][Q]);
78 t = 1.0 / (fabs(angle) + sqrt(1.0 + angle*angle));
79 if (angle < 0.0) t = -t;
80 }
81 c = 1.0 / sqrt(1 + t * t);
82 s = t * c;
83 tau = s / (1.0 + c);
84 h = t * a[P][Q];
85 z[P] = z[P] - h;
86 z[Q] = z[Q] + h;
87 d[P] = d[P] - h;
88 d[Q] = d[Q] + h;
89 a[P][Q] = 0.0;
90 for (j = 1; j <= P-1; j++) {
91 CYCLE(a, j, P, j, Q)
92 }
93 for (j = P+1; j <= Q-1; j++) {
94 CYCLE(a, P, j, j, Q)
95 }
96 for (j = Q+1; j <= n; j++) {
97 CYCLE(a, P, j, Q, j)
98 }
99 for (j = 1; j <= n; j++) {
100 CYCLE(v, j, P, j, Q)
101 }
102 }
103 }
104 }
105 for (P = 1; P <= n; P++) {
106 b[P] = b[P] + z[P];
107 d[P] = b[P];
108 z[P] = 0.0;
109 }
110 }
111
112 /* evidently it shouldn't get here */
113 DXSetError (ERROR_DATA_INVALID, "Error calculating eigenvectors.");
114 return ERROR;
115 }
116 #undef CYCLE
117
_dxfEigenVector(int start_vector,int end_vector)118 float *_dxfEigenVector(int start_vector, int end_vector)
119 {
120 float *v;
121
122 /* Allocates a memory for a vector */
123
124 v=(float *)DXAllocate((unsigned) (end_vector-start_vector+1)*sizeof(float));
125 if (!v)
126 return NULL;
127 return v-start_vector;
128 }
129
_dxfEigenMatrix(int start_matrix_row,int end_matrix_row,int start_matrix_col,int end_matrix_col)130 float **_dxfEigenMatrix(int start_matrix_row, int end_matrix_row, int start_matrix_col, int end_matrix_col)
131 {
132 int i;
133 float **m;
134
135 /* Allocates memory for a matrix */
136
137 m=(float **) DXAllocate((unsigned) (end_matrix_row-start_matrix_row+1)*sizeof(float*));
138 if (!m)
139 return NULL;
140 m -= start_matrix_row;
141
142 for(i=start_matrix_row;i<=end_matrix_row;i++) {
143 m[i]=(float *) DXAllocate((unsigned) (end_matrix_col-start_matrix_col+1)*sizeof(float));
144 if (!m[i])
145 return NULL;
146 m[i] -= start_matrix_col;
147 }
148 return m;
149 }
150
_dxfEigenFreeVector(float * vec,int start_vector)151 void _dxfEigenFreeVector(float *vec, int start_vector)
152 {
153 /* fress memory Allocated by _dxfEigenVector */
154 DXFree((char*) (vec+start_vector));
155 }
156
_dxfEigenFreeMatrix(float ** m,int start_matrix_row,int end_matrix_row,int start_matrix_col)157 void _dxfEigenFreeMatrix(float **m, int start_matrix_row, int end_matrix_row, int start_matrix_col)
158 {
159 int i;
160
161 /* frees memory Allocated by _dxfEigenMatrix */
162
163 for(i=end_matrix_row;i>=start_matrix_row;i--)
164 DXFree((char*) (m[i]+start_matrix_col));
165 DXFree((char*) (m+start_matrix_row));
166 }
167
_dxfEigenConvertMatrix(float * a,int start_matrix_row,int end_matrix_row,int start_matrix_col,int end_matrix_col)168 float **_dxfEigenConvertMatrix(float *a, int start_matrix_row, int end_matrix_row, int start_matrix_col, int end_matrix_col)
169 {
170 int i,j,nrow,ncol;
171 float **m;
172
173 /*
174 * Allocates a float matrix
175 * m[start_matrix_row..end_matrix_row][start_matrix_col..end_matrix_col]
176 *
177 */
178
179 nrow=end_matrix_row-start_matrix_row+1;
180 ncol=end_matrix_col-start_matrix_col+1;
181 m = (float **) DXAllocate((unsigned) (nrow)*sizeof(float*));
182 if (!m)
183 return NULL;
184 m -= start_matrix_row;
185 for(i=0,j=start_matrix_row;i<=nrow-1;i++,j++)
186 m[j]=a+ncol*i-start_matrix_col;
187 return m;
188 }
189
_dxfEigenFreeConvertMatrix(float ** c,int start_matrix_row)190 void _dxfEigenFreeConvertMatrix(float**c, int start_matrix_row)
191 {
192 /* free the memory Allocated by _dxfEigenConvertMatrix */
193 DXFree((char*) (c+start_matrix_row));
194 }
195