0

I am trying to speed up some of my methods and can't figure out, why two of my methods calculate correctly but fail to run several times.

The Setting: I am getting a (possibly huge) set of vectors (columns) in p and v, let's assume for simplicity just 3D unit vectors. For each pair of columns I would like to compute the corresponding q = p*cos(||v||) + v/||v|| * sin(||v||). For convenience (and my larger example) I would like to use Eigen and its MatrixXd class.

#include "mex.h"
#include "math.h"
#include <stdlib.h>
#include <Eigen/Dense>
using namespace std;
using namespace Eigen;

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[])   {
    const mxArray *IpA, *IvA; //Just names
    double *Ip, *Iv, *Oq;
    mwSize i,j;
    size_t ItemSize,N;
    // Inputs
    IpA = prhs[0];
    IvA = prhs[1];
    Ip = mxGetPr(IpA); // Base Points
    Iv = mxGetPr(IvA); // Directions

    ItemSize = mxGetM(IpA);
    N = mxGetN(IpA);
    if ( (ItemSize != mxGetM(IvA)) || (N != mxGetN(IvA)) )
        mexErrMsgTxt("p and v have to be of same size.");

    plhs[0] = mxCreateDoubleMatrix((mwSize)ItemSize,(mwSize)N,mxREAL);
    Oq = mxGetPr(plhs[0]);
    // Norm directional vectors V
    MatrixXd normV(ItemSize,1);
    for (i=0; i<N; i++) {
        normV(i) = 0;
        for (j=0; j<ItemSize; j++) {
            normV(i) += Iv[j + ItemSize*i]*Iv[j + ItemSize*i];
        }
    }
    //Compute result
    normV = normV.array().sqrt();
    for (i=0; i<N; i++) {
        for (j=0; j<ItemSize; j++) {
            if (normV(i)==0)
                Oq[j + ItemSize*i] = 0;
            else {
                Oq[j + ItemSize*i] = Ip[j + ItemSize*i]*cos(normV(i)) + Iv[j + ItemSize*i]*sin(normV(i))/normV(i);
            }
        }
    }
}

I'm following the tutorial by Pascal Getreuer.

The Problem: Nevertheless, if I am calling this compiled function (Let's name it SnExp I'm getting the correct result for single vectors or a few ones. If I am using a large number (3x8000which is not that large) or call the functions a few times (let's take SnExp(repmat([1;0;0],1,1000),repmat([1;0;0],1,1000)) which should return a matrix containing 1000 columns [1.3818;0;0].

Matlab just freezes or crashes and I can't figure out why. Can you confirm that crash and help me fix it? That would be awesome.

Ronny
  • 139
  • 9
  • You should opt to debug your code instead. Check out these instructions on how to enable debugging of your MEX code so you can step through it yourself and determine where and why it's crashing: http://stackoverflow.com/questions/23714141/preventing-a-mex-file-from-crashing-in-matlab/23714301#23714301 – rayryeng Apr 12 '15 at 02:07
  • Thanks for the link, you're right; when you're used to the Matlab debugger it's not easy to see the problem in C++. I'll check you're instructions. Running Mac OS here btw. – Ronny Apr 12 '15 at 04:38
  • Totally understand. MATLAB debugger is quite good. Debugging MEX files are a pain, and I'm sorry that you can't do it in MATLAB... if you're running in Mac OS, you'll have to source out the debugging to XCode and debug it there. That's really the only way you'll see what's going on. The link I posted has a Mac OS setup. Good luck! – rayryeng Apr 12 '15 at 05:03
  • 1
    Already running debugger now, thanks for the link, that was really helpful. Well - the functions already run in Matlab but they are too slow. That's at least why i know it's computing correctly - i have a reference implementation :) – Ronny Apr 12 '15 at 05:30
  • You seem to declare normV a size of ItemSizex1. But you seem to index it from 0 to N-1 in the first for loop. Is that correct? – Navan Apr 13 '15 at 13:37
  • No, that's actually the mistake made :) Noticed that yesterday, too thanks to the debugger. And that was the only mistake now it runs well. Thanks for the hint to the debugger and your hint to the mistake. – Ronny Apr 13 '15 at 13:58

0 Answers0