Re: memory-efficient circshift mex?

"Matt " <mjacobson.removethis@xxxxxxxxxxxxx> wrote in message <gbjirr$i5o$1@xxxxxxxxxxxxxxxxxx>...
I'm encountering an annoying difficulty with circhsift.m, and I'm wondering if anyone possesses an alternative mex implementation.

I have a large 3D array X (e.g. 700x700x350) and I want to circshift it in the 3rd dimension by a small number of slices, for example

circshift(X,[0 0 1])

The memory-efficient thing to do would be to

(1) Separately store only the slices that must be circulated to the bottom of the array, in this case
Y=X(:,:,end), then

(2) Loop backward through all the slices, overwriting
X(:,:,i) with X(:,:,i-1), then

(3) Overwrite X(:,:,1) with Y

In general, this only requires the pre-allocation of N slices where N is the amount of the shift.

However, the native circhsift.m naively allocates an additional array of size(X) regardless of N. This is of course difficult for me given the size of X. The fact that it is implemented as an m-file also doesn't help.

Does anyone have something that implements it the better way? I checked the FEX already to no avail.


I was inspired, so I wrote a mex file for you. It does circular shifting of a 3-dimensional or higher double array based on the last dimension. Does the shifting in-place, so only works if the input array is not shared. If you think this is worthwhile let me know and I will fix it up with more help & an associated m-file and post it to the FEX. Enjoy.

James Tursa

------------------------------------------------------

mex instructions if you don't know:

Name the file circshiftlast.c and place it somewhere on the MATLAB path.
Set your default directory to that directory.
Type the following commands:

mex -setup
(Enter)
(Select a C compiler, such as the built in lcc)
(Enter Enter etc.)
mex circshiftlast.c

Now you have a function circshiftlast that can be called.

-------------------------------------------------------

// circshiftlast, circular shift the last dimension in place
// Programmer: James Tursa
// Circular shift a 3-dimensional or higher matrix in the last dimension.
// 1st input argument must be a 3-dimensional or higher double array.
// 2nd input argument must be a scalar with integral value.
// negative = shift up, positive = shift down
// Does shifting in-place, so input argument cannot be shared.

#include <string.h>
#include "mex.h"

#ifndef mwSize
#define mwSize int
#endif

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwSize ndim, numel, k, slicenumel, slicebytes, kshift;
mwSize *dims;
double d;
double *dpr, *dpi, *dptemp;

if( nrhs != 2 )
mexErrMsgTxt("Must have 2 input arguments");
if( nlhs != 0 )
mexErrMsgTxt("Must have 0 output arguments");
if( !mxIsDouble(prhs[0]) )
mexErrMsgTxt("1st input argument must be double");
ndim = mxGetNumberOfDimensions(prhs[0]);
if( ndim < 3 )
mexErrMsgTxt("1st input argument must be 3-dimensional or higher");
if( !mxIsNumeric(prhs[1]) || mxGetNumberOfElements(prhs[1]) != 1 )
mexErrMsgTxt("2nd input argument must be numeric scalar");
d = mxGetScalar(prhs[1]);
if( mxIsInf(d) || mxIsNaN(d) )
mexErrMsgTxt("2nd input argument cannot be inf or NaN");
kshift = (int) d;
if( kshift != d )
mexErrMsgTxt("2nd input argument must be integral value");
dpr = mxGetPr(prhs[0]);
dpi = mxGetPi(prhs[0]);
dims = mxGetDimensions(prhs[0]);
numel = mxGetNumberOfElements(prhs[0]);
if( numel == 0 )
return;
k = dims[ndim-1]; // size of last dimension
slicenumel = numel / k; // number of elements in one slice
slicebytes = slicenumel * sizeof(double); // number of bytes in one slice
kshift = kshift % k; // number of slices to shift
if( kshift > k/2 ) { // modify kshift to minimize amount of temp memory copied
kshift -= k;
} else if( kshift < -k/2 ) {
kshift += k;
} else if( kshift == 0 ) {
return;
}
if( kshift < 0 ) { // shift up
dptemp = mxMalloc(-kshift*slicebytes);
memcpy(dptemp, dpr, -kshift*slicebytes);
memmove(dpr, dpr-kshift*slicenumel, (k+kshift)*slicebytes);
memcpy(dpr+(k+kshift)*slicenumel, dptemp, -kshift*slicebytes);
if( dpi != NULL ) {
memcpy(dptemp, dpi, -kshift*slicebytes);
memmove(dpi, dpi-kshift*slicenumel, (k+kshift)*slicebytes);
memcpy(dpi+(k+kshift)*slicenumel, dptemp, -kshift*slicebytes);
}
mxFree(dptemp);
} else { // shift down
dptemp = mxMalloc(kshift*slicebytes);
memcpy(dptemp, dpr+(k-kshift)*slicenumel, kshift*slicebytes);
memmove(dpr+kshift*slicenumel, dpr, (k-kshift)*slicebytes);
memcpy(dpr, dptemp, kshift*slicebytes);
if( dpi != NULL ) {
memcpy(dptemp, dpi+(k-kshift)*slicenumel, kshift*slicebytes);
memmove(dpi+kshift*slicenumel, dpi, (k-kshift)*slicebytes);
memcpy(dpi, dptemp, kshift*slicebytes);
}
mxFree(dptemp);
}
}


.

Relevant Pages

  • Re: memory-efficient circshift mex?
    ... Separately store only the slices that must be circulated to the bottom of the array, ... the native circhsift.m naively allocates an additional array of sizeregardless of N. This is of course difficult for me given the size of X. ... // Circular shift a 3-dimensional or higher matrix in the last dimension. ... mwSize ndim, numel, k, slicenumel, slicebytes, kshift; ...
    (comp.soft-sys.matlab)
  • Re: Numerical Python question
    ... you're now looping over ten slices of tens of thousands. ... >> original array. ... loop in Python. ...
    (comp.lang.python)
  • Re: Newbie Looking for Advice/Comments on Script
    ... > array, then write out the array in a specific order. ... Please check the Perl FAQ *before* posting: ... Google Groups to do so. ... or, preferred, don't interpolate the slices. ...
    (comp.lang.perl.misc)
  • Re: Buffering TCPIP data
    ... It would be very helpful if you could outline the benefits of a circular ... You'd implement the circular buffer with an underlying byte array - it ... with no extra copying involved - i.e. you don't copy the ...
    (microsoft.public.dotnet.languages.csharp)
  • Re: Surprise in array concatenation
    ... > what about a recursive subprogram? ... carries with it the original array. ... > actually I have found Ada's passing of slice bounds quite elegant ... > Ada's slices, i.e. strings include slices which are strings). ...
    (comp.lang.ada)