Буду признателен за помощь/руководство по преобразованию полной матрицы в разреженную матрицу.
Внутри mex-функции: выполнить некоторые вычисления ==> получить полную матрицу ==> Преобразовать полную матрицу в разреженную матрицу ==> выполнить еще несколько вычислений ==> вернуть результат в Matlab
1) Я не уверен, что я правильно преобразовываю полную матрицу в разреженную матрицу? Например, я получаю доступ к матрицам, используя одиночные линейные индексы, и я думаю, что это ничего не изменит в том, как матрица хранится. Я прав?
2) При преобразовании полной матрицы в разреженную матрицу (мой код основан на примере кода Matlab: fulltosparse.c) я получаю следующие ошибки:
normpdfDKU.cpp: In function ‘int Convert_Full#include "mex.h"
#include <math.h>
#include "/home/dkumar/libtsnnls-2.3.3/tsnnls/taucs_basic/taucs.h"
#include "/home/dkumar/armadillo-4.600.3/include/armadillo"
using namespace arma;
using namespace std;
#define PI (3.141592653589793)
#if defined(NAN_EQUALS_ZERO)
#define IsNonZero(d) ((d)!=0.0 || mxIsNaN(d))
#else
#define IsNonZero(d) ((d)!=0.0)
#endif
extern void _main();
const int numInputArgs = 3;
const int numOutputArgs = 2;
// Function declarations.
// -----------------------------------------------------------------
double getMatlabScalar (const mxArray* ptr);
double& createMatlabScalar (mxArray*& ptr);
int TestingLibraries() ; // declared and defined in In Include_4_TSNNLS.c and Include_4_TSNNLS.h
struct stColMajorMAt_EXT{
double* pColMajMat;
int nrows, ncols;
};
stColMajorMAt_EXT Convert_ARMA_RowMajorMAt_2_ColMajorMat(mat);
// Function definitions.
// -----------------------------------------------------------------
void mexFunction (int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]) {
// Check to see if we have the correct number of input and output
// arguments.
if (nrhs != numInputArgs)
mexErrMsgTxt(" DKU-1: Incorrect number of input arguments");
if (nlhs != numOutputArgs)
mexErrMsgTxt("DKU-2: Incorrect number of output arguments");
// Get the inputs.
double x = getMatlabScalar(prhs[0]);
double mu = getMatlabScalar(prhs[1]);
double v = getMatlabScalar(prhs[2]);
// Create the output. It is also a double-precision scalar.
double& p = createMatlabScalar(plhs[0]);
// Compute the value of the univariate Normal at x.
p = exp(-(x-mu)*(x-mu)/(2*v)) / sqrt(2*PI*v);
// CREATE ARMA::mat and print
mat B = eye<mat>(20,30);
//Print B
B.print();
mwSize sz[2];
sz[0] = B.n_rows ; // Matlab is row first
sz[1] = B.n_cols ;
//mxArray* pOUT = mxCreateNumericArray(2, sz, mxDOUBLE_CLASS, mxREAL);
plhs[1] = mxCreateNumericArray(2, sz, mxDOUBLE_CLASS, mxREAL);
//Get a pointer to pOUT
double* p2 = (double*)mxGetData(plhs[1]);
for (int i = 0; i<B.n_rows*B.n_cols; i++)
{
p2[i] = B[i];
}
// Convert B to a sparse column major matrix
/*p is already column major; but, not sparse*/
mxArray* M1 = (mxArray*)p2;
// I will convert full matrix to sparse matrix now
}
double getMatlabScalar (const mxArray* ptr) {
// Make sure the input argument is a scalar in double-precision.
if (!mxIsDouble(ptr) || mxGetNumberOfElements(ptr) != 1)
mexErrMsgTxt("The input argument must be a double-precision scalar");
return *mxGetPr(ptr);
}
double& createMatlabScalar (mxArray*& ptr) {
ptr = mxCreateDoubleMatrix(1,1,mxREAL);
return *mxGetPr(ptr);
}
mxArray* Convert_Full_2_Sparse(double* pColMajMat)
// MLoc is already a pointer to column major matrix
{
mxArray* M_Loc = (mxArray*)pColMajMat;
/* Declare variable */
mwSize m,n;
mwSize nzmax;
mwIndex *irs,*jcs,j,k;
int cmplx,isfull;
double *pr,*pi,*si,*sr;
double percent_sparse;
/* Get the size and pointers to input data */
m = mxGetM(M_Loc);
n = mxGetN(M_Loc);
pr = mxGetPr(M_Loc);
pi = mxGetPi(M_Loc);
cmplx = (pi==NULL ? 0 : 1);
/* Allocate space for sparse matrix
* NOTE: Assume at most 20% of the data is sparse. Use ceil
* to cause it to round up.
*/
percent_sparse = 0.2;
nzmax=(mwSize)ceil((double)m*(double)n*percent_sparse);
mxArray* pOUT = mxCreateSparse(m,n,nzmax,mxREAL);
sr = mxGetPr(pOUT);
si = mxGetPi(pOUT);
irs = mxGetIr(pOUT);
jcs = mxGetJc(pOUT);
/* Copy nonzeros */
k = 0;
isfull=0;
for (j=0; (j<n); j++) {
mwSize i;
jcs[j] = k;
for (i=0; (i<m ); i++) {
if (IsNonZero(pr[i])) {
/* Check to see if non-zero element will fit in
* allocated output array. If not, increase percent_sparse
* by 10%, recalculate nzmax, and augment the sparse array
*/
if (k>=nzmax){
mwSize oldnzmax = nzmax;
percent_sparse += 0.1;
nzmax = (mwSize)ceil((double)m*(double)n*percent_sparse);
/* make sure nzmax increases atleast by 1 */
if (oldnzmax == nzmax) nzmax++;
mxSetNzmax(pOUT, nzmax);
mxSetPr(pOUT, mxRealloc(sr, nzmax*sizeof(double)));
if(si != NULL) mxSetPi(pOUT, mxRealloc(si, nzmax*sizeof(double)));
mxSetIr(pOUT, mxRealloc(irs, nzmax*sizeof(mwIndex)));
sr = mxGetPr(pOUT);
si = mxGetPi(pOUT);
irs = mxGetIr(pOUT);
}
sr[k] = pr[i];
if (cmplx){
si[k]=pi[i];
}
irs[k] = i;
k++;
}
}
pr += m;
pi += m;
}
jcs[n] = k;
return pOUT;
}
Sparse(mxArray*)’:
normpdfDKU.cpp:216:69: error: invalid conversion from ‘void*’ to ‘double*’ [-fpermissive]
mxSetPr(pOUT, mxRealloc(sr, nzmax*sizeof(double)));
^
In file included from /usr/local/MATLAB/R2011b/extern/include/mex.h:58:0,
from normpdfDKU.cpp:4:
/usr/local/MATLAB/R2011b/extern/include/matrix.h:606:15: error: initializing argument 2 of ‘void mxSetPr(mxArray*, double*)’ [-fpermissive]
EXTERN_C void mxSetPr(
^
normpdfDKU.cpp:218:73: error: invalid conversion from ‘void*’ to ‘double*’ [-fpermissive]
mxSetPi(pOUT, mxRealloc(si, nzmax*sizeof(double)));
^
In file included from /usr/local/MATLAB/R2011b/extern/include/mex.h:58:0,
from normpdfDKU.cpp:4:
/usr/local/MATLAB/R2011b/extern/include/matrix.h:623:15: error: initializing argument 2 of ‘void mxSetPi(mxArray*, double*)’ [-fpermissive]
EXTERN_C void mxSetPi(
^
normpdfDKU.cpp:219:71: error: invalid conversion from ‘void*’ to ‘int*’ [-fpermissive]
mxSetIr(pOUT, mxRealloc(irs, nzmax*sizeof(mwIndex)));
^
In file included from /usr/local/MATLAB/R2011b/extern/include/mex.h:58:0,
from normpdfDKU.cpp:4:
/usr/local/MATLAB/R2011b/extern/include/matrix.h:1275:6: error: initializing argument 2 of ‘void mxSetIr_700(mxArray*, int*)’ [-fpermissive]
void mxSetIr_700(mxArray *, int *);
^
make: *** [normpdfDKU.o] Error 1
Код из файла mex:
#include "mex.h"
#include <math.h>
#include "/home/dkumar/libtsnnls-2.3.3/tsnnls/taucs_basic/taucs.h"
#include "/home/dkumar/armadillo-4.600.3/include/armadillo"
using namespace arma;
using namespace std;
#define PI (3.141592653589793)
#if defined(NAN_EQUALS_ZERO)
#define IsNonZero(d) ((d)!=0.0 || mxIsNaN(d))
#else
#define IsNonZero(d) ((d)!=0.0)
#endif
extern void _main();
const int numInputArgs = 3;
const int numOutputArgs = 2;
// Function declarations.
// -----------------------------------------------------------------
double getMatlabScalar (const mxArray* ptr);
double& createMatlabScalar (mxArray*& ptr);
int TestingLibraries() ; // declared and defined in In Include_4_TSNNLS.c and Include_4_TSNNLS.h
struct stColMajorMAt_EXT{
double* pColMajMat;
int nrows, ncols;
};
stColMajorMAt_EXT Convert_ARMA_RowMajorMAt_2_ColMajorMat(mat);
// Function definitions.
// -----------------------------------------------------------------
void mexFunction (int nlhs, mxArray *plhs[],
int nrhs, const mxArray *prhs[]) {
// Check to see if we have the correct number of input and output
// arguments.
if (nrhs != numInputArgs)
mexErrMsgTxt(" DKU-1: Incorrect number of input arguments");
if (nlhs != numOutputArgs)
mexErrMsgTxt("DKU-2: Incorrect number of output arguments");
// Get the inputs.
double x = getMatlabScalar(prhs[0]);
double mu = getMatlabScalar(prhs[1]);
double v = getMatlabScalar(prhs[2]);
// Create the output. It is also a double-precision scalar.
double& p = createMatlabScalar(plhs[0]);
// Compute the value of the univariate Normal at x.
p = exp(-(x-mu)*(x-mu)/(2*v)) / sqrt(2*PI*v);
// CREATE ARMA::mat and print
mat B = eye<mat>(20,30);
//Print B
B.print();
mwSize sz[2];
sz[0] = B.n_rows ; // Matlab is row first
sz[1] = B.n_cols ;
//mxArray* pOUT = mxCreateNumericArray(2, sz, mxDOUBLE_CLASS, mxREAL);
plhs[1] = mxCreateNumericArray(2, sz, mxDOUBLE_CLASS, mxREAL);
//Get a pointer to pOUT
double* p2 = (double*)mxGetData(plhs[1]);
for (int i = 0; i<B.n_rows*B.n_cols; i++)
{
p2[i] = B[i];
}
// Convert B to a sparse column major matrix
/*p is already column major; but, not sparse*/
mxArray* M1 = (mxArray*)p2;
// I will convert full matrix to sparse matrix now
}
double getMatlabScalar (const mxArray* ptr) {
// Make sure the input argument is a scalar in double-precision.
if (!mxIsDouble(ptr) || mxGetNumberOfElements(ptr) != 1)
mexErrMsgTxt("The input argument must be a double-precision scalar");
return *mxGetPr(ptr);
}
double& createMatlabScalar (mxArray*& ptr) {
ptr = mxCreateDoubleMatrix(1,1,mxREAL);
return *mxGetPr(ptr);
}
mxArray* Convert_Full_2_Sparse(double* pColMajMat)
// MLoc is already a pointer to column major matrix
{
mxArray* M_Loc = (mxArray*)pColMajMat;
/* Declare variable */
mwSize m,n;
mwSize nzmax;
mwIndex *irs,*jcs,j,k;
int cmplx,isfull;
double *pr,*pi,*si,*sr;
double percent_sparse;
/* Get the size and pointers to input data */
m = mxGetM(M_Loc);
n = mxGetN(M_Loc);
pr = mxGetPr(M_Loc);
pi = mxGetPi(M_Loc);
cmplx = (pi==NULL ? 0 : 1);
/* Allocate space for sparse matrix
* NOTE: Assume at most 20% of the data is sparse. Use ceil
* to cause it to round up.
*/
percent_sparse = 0.2;
nzmax=(mwSize)ceil((double)m*(double)n*percent_sparse);
mxArray* pOUT = mxCreateSparse(m,n,nzmax,mxREAL);
sr = mxGetPr(pOUT);
si = mxGetPi(pOUT);
irs = mxGetIr(pOUT);
jcs = mxGetJc(pOUT);
/* Copy nonzeros */
k = 0;
isfull=0;
for (j=0; (j<n); j++) {
mwSize i;
jcs[j] = k;
for (i=0; (i<m ); i++) {
if (IsNonZero(pr[i])) {
/* Check to see if non-zero element will fit in
* allocated output array. If not, increase percent_sparse
* by 10%, recalculate nzmax, and augment the sparse array
*/
if (k>=nzmax){
mwSize oldnzmax = nzmax;
percent_sparse += 0.1;
nzmax = (mwSize)ceil((double)m*(double)n*percent_sparse);
/* make sure nzmax increases atleast by 1 */
if (oldnzmax == nzmax) nzmax++;
mxSetNzmax(pOUT, nzmax);
mxSetPr(pOUT, mxRealloc(sr, nzmax*sizeof(double)));
if(si != NULL) mxSetPi(pOUT, mxRealloc(si, nzmax*sizeof(double)));
mxSetIr(pOUT, mxRealloc(irs, nzmax*sizeof(mwIndex)));
sr = mxGetPr(pOUT);
si = mxGetPi(pOUT);
irs = mxGetIr(pOUT);
}
sr[k] = pr[i];
if (cmplx){
si[k]=pi[i];
}
irs[k] = i;
k++;
}
}
pr += m;
pi += m;
}
jcs[n] = k;
return pOUT;
}
