Code generation for dynamic size arrays

[Morten Olsen Lysgaard]

I am using Sympy for code generation.

I would like to figure out the boundaries of what Sympy can do with respect to code generation for dynamically sized arrays.

Let's say I would like to generate a simple dot-product C-function.

void dot(int m, const double *A, const double *B, double *C) {
  for(int i=0; i<m; i++){
    C[i] = A[i]*B[i];
  }
}

Using `sympy.Matrix` I am able to generate a loop-unrolled version of this function for a fixed `m`, but this is not what I want. I want the general thing.

I have been trying to get this working using `IndexedBase` et. al. but no luck. This is the closest I have gotten:

import sys
import sympy as sp
import sympy.utilities.codegen
from sympy.utilities import codegen

m = sp.symbols('m', integer=True)
i = sp.Idx('i', m)
M = sp.IndexedBase('M', shape=(m))
N = sp.IndexedBase('N', shape=(m))
K = sp.IndexedBase('K', shape=(m))
dot =sp.Eq(K[i], M[i]*N[i])

rut_obj = codegen.make_routine('dot', dot, (m, M, N, K), language='C')
rcg = codegen.CCodeGen()
rcg.dump_c([rut_obj], sys.stdout, '')

But it outputs spurious loops which I can not explain, and will not compute the correct answer:

void dot(int m, double *M, double *N, double *K) {
   for (int i=0; i<m; i++){
      K[i] = 0;
   }
   for (int i=0; i<m; i++){
      for (int i=0; i<m; i++){
         K[i] = K[i] + M[i]*N[i];
      }
   }
}

Is there anyone that are able to generate a simple dot-product function like the one in the top of my email using Sympy code generation?

Any help would be greatly appreciated!

[Morten Olsen Lysgaard]

Ah, there is an inconsistency in the above email. I am saying I want dot-product, when in reality, what I want is an elementwise multiplication. The same question still applies though. Is this something Sympy can do?

[Aaron Meurer]

It looks like the code that automatically does contractions is messing
with what you are trying to do here.

What I would try doing is using the sympy.codegen classes to construct
the loop you want. You can then print this directly using
contract=False to disable the automatic contraction. For example

>>> from sympy.codegen import For, CodeBlock, Assignment
>>> print(ccode(For(i, Range(m), CodeBlock(Assignment(K[i], M[i]*N[i]))), contract=False))
for (i = 0; i < m; i += 1) {
K[i] = M[i]*N[i];
}

Unfortunately, the sympy.utilities.codegen classes currently do not
know about the sympy.codegen classes. The idea is to eventually
replace everything with the sympy.codegen ast nodes.

Aaron Meurer

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/7a5e17dd-91c2-4988-83d9-2a6e83b8b1b3n%40googlegroups.com.