Prev: 3d volume to 2d
Next: mex problem :: can't map file, errno=22
From: James Tursa on 22 Mar 2010 16:18
"Chad " <parishcm(a)ornl.gov> wrote in message <ho8eea$6vi$1(a)fred.mathworks.com>...
>
> James Tursa,
>
> Thank you for the offer! I hate to impose upon someones else's time, but if you're willing to do so, I'd be happy for your help. Real numbers, always.

Here it is. The routine is based on the sparse matrix cleaning code I wrote for my fast matrix multiply mtimesx FEX submission. I simply altered it to clean non-finite values from the array instead of cleaning 0's from the array. I will not try to explain the routine but simply ask you to accept it as a black box. Basically, I took some reasonably looking code (to someone who understands sparse matrix layout) and highly optimized it for speed ... ending up with a fairly ugly piece of code (even using a goto) but one that runs very fast. Be advised that it alters the sparse matrix in-place, so if you have any shared variables with this then it will cause some side-effects of altering those other variables. If that is the case let me know. The two main reasons for the speed are the in-place altering of the data and the fact that it traverses the array only once (unless there is a realloc,
in which case it traverses it twice). For the realloc I chose 1024 elements as the trigger, but you can change this to something else if you want (if the array drops by 1024 elements or more, then a realloc takes place to try to recover some memory). To use it simply put setnonfinite0.c on your path, make that directory your current directory, and then type the following:

mex -setup
(then select a C/C++ compiler from the menu, e.g., lcc)
mex setnonfinite0.c

Then you will have a function setnonfinite0 that you can call just like any other MATLAB function.

On a broader note, maybe a more generalized piece of code that sets sparse values meeting a certain criteria to another value (e.g., 0) would be useful. I will think about that one for my next project. (I am still working on updating mtimesx for 64-bit and unix support at the moment).

James Tursa

setnonfinite0.c
-----------------------------------------------------------------------------
// Clean a real sparse matrix of non-finite values in-place
// James Tursa
#include "mex.h"
#define REALLOCTOL 1024 // Number of elements dropped to trigger a realloc

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwIndex *ir;
mwIndex k, nzmax;
double *pr;

if( nlhs != 0 ) {
mexErrMsgTxt("Too many outputs.");
}
if( nrhs != 1 || !mxIsSparse(prhs[0]) || mxIsComplex(prhs[0]) ) {
mexErrMsgTxt("Need exactly one real sparse matrix input.");
}
pr = mxGetPr(prhs[0]);
ir = mxGetIr(prhs[0]);
nzmax = mxGetNzmax(prhs[0]);
k = spcleanfinite(prhs[0]);
if( k == 0 ) k = 1;
if( nzmax - k > REALLOCTOL ) {
mxSetPr(prhs[0], mxRealloc(pr, k * sizeof(*pr)));
mxSetIr(prhs[0], mxRealloc(ir, k * sizeof(*ir)));
mxSetNzmax(prhs[0], k);
}
}

//-------------------------------------------------------------------------------
//
// Clean a real sparse matrix of non-finite values. The code is a bit ugly
// because it was highly optimized for speed.
//
//-------------------------------------------------------------------------------

mwIndex spcleanfinite(mxArray *mx)
{
mwSize n, nrow;
mwIndex *ir, *jc, *iz;
mwIndex j, x, y, diff;
double *pr, *pz;

n = mxGetN(mx);
pz = pr = mxGetData(mx);
ir = mxGetIr(mx);
jc = mxGetJc(mx);

diff = 0;
for( y=0; y<n; y++ ) {
nrow = jc[y+1] - jc[y];
for( x=0; x<nrow; x++ ) {
if( mxIsInf(*pr) || mxIsNaN(*pr) ) {
iz = (ir += (pr - pz));
pz = pr;
j = y + 1;
goto copydata;
}
pr++;
}
}
return jc[n];

for( y=0; y<n; y++ ) {
j = y + 1;
nrow = (jc[j] - diff) - jc[y];
for( x=0; x<nrow; x++ ) {
if( !(mxIsInf(*pr) || mxIsNaN(*pr)) ) {
*pz++ = *pr;
*iz++ = *ir;
} else {
copydata:
diff++;
}
pr++;
ir++;
}
jc[j] -= diff;
}
return jc[n];
}
From: James Tursa on 22 Mar 2010 16:42
Of course, right after posting I remembered there is an isfinite function available in the API ... no need to resort to isinf and isnan separately. This version looks a bit cleaner but isn't going to run significantly faster. I had to post it anyway ...

James Tursa

setnonfinite0.c
-----------------------------------------------------------------------------
// Clean a real sparse matrix of non-finite values in-place
// James Tursa
#include "mex.h"
#define REALLOCTOL 1024 // Number of elements dropped to trigger a realloc

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
mwIndex *ir;
mwIndex k, nzmax;
double *pr;

if( nlhs != 0 ) {
mexErrMsgTxt("Too many outputs.");
}
if( nrhs != 1 || !mxIsSparse(prhs[0]) || mxIsComplex(prhs[0]) ) {
mexErrMsgTxt("Need exactly one real sparse matrix input.");
}
pr = mxGetPr(prhs[0]);
ir = mxGetIr(prhs[0]);
nzmax = mxGetNzmax(prhs[0]);
k = spcleanfinite(prhs[0]);
if( k == 0 ) k = 1;
if( nzmax - k > REALLOCTOL ) {
mxSetPr(prhs[0], mxRealloc(pr, k * sizeof(*pr)));
mxSetIr(prhs[0], mxRealloc(ir, k * sizeof(*ir)));
mxSetNzmax(prhs[0], k);
}
}

//-------------------------------------------------------------------------------
//
// Clean a real sparse matrix of non-finite values. The code is a bit ugly
// because it was highly optimized for speed.
//
//-------------------------------------------------------------------------------

mwIndex spcleanfinite(mxArray *mx)
{
mwSize n, nrow;
mwIndex *ir, *jc, *iz;
mwIndex j, x, y, diff;
double *pr, *pz;

n = mxGetN(mx);
pz = pr = mxGetData(mx);
ir = mxGetIr(mx);
jc = mxGetJc(mx);

diff = 0;
for( y=0; y<n; y++ ) {
nrow = jc[y+1] - jc[y];
for( x=0; x<nrow; x++ ) {
if( !mxIsFinite(*pr) ) {
iz = (ir += (pr - pz));
pz = pr;
j = y + 1;
goto copydata;
}
pr++;
}
}
return jc[n];

for( y=0; y<n; y++ ) {
j = y + 1;
nrow = (jc[j] - diff) - jc[y];
for( x=0; x<nrow; x++ ) {
if( mxIsFinite(*pr) ) {
*pz++ = *pr;
*iz++ = *ir;
} else {
copydata:
diff++;
}
pr++;
ir++;
}
jc[j] -= diff;
}
return jc[n];
}
From: Chad on 23 Mar 2010 09:42
Wow, thank you! This is great stuff, I appreciate the help!
From: Bruno Luong on 24 Mar 2010 07:30
"Pat Quillen" <pquillen(a)mathworks.com> wrote in message <hobruk$n97$1(a)fred.mathworks.com>...
> "Bruno Luong" <b.luong(a)fogale.findmycountry> wrote in message <ho8f2t$imm$1(a)fred.mathworks.com>...
> <snip>
> > Indeed, the bug is still not fixed in 2010A:
> >
> > >> D=sparse(1e5,1e5,Inf)
> >
> > D =
> >
> > (100000,100000) Inf
> >
> > >> D(isinf(D))=0
> >
> > D =
> >
> > (100000,100000) Inf
> <snip>
>
> Hi Bruno.
>
> This bug is fixed in R2010a in the sense that an error is thrown now instead of silently assigning into incorrect entries when trying to execute D(isinf(D)) = 0.
>
> >> version
> ans =
> 7.10.0.499 (R2010a)

Thank you Pat, mine is 2010A prereleased (7.10.0.59). Glad it has been "half corrected".

Bruno
First  |  Prev  | 
Pages: 1 2
Prev: 3d volume to 2d
Next: mex problem :: can't map file, errno=22