eigen.cpp

   1 //-----------------------------------------------------------------------------
   2 //
   3 //      Compute eigenvalues and eigenvectors
   4 //
   5 //      Input:          stack[tos - 1]          symmetric matrix
   6 //
   7 //      Output:         D                       diagnonal matrix
   8 //
   9 //                      Q                       eigenvector matrix
  10 //
  11 //      D and Q have the property that
  12 //
  13 //              A == dot(transpose(Q),D,Q)
  14 //
  15 //      where A is the original matrix.
  16 //
  17 //      The eigenvalues are on the diagonal of D.
  18 //
  19 //      The eigenvectors are row vectors in Q.
  20 //
  21 //      The eigenvalue relation
  22 //
  23 //              A X = lambda X
  24 //
  25 //      can be checked as follows:
  26 //
  27 //              lambda = D[1,1]
  28 //
  29 //              X = Q[1]
  30 //
  31 //              dot(A,X) - lambda X
  32 //
  33 //-----------------------------------------------------------------------------
  34
  35 #include "stdafx.h"
  36 #include "defs.h"
  37
  38 #define D(i, j) yydd[n * (i) + (j)]
  39 #define Q(i, j) yyqq[n * (i) + (j)]
  40
  41 extern void copy_tensor(void);
  42 static void eigen(int);
  43 static int check_arg(void);
  44 static int step(void);
  45 static void step2(int, int);
  46 static int n;
  47 static double *yydd, *yyqq;
  48
  49 void
  50 eval_eigen(void)
  51 {
  52         if (check_arg() == 0)
  53                 stop("eigen: argument is not a square matrix");
  54
  55         eigen(EIGEN);
  56
  57         p1 = usr_symbol("D");
  58         set_binding(p1, p2);
  59
  60         p1 = usr_symbol("Q");
  61         set_binding(p1, p3);
  62
  63         push(symbol(NIL));
  64 }
  65
  66 void
  67 eval_eigenval(void)
  68 {
  69         if (check_arg() == 0) {
  70                 push_symbol(EIGENVAL);
  71                 push(p1);
  72                 list(2);
  73                 return;
  74         }
  75
  76         eigen(EIGENVAL);
  77
  78         push(p2);
  79 }
  80
  81 void
  82 eval_eigenvec(void)
  83 {
  84         if (check_arg() == 0) {
  85                 push_symbol(EIGENVEC);
  86                 push(p1);
  87                 list(2);
  88                 return;
  89         }
  90
  91         eigen(EIGENVEC);
  92
  93         push(p3);
  94 }
  95
  96 static int
  97 check_arg(void)
  98 {
  99         int i, j;
 100
 101         push(cadr(p1));
 102         eval();
 103         yyfloat();
 104         eval();
 105         p1 = pop();
 106
 107         if (!istensor(p1))
 108                 return 0;
 109
 110         if (p1->u.tensor->ndim != 2 || p1->u.tensor->dim[0] != p1->u.tensor->dim[1])
 111                 stop("eigen: argument is not a square matrix");
 112
 113         n = p1->u.tensor->dim[0];
 114
 115         for (i = 0; i < n; i++)
 116                 for (j = 0; j < n; j++)
 117                         if (!isdouble(p1->u.tensor->elem[n * i + j]))
 118                                 stop("eigen: matrix is not numerical");
 119
 120         for (i = 0; i < n - 1; i++)
 121                 for (j = i + 1; j < n; j++)
 122                         if (fabs(p1->u.tensor->elem[n * i + j]->u.d - p1->u.tensor->elem[n * j + i]->u.d) > 1e-10)
 123                                 stop("eigen: matrix is not symmetrical");
 124
 125         return 1;
 126 }
 127
 128 //-----------------------------------------------------------------------------
 129 //
 130 //      Input:          p1              matrix
 131 //
 132 //      Output:         p2              eigenvalues
 133 //
 134 //                      p3              eigenvectors
 135 //
 136 //-----------------------------------------------------------------------------
 137
 138 static void
 139 eigen(int op)
 140 {
 141         int i, j;
 142
 143         // malloc working vars
 144
 145         yydd = (double *) malloc(n * n * sizeof (double));
 146
 147         if (yydd == NULL)
 148                 stop("malloc failure");
 149
 150         yyqq = (double *) malloc(n * n * sizeof (double));
 151
 152         if (yyqq == NULL)
 153                 stop("malloc failure");
 154
 155         // initialize D
 156
 157         for (i = 0; i < n; i++) {
 158                 D(i, i) = p1->u.tensor->elem[n * i + i]->u.d;
 159                 for (j = i + 1; j < n; j++) {
 160                         D(i, j) = p1->u.tensor->elem[n * i + j]->u.d;
 161                         D(j, i) = p1->u.tensor->elem[n * i + j]->u.d;
 162                 }
 163         }
 164
 165         // initialize Q
 166
 167         for (i = 0; i < n; i++) {
 168                 Q(i, i) = 1.0;
 169                 for (j = i + 1; j < n; j++) {
 170                         Q(i, j) = 0.0;
 171                         Q(j, i) = 0.0;
 172                 }
 173         }
 174
 175         // step up to 100 times
 176
 177         for (i = 0; i < 100; i++)
 178                 if (step() == 0)
 179                         break;
 180
 181         if (i == 100)
 182                 printstr("\nnote: eigen did not converge");
 183
 184         // p2 = D
 185
 186         if (op == EIGEN || op == EIGENVAL) {
 187
 188                 push(p1);
 189                 copy_tensor();
 190                 p2 = pop();
 191
 192                 for (i = 0; i < n; i++) {
 193                         for (j = 0; j < n; j++) {
 194                                 push_double(D(i, j));
 195                                 p2->u.tensor->elem[n * i + j] = pop();
 196                         }
 197                 }
 198         }
 199
 200         // p3 = Q
 201
 202         if (op == EIGEN || op == EIGENVEC) {
 203
 204                 push(p1);
 205                 copy_tensor();
 206                 p3 = pop();
 207
 208                 for (i = 0; i < n; i++) {
 209                         for (j = 0; j < n; j++) {
 210                                 push_double(Q(i, j));
 211                                 p3->u.tensor->elem[n * i + j] = pop();
 212                         }
 213                 }
 214         }
 215
 216         // free working vars
 217
 218         free(yydd);
 219         free(yyqq);
 220 }
 221
 222 //-----------------------------------------------------------------------------
 223 //
 224 //      Example: p = 1, q = 3
 225 //
 226 //              c       0       s       0
 227 //
 228 //              0       1       0       0
 229 //      G =
 230 //              -s      0       c       0
 231 //
 232 //              0       0       0       1
 233 //
 234 //      The effect of multiplying G times A is...
 235 //
 236 //      row 1 of A    = c (row 1 of A ) + s (row 3 of A )
 237 //                n+1                n                 n
 238 //
 239 //      row 3 of A    = c (row 3 of A ) - s (row 1 of A )
 240 //                n+1                n                 n
 241 //
 242 //      In terms of components the overall effect is...
 243 //
 244 //      row 1 = c row 1 + s row 3
 245 //
 246 //              A[1,1] = c A[1,1] + s A[3,1]
 247 //
 248 //              A[1,2] = c A[1,2] + s A[3,2]
 249 //
 250 //              A[1,3] = c A[1,3] + s A[3,3]
 251 //
 252 //              A[1,4] = c A[1,4] + s A[3,4]
 253 //
 254 //      row 3 = c row 3 - s row 1
 255 //
 256 //              A[3,1] = c A[3,1] - s A[1,1]
 257 //
 258 //              A[3,2] = c A[3,2] - s A[1,2]
 259 //
 260 //              A[3,3] = c A[3,3] - s A[1,3]
 261 //
 262 //              A[3,4] = c A[3,4] - s A[1,4]
 263 //
 264 //                                         T
 265 //      The effect of multiplying A times G  is...
 266 //
 267 //      col 1 of A    = c (col 1 of A ) + s (col 3 of A )
 268 //                n+1                n                 n
 269 //
 270 //      col 3 of A    = c (col 3 of A ) - s (col 1 of A )
 271 //                n+1                n                 n
 272 //
 273 //      In terms of components the overall effect is...
 274 //
 275 //      col 1 = c col 1 + s col 3
 276 //
 277 //              A[1,1] = c A[1,1] + s A[1,3]
 278 //
 279 //              A[2,1] = c A[2,1] + s A[2,3]
 280 //
 281 //              A[3,1] = c A[3,1] + s A[3,3]
 282 //
 283 //              A[4,1] = c A[4,1] + s A[4,3]
 284 //
 285 //      col 3 = c col 3 - s col 1
 286 //
 287 //              A[1,3] = c A[1,3] - s A[1,1]
 288 //
 289 //              A[2,3] = c A[2,3] - s A[2,1]
 290 //
 291 //              A[3,3] = c A[3,3] - s A[3,1]
 292 //
 293 //              A[4,3] = c A[4,3] - s A[4,1]
 294 //
 295 //      What we want to do is just compute the upper triangle of A since we
 296 //      know the lower triangle is identical.
 297 //
 298 //      In other words, we just want to update components A[i,j] where i < j.
 299 //
 300 //-----------------------------------------------------------------------------
 301 //
 302 //      Example: p = 2, q = 5
 303 //
 304 //                              p                       q
 305 //
 306 //                      j=1     j=2     j=3     j=4     j=5     j=6
 307 //
 308 //              i=1     .       A[1,2]  .       .       A[1,5]  .
 309 //
 310 //      p       i=2     A[2,1]  A[2,2]  A[2,3]  A[2,4]  A[2,5]  A[2,6]
 311 //
 312 //              i=3     .       A[3,2]  .       .       A[3,5]  .
 313 //
 314 //              i=4     .       A[4,2]  .       .       A[4,5]  .
 315 //
 316 //      q       i=5     A[5,1]  A[5,2]  A[5,3]  A[5,4]  A[5,5]  A[5,6]
 317 //
 318 //              i=6     .       A[6,2]  .       .       A[6,5]  .
 319 //
 320 //-----------------------------------------------------------------------------
 321 //
 322 //      This is what B = GA does:
 323 //
 324 //      row 2 = c row 2 + s row 5
 325 //
 326 //              B[2,1] = c * A[2,1] + s * A[5,1]
 327 //              B[2,2] = c * A[2,2] + s * A[5,2]
 328 //              B[2,3] = c * A[2,3] + s * A[5,3]
 329 //              B[2,4] = c * A[2,4] + s * A[5,4]
 330 //              B[2,5] = c * A[2,5] + s * A[5,5]
 331 //              B[2,6] = c * A[2,6] + s * A[5,6]
 332 //
 333 //      row 5 = c row 5 - s row 2
 334 //
 335 //              B[5,1] = c * A[5,1] + s * A[2,1]
 336 //              B[5,2] = c * A[5,2] + s * A[2,2]
 337 //              B[5,3] = c * A[5,3] + s * A[2,3]
 338 //              B[5,4] = c * A[5,4] + s * A[2,4]
 339 //              B[5,5] = c * A[5,5] + s * A[2,5]
 340 //              B[5,6] = c * A[5,6] + s * A[2,6]
 341 //
 342 //                     T
 343 //      This is what BG  does:
 344 //
 345 //      col 2 = c col 2 + s col 5
 346 //
 347 //              B[1,2] = c * A[1,2] + s * A[1,5]
 348 //              B[2,2] = c * A[2,2] + s * A[2,5]
 349 //              B[3,2] = c * A[3,2] + s * A[3,5]
 350 //              B[4,2] = c * A[4,2] + s * A[4,5]
 351 //              B[5,2] = c * A[5,2] + s * A[5,5]
 352 //              B[6,2] = c * A[6,2] + s * A[6,5]
 353 //
 354 //      col 5 = c col 5 - s col 2
 355 //
 356 //              B[1,5] = c * A[1,5] - s * A[1,2]
 357 //              B[2,5] = c * A[2,5] - s * A[2,2]
 358 //              B[3,5] = c * A[3,5] - s * A[3,2]
 359 //              B[4,5] = c * A[4,5] - s * A[4,2]
 360 //              B[5,5] = c * A[5,5] - s * A[5,2]
 361 //              B[6,5] = c * A[6,5] - s * A[6,2]
 362 //
 363 //-----------------------------------------------------------------------------
 364 //
 365 //      Step 1: Just do upper triangle (i < j), B[2,5] = 0
 366 //
 367 //              B[1,2] = c * A[1,2] + s * A[1,5]
 368 //
 369 //              B[2,3] = c * A[2,3] + s * A[5,3]
 370 //              B[2,4] = c * A[2,4] + s * A[5,4]
 371 //              B[2,6] = c * A[2,6] + s * A[5,6]
 372 //
 373 //              B[1,5] = c * A[1,5] - s * A[1,2]
 374 //              B[3,5] = c * A[3,5] - s * A[3,2]
 375 //              B[4,5] = c * A[4,5] - s * A[4,2]
 376 //
 377 //              B[5,6] = c * A[5,6] + s * A[2,6]
 378 //
 379 //-----------------------------------------------------------------------------
 380 //
 381 //      Step 2: Transpose where i > j since A[i,j] == A[j,i]
 382 //
 383 //              B[1,2] = c * A[1,2] + s * A[1,5]
 384 //
 385 //              B[2,3] = c * A[2,3] + s * A[3,5]
 386 //              B[2,4] = c * A[2,4] + s * A[4,5]
 387 //              B[2,6] = c * A[2,6] + s * A[5,6]
 388 //
 389 //              B[1,5] = c * A[1,5] - s * A[1,2]
 390 //              B[3,5] = c * A[3,5] - s * A[2,3]
 391 //              B[4,5] = c * A[4,5] - s * A[2,4]
 392 //
 393 //              B[5,6] = c * A[5,6] + s * A[2,6]
 394 //
 395 //-----------------------------------------------------------------------------
 396 //
 397 //      Step 3: Same as above except reorder
 398 //
 399 //      k < p           (k = 1)
 400 //
 401 //              A[1,2] = c * A[1,2] + s * A[1,5]
 402 //              A[1,5] = c * A[1,5] - s * A[1,2]
 403 //
 404 //      p < k < q       (k = 3..4)
 405 //
 406 //              A[2,3] = c * A[2,3] + s * A[3,5]
 407 //              A[3,5] = c * A[3,5] - s * A[2,3]
 408 //
 409 //              A[2,4] = c * A[2,4] + s * A[4,5]
 410 //              A[4,5] = c * A[4,5] - s * A[2,4]
 411 //
 412 //      q < k           (k = 6)
 413 //
 414 //              A[2,6] = c * A[2,6] + s * A[5,6]
 415 //              A[5,6] = c * A[5,6] - s * A[2,6]
 416 //
 417 //-----------------------------------------------------------------------------
 418
 419 static int
 420 step(void)
 421 {
 422         int count, i, j;
 423
 424         count = 0;
 425
 426         // for each upper triangle "off-diagonal" component do step2
 427
 428         for (i = 0; i < n - 1; i++) {
 429                 for (j = i + 1; j < n; j++) {
 430                         if (D(i, j) != 0.0) {
 431                                 step2(i, j);
 432                                 count++;
 433                         }
 434                 }
 435         }
 436
 437         return count;
 438 }
 439
 440 static void
 441 step2(int p, int q)
 442 {
 443         int k;
 444         double t, theta;
 445         double c, cc, s, ss;
 446
 447         // compute c and s
 448
 449         // from Numerical Recipes (except they have a_qq - a_pp)
 450
 451         theta = 0.5 * (D(p, p) - D(q, q)) / D(p, q);
 452
 453         t = 1.0 / (fabs(theta) + sqrt(theta * theta + 1.0));
 454
 455         if (theta < 0.0)
 456                 t = -t;
 457
 458         c = 1.0 / sqrt(t * t + 1.0);
 459
 460         s = t * c;
 461
 462         // D = GD
 463
 464         // which means "add rows"
 465
 466         for (k = 0; k < n; k++) {
 467                 cc = D(p, k);
 468                 ss = D(q, k);
 469                 D(p, k) = c * cc + s * ss;
 470                 D(q, k) = c * ss - s * cc;
 471         }
 472
 473         // D = D transpose(G)
 474
 475         // which means "add columns"
 476
 477         for (k = 0; k < n; k++) {
 478                 cc = D(k, p);
 479                 ss = D(k, q);
 480                 D(k, p) = c * cc + s * ss;
 481                 D(k, q) = c * ss - s * cc;
 482         }
 483
 484         // Q = GQ
 485
 486         // which means "add rows"
 487
 488         for (k = 0; k < n; k++) {
 489                 cc = Q(p, k);
 490                 ss = Q(q, k);
 491                 Q(p, k) = c * cc + s * ss;
 492                 Q(q, k) = c * ss - s * cc;
 493         }
 494
 495         D(p, q) = 0.0;
 496         D(q, p) = 0.0;
 497 }