117 lines
1.6 KiB
C++
117 lines
1.6 KiB
C++
// Outer product of tensors
|
|
|
|
#include "stdafx.h"
|
|
#include "defs.h"
|
|
|
|
void
|
|
eval_outer(void)
|
|
{
|
|
p1 = cdr(p1);
|
|
push(car(p1));
|
|
eval();
|
|
p1 = cdr(p1);
|
|
while (iscons(p1)) {
|
|
push(car(p1));
|
|
eval();
|
|
outer();
|
|
p1 = cdr(p1);
|
|
}
|
|
}
|
|
|
|
void
|
|
outer(void)
|
|
{
|
|
save();
|
|
p2 = pop();
|
|
p1 = pop();
|
|
if (istensor(p1) && istensor(p2))
|
|
yyouter();
|
|
else {
|
|
push(p1);
|
|
push(p2);
|
|
if (istensor(p1))
|
|
tensor_times_scalar();
|
|
else if (istensor(p2))
|
|
scalar_times_tensor();
|
|
else
|
|
multiply();
|
|
}
|
|
restore();
|
|
}
|
|
|
|
void
|
|
yyouter(void)
|
|
{
|
|
int i, j, k, ndim, nelem;
|
|
|
|
ndim = p1->u.tensor->ndim + p2->u.tensor->ndim;
|
|
|
|
if (ndim > MAXDIM)
|
|
stop("outer: rank of result exceeds maximum");
|
|
|
|
nelem = p1->u.tensor->nelem * p2->u.tensor->nelem;
|
|
|
|
p3 = alloc_tensor(nelem);
|
|
|
|
p3->u.tensor->ndim = ndim;
|
|
|
|
for (i = 0; i < p1->u.tensor->ndim; i++)
|
|
p3->u.tensor->dim[i] = p1->u.tensor->dim[i];
|
|
|
|
j = i;
|
|
|
|
for (i = 0; i < p2->u.tensor->ndim; i++)
|
|
p3->u.tensor->dim[j + i] = p2->u.tensor->dim[i];
|
|
|
|
k = 0;
|
|
|
|
for (i = 0; i < p1->u.tensor->nelem; i++)
|
|
for (j = 0; j < p2->u.tensor->nelem; j++) {
|
|
push(p1->u.tensor->elem[i]);
|
|
push(p2->u.tensor->elem[j]);
|
|
multiply();
|
|
p3->u.tensor->elem[k++] = pop();
|
|
}
|
|
|
|
push(p3);
|
|
}
|
|
|
|
static char *s[] = {
|
|
|
|
"outer(a,b)",
|
|
"a*b",
|
|
|
|
"outer(a,(b1,b2))",
|
|
"(a*b1,a*b2)",
|
|
|
|
"outer((a1,a2),b)",
|
|
"(a1*b,a2*b)",
|
|
|
|
"H33=hilbert(3)",
|
|
"",
|
|
|
|
"H44=hilbert(4)",
|
|
"",
|
|
|
|
"H55=hilbert(5)",
|
|
"",
|
|
|
|
"H3344=outer(H33,H44)",
|
|
"",
|
|
|
|
"H4455=outer(H44,H55)",
|
|
"",
|
|
|
|
"H33444455=outer(H33,H44,H44,H55)",
|
|
"",
|
|
|
|
"simplify(inner(H3344,H4455)-contract(H33444455,4,5))",
|
|
"0",
|
|
};
|
|
|
|
void
|
|
test_outer(void)
|
|
{
|
|
test(__FILE__, s, sizeof s / sizeof (char *));
|
|
}
|