codegen C output of MatrixSymbol writes undefined behavior

[stewpend0us]

from sympy import MatrixSymbol
from sympy.utilities.codegen import codegen

M = MatrixSymbol('M',2,1)
F = 2*M
codegen( ('example', F), "C", "bad_function", to_files=True )

generates

double example(double *M) {

   double example_result;
   example_result[0] = 2*M[0];
   example_result[1] = 2*M[1];
   return example_result;

}

which is not ok in c as far as I know...I would expect something more like

void example(double * example_result, double *M) {

   example_result[0] = 2*M[0];
   example_result[1] = 2*M[1];

}

[asmeurer]

@bjodah @moorepants

[moorepants]

I don't think that codegen works with MatrixSymbols yet in this way. You'd need to do something like:

from sympy import MatrixSymbol, Matrix
from sympy.utilities.codegen import codegen

M = MatrixSymbol('M', 2, 1)
F = 2 * M
G = Matrix([F[0, 0], F[1, 0]])
stuff = codegen( ('example', F), "C", "bad_function")
print(stuff[0][1])

I think there is an open issue for this already.

[moorepants]

Master produces this:

from sympy import MatrixSymbol, Matrix
from sympy.utilities.codegen import codegen

M = MatrixSymbol('M', 2, 1)
F = 2 * M
G = Matrix([F[0, 0], F[1, 0]])
stuff = codegen( ('example', F), "C", "bad_function")
print(stuff[0][1])

## -- End pasted text --
/******************************************************************************
 *                    Code generated with sympy 1.1.2.dev                     *
 *                                                                            *
 *              See http://www.sympy.org/ for more information.               *
 *                                                                            *
 *                       This file is part of 'project'                       *
 ******************************************************************************/
#include "bad_function.h"
#include <math.h>

double example(double *M) {

   double example_result;
   example_result[0] = 2*M[0];
   example_result[1] = 2*M[1];
   return example_result;

}

[moorepants]

This all seems to work here: http://www.sympy.org/scipy-2017-codegen-tutorial/notebooks/08-cythonizing.html#2.-Generating-C-Code-with-SymPy's-codegen()

Not quite sure what I'm missing that is different here though.

[moorepants]

This works:

In [8]: import sympy as sm

In [9]: M = MatrixSymbol('M', 2, 1)

In [10]: R = MatrixSymbol('R', 2, 1)

In [11]: code = codegen(('evaluate', sm.Eq(R, 2 * M)), language='c')

In [12]: print(code[0][1])
/******************************************************************************
 *                    Code generated with sympy 1.1.2.dev                     *
 *                                                                            *
 *              See http://www.sympy.org/ for more information.               *
 *                                                                            *
 *                       This file is part of 'project'                       *
 ******************************************************************************/
#include "evaluate.h"
#include <math.h>

void evaluate(double *M, double *R) {

   R[0] = 2*M[0];
   R[1] = 2*M[1];

}

[moorepants]

I'm thinking that it misinterprets the MatMul object.

This works too:

In [16]: B = sm.Matrix([2 * M[0, 0], 3 * M[1, 0]])

In [17]: B
Out[17]: 
Matrix([
[2*M[0, 0]],
[3*M[1, 0]]])

In [18]: code = codegen(('evaluate', B), language='c')

In [19]: print(code[0][1])
/******************************************************************************
 *                    Code generated with sympy 1.1.2.dev                     *
 *                                                                            *
 *              See http://www.sympy.org/ for more information.               *
 *                                                                            *
 *                       This file is part of 'project'                       *
 ******************************************************************************/
#include "evaluate.h"
#include <math.h>

void evaluate(double *M, double *out_4466384987431072851) {

   out_4466384987431072851[0] = 2*M[0];
   out_4466384987431072851[1] = 3*M[1];

}

[stewpend0us]

This seems to work also:

F = Matrix(2*M)