* mex.oct is a shared library.  We should rename it libmex.so and link to
it instead of including a separate copy of mex.o in each oct file.  Only
thing is, how do we install it beside version of liboctave we are using,
and how do we make the loader find the correct version if we have multiple
versions of octave installed.

* More complete testing
* Better error handling in mex
* Use tempory files for Octave glue?
* mex should translate -V4 to -DV4
* mex should remove mex_*.cc when done
* clean up warnings in myset.c

* If you want to support Octave and Matlab with the same code base
without penalizing either, you will need to write a high level layer
which supports both.  

Particularly, you must prohibit direct modification of the fields of 
the mxArray type, which means eliminating mxSet{Pi,Pr,N,M} from your
matlab code.  You can get the same functionality by adding the routines:
	mxResize (to add/remove rows and columns), 
	mxReshape (to change the shape without changing the data),
	mxMakeComplex (to change a double matrix into a complex matrix),
	mxMakeReal (to change a complex matrix into a double matrix),
	mxAssign (to replace an array with a new one)
These are all easy enough to define in Matlab. The octave routines
will represent the data using either ComplexMatrix or Matrix so that
an octave value can easily be constructed when it is needed.  The
matlab side can continue to use mxArray.

The function mxGetPr is problematic for values passed from octave. 
Because you cannot get a reference to the Matrix in the ocave_value,
you must assign it.  Since mxGetPr returns a pointer to the contents
of the array, you must use tmp.fortran_vec().  Since tmp is assigned,
and arrays are copy-on-modify, the simple act of mxGetPr will trigger
a copy.  A solution is to trick the compiler: return the const pointer
given by tmp.data(), but declare the function header included in the
mex files as returning a nonconst pointer.

You will need to hide the difference in storage between complex
arrays, which matlab defines as two matrices (one for real and one for
imaginary), and Octave defines as one matrix (with real and imaginary
values alternating).  This is easy to handle by returning Pi as Pr+1,
and multiplying the current matlab index by 2. In octave this would be:
  #define mxStride 2
  double *mxGetPr(mxArray *m) { return m->data(); }
  double *mxGetPzr(mxArray *m) { return (double *)(cm->data()); }
  double *mxGetPzi(mxArray *m) { return (double *)(cm->data()) + 1; }
and for matlab
  #define mxStride 1
  #define mxGetPzr mxGetPr
  #define mxGetPzi mxGetPi
Then you could implement abs(z) as something like:
  int M, N;
  double *pr, *pi, *pd;
  M = mxGetM(prhs[0]);
  N = mxGetN(prhs[0]);
  pr = mxGetPzr(prhs[0]);
  pi = mxGetPzi(prhs[0]);
  plhs[0] = mxCreateDoubleMatrix(M,N);
  pd = mxGetPr(plhs[0]);
  for (j=0; j < M*N; j++)
    pd[j] = sqrt(pr[j*mxStride]*pr[j*mxStride]+pi[j*mxStride]*pi[j*mxStrde]);

If you have any functions which take real and imaginary parts as
separate arrays, these should be rewritten to accept an mxArray
applying the solution above, otherwise you will have to do something
ugly like:

#ifdef HAVE_OCTAVE
  double *pr, *pi, *newr, *newi;
  pr = mxGetPzr(m);
  pi = mxGetPzi(m);
  newr = mxMalloc(M*N*sizeof(double));
  newi = mxMalloc(M*N*sizeof(double));
  for (j=0; j < M*N*mxStride; j+=mxStride) newr[j]=pr[j], newi[j]=pi[j];
  complex_fun(newr,newi);
  for (j=0; j < M*N*mxStride; j+=mxStride) pr[j]=newr[j], pi[j]=newi[j];
  mxFree(newr);
  mxFree(newi);
#else
   complex_fun(mxGetPr(m),mxGetPi(m));
#endif