outer.cpp

   1 // Outer product of tensors
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 void
   7 eval_outer(void)
   8 {
   9         p1 = cdr(p1);
  10         push(car(p1));
  11         eval();
  12         p1 = cdr(p1);
  13         while (iscons(p1)) {
  14                 push(car(p1));
  15                 eval();
  16                 outer();
  17                 p1 = cdr(p1);
  18         }
  19 }
  20
  21 void
  22 outer(void)
  23 {
  24         save();
  25         p2 = pop();
  26         p1 = pop();
  27         if (istensor(p1) && istensor(p2))
  28                 yyouter();
  29         else {
  30                 push(p1);
  31                 push(p2);
  32                 if (istensor(p1))
  33                         tensor_times_scalar();
  34                 else if (istensor(p2))
  35                         scalar_times_tensor();
  36                 else
  37                         multiply();
  38         }
  39         restore();
  40 }
  41
  42 void
  43 yyouter(void)
  44 {
  45         int i, j, k, ndim, nelem;
  46
  47         ndim = p1->u.tensor->ndim + p2->u.tensor->ndim;
  48
  49         if (ndim > MAXDIM)
  50                 stop("outer: rank of result exceeds maximum");
  51
  52         nelem = p1->u.tensor->nelem * p2->u.tensor->nelem;
  53
  54         p3 = alloc_tensor(nelem);
  55
  56         p3->u.tensor->ndim = ndim;
  57
  58         for (i = 0; i < p1->u.tensor->ndim; i++)
  59                 p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
  60
  61         j = i;
  62
  63         for (i = 0; i < p2->u.tensor->ndim; i++)
  64                 p3->u.tensor->dim[j + i] = p2->u.tensor->dim[i];
  65
  66         k = 0;
  67
  68         for (i = 0; i < p1->u.tensor->nelem; i++)
  69                 for (j = 0; j < p2->u.tensor->nelem; j++) {
  70                         push(p1->u.tensor->elem[i]);
  71                         push(p2->u.tensor->elem[j]);
  72                         multiply();
  73                         p3->u.tensor->elem[k++] = pop();
  74                 }
  75
  76         push(p3);
  77 }
  78
  79 #if SELFTEST
  80
  81 static char *s[] = {
  82
  83         "outer(a,b)",
  84         "a*b",
  85
  86         "outer(a,(b1,b2))",
  87         "(a*b1,a*b2)",
  88
  89         "outer((a1,a2),b)",
  90         "(a1*b,a2*b)",
  91
  92         "H33=hilbert(3)",
  93         "",
  94
  95         "H44=hilbert(4)",
  96         "",
  97
  98         "H55=hilbert(5)",
  99         "",
 100
 101         "H3344=outer(H33,H44)",
 102         "",
 103
 104         "H4455=outer(H44,H55)",
 105         "",
 106
 107         "H33444455=outer(H33,H44,H44,H55)",
 108         "",
 109
 110         "simplify(inner(H3344,H4455)-contract(H33444455,4,5))",
 111         "0",
 112 };
 113
 114 void
 115 test_outer(void)
 116 {
 117         test(__FILE__, s, sizeof s / sizeof (char *));
 118 }
 119
 120 #endif