| 1237 | } |
| 1238 | |
| 1239 | void eigenVector_d(float cov[MAX_DIMENSION_BIG][MAX_DIMENSION_BIG], float vector[MAX_DIMENSION_BIG], int dimension) { |
| 1240 | // calculate an eigenvecto corresponding to a biggest eigenvalue |
| 1241 | // will work for non-zero non-negative matricies only |
| 1242 | |
| 1243 | #define EV_ITERATION_NUMBER 20 |
| 1244 | #define EV_SLACK 2 /* additive for exp base 2)*/ |
| 1245 | |
| 1246 | |
| 1247 | int i, j, k, l, m, n, p, q; |
| 1248 | float c[2][MAX_DIMENSION_BIG][MAX_DIMENSION_BIG]; |
| 1249 | float maxDiag; |
| 1250 | |
| 1251 | for (i = 0; i<dimension; i++) |
| 1252 | for (j = 0; j<dimension; j++) |
| 1253 | c[0][i][j] = cov[i][j]; |
| 1254 | |
| 1255 | p = (int)floorf(logf((HDR_FLT_MAX_EXP - EV_SLACK) / ceilf(logf((float)dimension) / logf(2.0f))) / logf(2.0f)); |
| 1256 | |
| 1257 | //assert(p>0); |
| 1258 | |
| 1259 | p = p >0 ? p : 1; |
| 1260 | |
| 1261 | q = (EV_ITERATION_NUMBER + p - 1) / p; |
| 1262 | |
| 1263 | l = 0; |
| 1264 | |
| 1265 | for (n = 0; n<q; n++) { |
| 1266 | maxDiag = 0; |
| 1267 | |
| 1268 | for (i = 0; i<dimension; i++) |
| 1269 | maxDiag = c[l][i][i] > maxDiag ? c[l][i][i] : maxDiag; |
| 1270 | |
| 1271 | if (maxDiag <= 0) { |
| 1272 | return; |
| 1273 | } |
| 1274 | |
| 1275 | //assert(maxDiag >0); |
| 1276 | |
| 1277 | for (i = 0; i<dimension; i++) |
| 1278 | for (j = 0; j<dimension; j++) |
| 1279 | c[l][i][j] /= maxDiag; |
| 1280 | |
| 1281 | for (m = 0; m<p; m++) { |
| 1282 | for (i = 0; i<dimension; i++) |
| 1283 | for (j = 0; j<dimension; j++) { |
| 1284 | float temp = 0; |
| 1285 | for (k = 0; k<dimension; k++) { |
| 1286 | // Notes: |
| 1287 | // This is the most consuming portion of the code and needs optimizing for perfromance |
| 1288 | temp += c[l][i][k] * c[l][k][j]; |
| 1289 | } |
| 1290 | c[1 - l][i][j] = temp; |
| 1291 | } |
| 1292 | l = 1 - l; |
| 1293 | } |
| 1294 | } |
| 1295 | |
| 1296 | maxDiag = 0; |