Search All of the Math Forum:

Views expressed in these public forums are not endorsed by NCTM or The Math Forum.

Topic: second largest element in a matrix
Replies: 32   Last Post: May 5, 2013 10:31 AM

 Search Thread: Advanced Search

 Messages: [ Previous | Next ]
 Bruno Luong Posts: 9,822 Registered: 7/26/08
Re: second largest element in a matrix
Posted: Apr 6, 2009 8:04 PM
 Plain Text Reply

I have implemented the MIN by partial quicksort in MEX. Any comment is welcome. As soon as I get the MAX going, it will be on FEX.

Thanks,

Bruno

/**************************************************************************
* function [res loc] = minkmex(list, k)
* Matlab C-Mex
* Purpose: Same as MINK, i.e.,
* Return in RES the K smallest elements of LIST
* RES is sorted in ascending order [res loc] = mink(...)
* Location of the smallest: RES=LIST(LOC)
* But this MEX works on double, and output RES is unsorted
* Algorithm according to http://en.wikipedia.org/wiki/Selection_algorithm
* Compilation: mex -O -v minkmex.c
* Author Bruno Luong <brunoluong@yahoo.com>
* Last update: 06/April/2009
*************************************************************************/

#include "mex.h"
#include "matrix.h"

#define _DEBUG

/* Global variables, used to avoid stacking them during recusive call since
they do not change */
mwSize k;
mwSize *pos;
double *list;

/* Partitioning the list around pivot pivotValue := l[pivotIndex];
* After runing, at exit we obtain:
l[left]...l[index-1] < pivotValue <= l[index] ... l[right]
where l[i] := list[pos[i]] for all i */
mwSize partition(mwSize left, mwSize right, mwSize pivotIndex) {

double pivotValue;
mwSize *pindex, *pi, *pright;
mwSize tmp;

pright=pos+right;
pindex=pos+pivotIndex;
pivotValue = list[tmp = *pindex];
/* Swap pos[pivotIndex] with pos[right] */
*pindex = *pright;
*pright = tmp;

pindex=pos+left;
for (pi=pindex; pi<pright; pi++)
/* Compare with pivotValue */
if (list[*pi] < pivotValue) {
/* if smaller; Swap pos[index] with p[i] */
tmp = *pindex;
*pindex = *pi;
*pi = tmp;
pindex++;
}

/* Swap pos[index] with p[right] */
tmp = *pindex;
*pindex = *pright;
*pright = tmp;

return (pindex-pos); /* Pointer arithmetic */
} /* Partition */

/* Recursive engine (partial quicksort) */
void findFirstK(mwSize left, mwSize right) {

mwSize pivotIndex;

if (right > left) {
pivotIndex = partition(left, right, (left+right+1)/2);
if (pivotIndex > k)
findFirstK(left, pivotIndex-1);
else if (pivotIndex < k)
findFirstK(pivotIndex+1, right);
}

return;
} /* findFirstK */

/* Gateway of minkmex */
void mexFunction(int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]) {

mwSize l, i;
double *data;
mwSize dims[2];

/* Check arguments */
if (nrhs<2)
mexErrMsgTxt("MINKMEX: Two input arguments required.");

if (!mxIsNumeric(prhs[0]))
mexErrMsgTxt("MINKMEX: First input LIST argument must be numeric.");

if (!mxIsNumeric(prhs[1]))
mexErrMsgTxt("MINKMEX: Second input K must be numeric.");

/* Get the number of elements of the list of subindexes */
if (mxGetM(prhs[0])==1)
l = mxGetN(prhs[0]);
else if (mxGetN(prhs[0])==1)
l = mxGetM(prhs[0]);
else
mexErrMsgTxt("MINKMEX: First input LIST must be a vector.");

if (mxGetClassID(prhs[0]) != mxDOUBLE_CLASS)
mexErrMsgTxt("MINKMEX: First input LIST must be a double.");

/* Get the number of elements of the list of subindexes */
if (mxGetM(prhs[1])!=1 || mxGetN(prhs[1])!=1)
mexErrMsgTxt("MINKMEX: Second input K must be a scalar.");

if (mxGetClassID(prhs[1]) != mxDOUBLE_CLASS)
mexErrMsgTxt("MINKMEX: Second input K must be a double.");

k = (int)(*mxGetPr(prhs[1]));
if (k<0)
mexErrMsgTxt("MINKMEX: K must be non-negative integer.");

/* Get a data pointer */
list = mxGetPr(prhs[0]);

/* Clip k */
if (k>l) k=l;

/* Clean programming */
pos=NULL;

/* Work for non-empty array */
if (l>0) {
/* Vector of index */
pos = mxMalloc(sizeof(mwSize)*l);
if (pos==NULL)
mexErrMsgTxt("Out of memory.");
/* Initialize the array of position (zero-based index) */
for (i=0; i<l; i++) pos[i]=i;

/* Call the recursive engine */
k--; /* because we work on zero-based */
findFirstK(0, l-1);
k++; /* Restore one-based index */
} /* if (l>0) */

/* Create the Matrix result (first output) */
dims[0] = 1; dims[1] = k;
plhs[0] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
if (plhs[0] == NULL)
mexErrMsgTxt("Out of memory.");
data = mxGetPr(plhs[0]);
for (i=0; i<k; i++) data[i]=list[pos[i]];

/* Create the Matrix position (second output) */
if (nlhs>=2)
{
dims[0] = 1; dims[1] = k;
plhs[1] = mxCreateNumericArray(2, dims, mxDOUBLE_CLASS, mxREAL);
if (plhs[1] == NULL)
mexErrMsgTxt("Out of memory.");
data = mxGetPr(plhs[1]);
for (i=0; i<k; i++) data[i]=(double)(pos[i]+1); /* one-based */
}

/* Free the array of position */
if (pos) mxFree(pos);
pos = NULL; /* clean programming */

return;

} /* Gateway of minkmex.c */

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [res loc] = mink(list, k)
% function res = mink(list, k)
%
% Return in RES the K smallest elements of LIST
% RES is sorted in ascending order
% [res loc] = mink(...)
% Location of the smallest: RES=LIST(LOC)
%
% Author Bruno Luong <brunoluong@yahoo.com>
% Last update: 06/April/2009

clist=class(list);
% Mex file requires double
if ~strcmpi(clist,'double')
list=double(list);
end
k=double(k);

[res loc] = minkmex(list(:),k);
[res is] = sort(res);
loc = loc(is);

% Cast to original class
if ~strcmpi(clist,'double')
res=feval(clist,res);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Test on command line

clear
k=10;
n=1e6;
ntest=50;
tmink=zeros(1,ntest);
tsort=zeros(1,ntest);
for i=1:ntest
list=rand(1,n);

tic
m=mink(list,k);
tmink(i)=toc;

tic
s=sort(list);
s=s(1:k);
tsort(i)=toc;

if ~isequal(m,s)
keyboard;
end
end
disp('mink')
disp(median(tmink))
disp('sort:')
disp(median(tsort))

Date Subject Author
9/18/08 Oluwa KuIse
9/19/08 stephanie
9/19/08 Steven Lord
9/19/08 Oluwa KuIse
9/19/08 Pekka
9/19/08 Oluwa KuIse
9/21/08 Greg Heath
9/24/08 Walter Roberson
3/4/09 Justin
4/6/09 Bruno Luong
4/6/09 Jos
4/6/09 Bruno Luong
4/6/09 Bruno Luong
4/6/09 mike zander
4/7/09 Bruno Luong
4/7/09 Bruno Luong
4/3/09 Paul
11/7/12 sssbi2009@gmail.com
11/19/12 dpb
4/3/09 Paul
4/3/09 mike zander
4/3/09 mike zander
4/3/09 mike zander
4/3/09 Walter Roberson
4/3/09 mike zander
5/4/13 tilindg1@gmail.com
5/4/13 Nasser Abbasi
5/4/13 dpb
5/5/13 Nasser Abbasi
5/5/13 dpb
5/5/13 Bruno Luong
5/5/13 Nasser Abbasi
5/5/13 Bruno Luong

© The Math Forum at NCTM 1994-2017. All Rights Reserved.