35605 # MATLAB: MEX matrix division gives different result than m-file

I've used MATLAB's coder tool to create a MEX version of the matrix exponential function, to be used in another set of functions. The issue is, the MEX version gives different results than the original m-file.

After debugging it, I believe that the reason this is, is because the MEX file and the m-file do not do matrix division (\) the same. Or the MEX file has issues with it in the first place. All the variables leading up to the line where the matrix division occurs are equivalent on both sides.

This is the line where the issue occurs:

```F = (V-U)\(2*U) + I ```

Where I is the identity matrix of the size of V and U.

What is the reason behind the discrepancy when a MEX file does matrix division, and how can I fix this issue? Can this line of code be re-written without the division?

I have no problem generating C code from such an operation.

Here is a test function I tried:

### myfcn.m

```function F = myfcn(U,V) I = eye(size(U)); F = (V-U)\(2*U) + I; end ```

here's a test script we'll use to verify the results:

### test_myfcn.m

```U = rand(5); V = rand(5); F = myfcn(U,V); ```

I start by launching the code generation tool (`ccoder`), create a new project set to produce a MEX-file, then add the `myfnc.m` function from before as entry-point. Then I define both input variables types as:

```double (:Inf x :Inf) ```

which specifies an MxN matrix of unbounded size of type double.

Finally we can build the project. This produces `myfcn_mex.mexw64`.

Testing both the original M-function and the generated MEX-function, I get pretty much the same result (difference is in the order of machine epsilon):

```>> F = myfnc(U,V); >> FF = myfcn_mex(U,V); >> norm(F-FF) ans = 1.4391e-14 ```

MATLAB made a minor change to the EXPM algorithm in the 2014a release. MATLAB Coder implementations are separate, and the corresponding change has not yet been made in the algorithm for code generation, so some discrepancies are possible there. The new implementation does `(V - U)\(2*U) + I`, whereas the old one does `(V - U)\(V + U)`. These are mathematically equivalent but will give different round-off behaviors in general.