tensor.cpp

   1 #include "stdafx.h"
   2
   3 #include "defs.h"
   4
   5 static void promote_tensor(void);
   6 static int compatible(U *, U *);
   7
   8 //-----------------------------------------------------------------------------
   9 //
  10 //      Called from the "eval" module to evaluate tensor elements.
  11 //
  12 //      p1 points to the tensor operand.
  13 //
  14 //-----------------------------------------------------------------------------
  15
  16 void
  17 eval_tensor(void)
  18 {
  19         int i, ndim, nelem;
  20         U **a, **b;
  21
  22         //---------------------------------------------------------------------
  23         //
  24         //      create a new tensor for the result
  25         //
  26         //---------------------------------------------------------------------
  27
  28         nelem = p1->u.tensor->nelem;
  29
  30         ndim = p1->u.tensor->ndim;
  31
  32         p2 = alloc_tensor(nelem);
  33
  34         p2->u.tensor->ndim = ndim;
  35
  36         for (i = 0; i < ndim; i++)
  37                 p2->u.tensor->dim[i] = p1->u.tensor->dim[i];
  38
  39         //---------------------------------------------------------------------
  40         //
  41         //      b = eval(a)
  42         //
  43         //---------------------------------------------------------------------
  44
  45         a = p1->u.tensor->elem;
  46         b = p2->u.tensor->elem;
  47
  48         for (i = 0; i < nelem; i++) {
  49                 push(a[i]);
  50                 eval();
  51                 b[i] = pop();
  52         }
  53
  54         //---------------------------------------------------------------------
  55         //
  56         //      push the result
  57         //
  58         //---------------------------------------------------------------------
  59
  60         push(p2);
  61
  62         promote_tensor();
  63 }
  64
  65 //-----------------------------------------------------------------------------
  66 //
  67 //      Add tensors
  68 //
  69 //      Input:          Operands on stack
  70 //
  71 //      Output:         Result on stack
  72 //
  73 //-----------------------------------------------------------------------------
  74
  75 void
  76 tensor_plus_tensor(void)
  77 {
  78         int i, ndim, nelem;
  79         U **a, **b, **c;
  80
  81         save();
  82
  83         p2 = pop();
  84         p1 = pop();
  85
  86         // are the dimension lists equal?
  87
  88         ndim = p1->u.tensor->ndim;
  89
  90         if (ndim != p2->u.tensor->ndim) {
  91                 push(symbol(NIL));
  92                 restore();
  93                 return;
  94         }
  95
  96         for (i = 0; i < ndim; i++)
  97                 if (p1->u.tensor->dim[i] != p2->u.tensor->dim[i]) {
  98                         push(symbol(NIL));
  99                         restore();
 100                         return;
 101                 }
 102
 103         // create a new tensor for the result
 104
 105         nelem = p1->u.tensor->nelem;
 106
 107         p3 = alloc_tensor(nelem);
 108
 109         p3->u.tensor->ndim = ndim;
 110
 111         for (i = 0; i < ndim; i++)
 112                 p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
 113
 114         // c = a + b
 115
 116         a = p1->u.tensor->elem;
 117         b = p2->u.tensor->elem;
 118         c = p3->u.tensor->elem;
 119
 120         for (i = 0; i < nelem; i++) {
 121                 push(a[i]);
 122                 push(b[i]);
 123                 add();
 124                 c[i] = pop();
 125         }
 126
 127         // push the result
 128
 129         push(p3);
 130
 131         restore();
 132 }
 133
 134 //-----------------------------------------------------------------------------
 135 //
 136 //      careful not to reorder factors
 137 //
 138 //-----------------------------------------------------------------------------
 139
 140 void
 141 tensor_times_scalar(void)
 142 {
 143         int i, ndim, nelem;
 144         U **a, **b;
 145
 146         save();
 147
 148         p2 = pop();
 149         p1 = pop();
 150
 151         ndim = p1->u.tensor->ndim;
 152         nelem = p1->u.tensor->nelem;
 153
 154         p3 = alloc_tensor(nelem);
 155
 156         p3->u.tensor->ndim = ndim;
 157
 158         for (i = 0; i < ndim; i++)
 159                 p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
 160
 161         a = p1->u.tensor->elem;
 162         b = p3->u.tensor->elem;
 163
 164         for (i = 0; i < nelem; i++) {
 165                 push(a[i]);
 166                 push(p2);
 167                 multiply();
 168                 b[i] = pop();
 169         }
 170
 171         push(p3);
 172         restore();
 173 }
 174
 175 void
 176 scalar_times_tensor(void)
 177 {
 178         int i, ndim, nelem;
 179         U **a, **b;
 180
 181         save();
 182
 183         p2 = pop();
 184         p1 = pop();
 185
 186         ndim = p2->u.tensor->ndim;
 187         nelem = p2->u.tensor->nelem;
 188
 189         p3 = alloc_tensor(nelem);
 190
 191         p3->u.tensor->ndim = ndim;
 192
 193         for (i = 0; i < ndim; i++)
 194                 p3->u.tensor->dim[i] = p2->u.tensor->dim[i];
 195
 196         a = p2->u.tensor->elem;
 197         b = p3->u.tensor->elem;
 198
 199         for (i = 0; i < nelem; i++) {
 200                 push(p1);
 201                 push(a[i]);
 202                 multiply();
 203                 b[i] = pop();
 204         }
 205
 206         push(p3);
 207
 208         restore();
 209 }
 210
 211 int
 212 is_square_matrix(U *p)
 213 {
 214         if (istensor(p) && p->u.tensor->ndim == 2 && p->u.tensor->dim[0] == p->u.tensor->dim[1])
 215                 return 1;
 216         else
 217                 return 0;
 218 }
 219
 220 //-----------------------------------------------------------------------------
 221 //
 222 //      gradient of tensor
 223 //
 224 //-----------------------------------------------------------------------------
 225
 226 void
 227 d_tensor_tensor(void)
 228 {
 229         int i, j, ndim, nelem;
 230         U **a, **b, **c;
 231
 232         ndim = p1->u.tensor->ndim;
 233         nelem = p1->u.tensor->nelem;
 234
 235         if (ndim + 1 >= MAXDIM)
 236                 goto dont_evaluate;
 237
 238         p3 = alloc_tensor(nelem * p2->u.tensor->nelem);
 239
 240         p3->u.tensor->ndim = ndim + 1;
 241
 242         for (i = 0; i < ndim; i++)
 243                 p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
 244
 245         p3->u.tensor->dim[ndim] = p2->u.tensor->dim[0];
 246
 247         a = p1->u.tensor->elem;
 248         b = p2->u.tensor->elem;
 249         c = p3->u.tensor->elem;
 250
 251         for (i = 0; i < nelem; i++) {
 252                 for (j = 0; j < p2->u.tensor->nelem; j++) {
 253                         push(a[i]);
 254                         push(b[j]);
 255                         derivative();
 256                         c[i * p2->u.tensor->nelem + j] = pop();
 257                 }
 258         }
 259
 260         push(p3);
 261
 262         return;
 263
 264 dont_evaluate:
 265
 266         push_symbol(DERIVATIVE);
 267         push(p1);
 268         push(p2);
 269         list(3);
 270 }
 271
 272 //-----------------------------------------------------------------------------
 273 //
 274 //      gradient of scalar
 275 //
 276 //-----------------------------------------------------------------------------
 277
 278 void
 279 d_scalar_tensor(void)
 280 {
 281         int i;
 282         U **a, **b;
 283
 284         p3 = alloc_tensor(p2->u.tensor->nelem);
 285
 286         p3->u.tensor->ndim = 1;
 287
 288         p3->u.tensor->dim[0] = p2->u.tensor->dim[0];
 289
 290         a = p2->u.tensor->elem;
 291         b = p3->u.tensor->elem;
 292
 293         for (i = 0; i < p2->u.tensor->nelem; i++) {
 294                 push(p1);
 295                 push(a[i]);
 296                 derivative();
 297                 b[i] = pop();
 298         }
 299
 300         push(p3);
 301 }
 302
 303 //-----------------------------------------------------------------------------
 304 //
 305 //      Derivative of tensor
 306 //
 307 //-----------------------------------------------------------------------------
 308
 309 void
 310 d_tensor_scalar(void)
 311 {
 312         int i;
 313         U **a, **b;
 314
 315         p3 = alloc_tensor(p1->u.tensor->nelem);
 316
 317         p3->u.tensor->ndim = p1->u.tensor->ndim;
 318
 319         for (i = 0; i < p1->u.tensor->ndim; i++)
 320                 p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
 321
 322         a = p1->u.tensor->elem;
 323         b = p3->u.tensor->elem;
 324
 325         for (i = 0; i < p1->u.tensor->nelem; i++) {
 326                 push(a[i]);
 327                 push(p2);
 328                 derivative();
 329                 b[i] = pop();
 330         }
 331
 332         push(p3);
 333 }
 334
 335 int
 336 compare_tensors(U *p1, U *p2)
 337 {
 338         int i;
 339
 340         if (p1->u.tensor->ndim < p2->u.tensor->ndim)
 341                 return -1;
 342
 343         if (p1->u.tensor->ndim > p2->u.tensor->ndim)
 344                 return 1;
 345
 346         for (i = 0; i < p1->u.tensor->ndim; i++) {
 347                 if (p1->u.tensor->dim[i] < p2->u.tensor->dim[i])
 348                         return -1;
 349                 if (p1->u.tensor->dim[i] > p2->u.tensor->dim[i])
 350                         return 1;
 351         }
 352
 353         for (i = 0; i < p1->u.tensor->nelem; i++) {
 354                 if (equal(p1->u.tensor->elem[i], p2->u.tensor->elem[i]))
 355                         continue;
 356                 if (lessp(p1->u.tensor->elem[i], p2->u.tensor->elem[i]))
 357                         return -1;
 358                 else
 359                         return 1;
 360         }
 361
 362         return 0;
 363 }
 364
 365 //-----------------------------------------------------------------------------
 366 //
 367 //      Raise a tensor to a power
 368 //
 369 //      Input:          p1      tensor
 370 //
 371 //                      p2      exponent
 372 //
 373 //      Output:         Result on stack
 374 //
 375 //-----------------------------------------------------------------------------
 376
 377 void
 378 power_tensor(void)
 379 {
 380         int i, k, n;
 381
 382         // first and last dims must be equal
 383
 384         k = p1->u.tensor->ndim - 1;
 385
 386         if (p1->u.tensor->dim[0] != p1->u.tensor->dim[k]) {
 387                 push_symbol(POWER);
 388                 push(p1);
 389                 push(p2);
 390                 list(3);
 391                 return;
 392         }
 393
 394         push(p2);
 395
 396         n = pop_integer();
 397
 398         if (n == (int) 0x80000000) {
 399                 push_symbol(POWER);
 400                 push(p1);
 401                 push(p2);
 402                 list(3);
 403                 return;
 404         }
 405
 406         if (n == 0) {
 407                 if (p1->u.tensor->ndim != 2)
 408                         stop("power(tensor,0) with tensor rank not equal to 2");
 409                 n = p1->u.tensor->dim[0];
 410                 p1 = alloc_tensor(n * n);
 411                 p1->u.tensor->ndim = 2;
 412                 p1->u.tensor->dim[0] = n;
 413                 p1->u.tensor->dim[1] = n;
 414                 for (i = 0; i < n; i++)
 415                         p1->u.tensor->elem[n * i + i] = one;
 416                 push(p1);
 417                 return;
 418         }
 419
 420         if (n < 0) {
 421                 n = -n;
 422                 push(p1);
 423                 inv();
 424                 p1 = pop();
 425         }
 426
 427         push(p1);
 428
 429         for (i = 1; i < n; i++) {
 430                 push(p1);
 431                 inner();
 432                 if (iszero(stack[tos - 1]))
 433                         break;
 434         }
 435 }
 436
 437 void
 438 copy_tensor(void)
 439 {
 440         int i;
 441
 442         save();
 443
 444         p1 = pop();
 445
 446         p2 = alloc_tensor(p1->u.tensor->nelem);
 447
 448         p2->u.tensor->ndim = p1->u.tensor->ndim;
 449
 450         for (i = 0; i < p1->u.tensor->ndim; i++)
 451                 p2->u.tensor->dim[i] = p1->u.tensor->dim[i];
 452
 453         for (i = 0; i < p1->u.tensor->nelem; i++)
 454                 p2->u.tensor->elem[i] = p1->u.tensor->elem[i];
 455
 456         push(p2);
 457
 458         restore();
 459 }
 460
 461 // Tensors with elements that are also tensors get promoted to a higher rank.
 462
 463 static void
 464 promote_tensor(void)
 465 {
 466         int i, j, k, nelem, ndim;
 467
 468         save();
 469
 470         p1 = pop();
 471
 472         if (!istensor(p1)) {
 473                 push(p1);
 474                 restore();
 475                 return;
 476         }
 477
 478         p2 = p1->u.tensor->elem[0];
 479
 480         for (i = 1; i < p1->u.tensor->nelem; i++)
 481                 if (!compatible(p2, p1->u.tensor->elem[i]))
 482                         stop("Cannot promote tensor due to inconsistent tensor components.");
 483
 484         if (!istensor(p2)) {
 485                 push(p1);
 486                 restore();
 487                 return;
 488         }
 489
 490         ndim = p1->u.tensor->ndim + p2->u.tensor->ndim;
 491
 492         if (ndim > MAXDIM)
 493                 stop("tensor rank > 24");
 494
 495         nelem = p1->u.tensor->nelem * p2->u.tensor->nelem;
 496
 497         p3 = alloc_tensor(nelem);
 498
 499         p3->u.tensor->ndim = ndim;
 500
 501         for (i = 0; i < p1->u.tensor->ndim; i++)
 502                 p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
 503
 504         for (j = 0; j < p2->u.tensor->ndim; j++)
 505                 p3->u.tensor->dim[i + j] = p2->u.tensor->dim[j];
 506
 507         k = 0;
 508
 509         for (i = 0; i < p1->u.tensor->nelem; i++) {
 510                 p2 = p1->u.tensor->elem[i];
 511                 for (j = 0; j < p2->u.tensor->nelem; j++)
 512                         p3->u.tensor->elem[k++] = p2->u.tensor->elem[j];
 513         }
 514
 515         push(p3);
 516
 517         restore();
 518 }
 519
 520 static int
 521 compatible(U *p, U *q)
 522 {
 523         int i;
 524
 525         if (!istensor(p) && !istensor(q))
 526                 return 1;
 527
 528         if (!istensor(p) || !istensor(q))
 529                 return 0;
 530
 531         if (p->u.tensor->ndim != q->u.tensor->ndim)
 532                 return 0;
 533
 534         for (i = 0; i < p->u.tensor->ndim; i++)
 535                 if (p->u.tensor->dim[i] != q->u.tensor->dim[i])
 536                         return 0;
 537
 538         return 1;
 539 }
 540
 541 #if SELFTEST
 542
 543 static char *s[] = {
 544
 545         "#test_tensor",
 546
 547         "a=(1,2,3)",
 548         "",
 549
 550         "b=(4,5,6)",
 551         "",
 552
 553         "c=(7,8,9)",
 554         "",
 555
 556         "rank((a,b,c))",
 557         "2",
 558
 559         "(a,b,c)",
 560         "((1,2,3),(4,5,6),(7,8,9))",
 561
 562         // check tensor promotion
 563
 564         "((1,0),(0,0))",
 565         "((1,0),(0,0))",
 566
 567         "a=quote(a)",
 568         "",
 569
 570         "b=quote(b)",
 571         "",
 572
 573         "c=quote(c)",
 574         "",
 575 };
 576
 577 void
 578 test_tensor(void)
 579 {
 580         test(__FILE__, s, sizeof s / sizeof (char *));
 581 }
 582
 583 #endif