2

I have struggling with an implementation of MEX file, I get run time error! I would like to given two inputs : inMatToInv, inMatToMul I would like to compute : inv(inMatToInv)*inMatToMul, I assume that the user enters only 3x3 matrices. I also assume that the inMatToInv in invertible.

    /*==========================================================
 * inv_and_mul_3by3.c 
 * inverse 3x3 matrix and multiply it by another 3x3 matrix
 * The calling syntax is: outMatrix = (mat_to_inv,mat_to_mul)
 *========================================================*/

#include "mex.h"

/* The computational routine */
void inv_and_mul_3by3(double *mat_to_inv, double *mat_to_mul, double *out)
{
    // Description : out = inv(mat_to_inv)*mat_to_mul
    double *inversed;
    double det;

    /* Calculate the matrix deteminant */
    det = mat_to_inv[0]*(mat_to_inv[4]*mat_to_inv[8]-mat_to_inv[7]*mat_to_inv[5]) - mat_to_inv[3]*(mat_to_inv[1]*mat_to_inv[8]-mat_to_inv[7]*mat_to_inv[2])+mat_to_inv[6]*(mat_to_inv[1]*mat_to_inv[5]-mat_to_inv[4]*mat_to_inv[2]);

    // Calcualte the inversed matrix
    inversed[0] = (mat_to_inv[4]*mat_to_inv[8] - mat_to_inv[7]*mat_to_inv[5])/det;
    inversed[3] = (mat_to_inv[6]*mat_to_inv[5] - mat_to_inv[3]*mat_to_inv[8])/det;
    inversed[6] = (mat_to_inv[3]*mat_to_inv[7] - mat_to_inv[6]*mat_to_inv[4])/det;
    inversed[1] = (mat_to_inv[7]*mat_to_inv[2] - mat_to_inv[1]*mat_to_inv[8])/det;
    inversed[4] = (mat_to_inv[0]*mat_to_inv[8] - mat_to_inv[6]*mat_to_inv[2])/det;
    inversed[7] = (mat_to_inv[6]*mat_to_inv[1] - mat_to_inv[0]*mat_to_inv[7])/det;
    inversed[2] = (mat_to_inv[1]*mat_to_inv[5] - mat_to_inv[4]*mat_to_inv[2])/det;
    inversed[5] = (mat_to_inv[3]*mat_to_inv[2] - mat_to_inv[0]*mat_to_inv[5])/det;
    inversed[8] = (mat_to_inv[0]*mat_to_inv[4] - mat_to_inv[3]*mat_to_inv[1])/det;

    out[0] = inversed[0]*mat_to_mul[0] + inversed[3]*mat_to_mul[1] + inversed[6]*mat_to_mul[2];
    out[1] = inversed[1]*mat_to_mul[0] + inversed[4]*mat_to_mul[1] + inversed[7]*mat_to_mul[2];
    out[2] = inversed[2]*mat_to_mul[0] + inversed[5]*mat_to_mul[1] + inversed[8]*mat_to_mul[2];
    out[3] = inversed[0]*mat_to_mul[3] + inversed[3]*mat_to_mul[4] + inversed[6]*mat_to_mul[5];
    out[4] = inversed[1]*mat_to_mul[3] + inversed[4]*mat_to_mul[4] + inversed[7]*mat_to_mul[5];
    out[5] = inversed[2]*mat_to_mul[3] + inversed[5]*mat_to_mul[4] + inversed[8]*mat_to_mul[5];
    out[6] = inversed[0]*mat_to_mul[6] + inversed[3]*mat_to_mul[7] + inversed[6]*mat_to_mul[8];
    out[7] = inversed[1]*mat_to_mul[6] + inversed[4]*mat_to_mul[7] + inversed[7]*mat_to_mul[8];
    out[8] = inversed[2]*mat_to_mul[6] + inversed[5]*mat_to_mul[7] + inversed[8]*mat_to_mul[8];
}

/* The gateway function */
void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])
{
    double *inMatToInv;             /* 3x3 input matrix that is inversed */
    double *inMatToMul;             /* 3x3 input matrix that multiply the inversed matrix*/ 
    double *outMatrix;              /* 3x3 output matrix */

    /* create pointers to the real data in the input matrix  */
    inMatToInv = mxGetPr(prhs[0]);
    inMatToMul = mxGetPr(prhs[1]);

    /* create the output matrix */
    plhs[0] = mxCreateDoubleMatrix(3,3,mxREAL);

    /* get a pointer to the real data in the output matrix */
    outMatrix = mxGetPr(plhs[0]);

    /* call the computational routine */
    inv_and_mul_3by3(inMatToInv,inMatToMul,outMatrix);
}

any help will be much appreciated. Thanks, Tamir

Daniel
  • 36,610
  • 3
  • 36
  • 69
Tamir Einy
  • 164
  • 1
  • 12

2 Answers2

0

You declare a pointer double *inversed, but you do not allocate any memory for it. So when you access inversed[0] you get an error. Try:

double inversed[9];
Itamar Katz
  • 9,544
  • 5
  • 42
  • 74
0

I'd recommend using the blas and lapack interface for matrix operations in a C mex file. Have a look at Mathworks documentation on doing this here.

You would be interested in the lapack routines dgetrf and dgetri for the matrix inversion and dgemm for matrix multiplication. The lapack interface can be really messy and difficult to grasp but there are lots of references online about them. For example, this stack overflow Q/A.

Community
  • 1
  • 1
ChisholmKyle
  • 449
  • 3
  • 10