Hello everyone. I'm still working with sparse matrices in oct-file, but i've got some problems: i'm doing some multithread and as the octave libraries revealed to be not really thread safe (some unexpected crashed ofter repeating executions) i did the dirty work by hand: i copied the 3 vectors of the matrix B (the input matrix) and copied the full matrix A into a vector. Then did what i have to do and then gave back the 3 vectors for the matrix C, then create the sparsematrix out of the 3 vectors and return from the function. something in creating the c sparse matrix out of the 3 vector is going wrong, and the documentation about this seem to be poor. For example, why the 3 vector has to be the same lenght? the cidx is always really shorter than the other 2. (i solved the problem filling with zeroes but i don't know if it's the correct why). Another problem is that sometimes when creating the matrix it says i'm referencing an item that's out of the matrix bounds (even if the 3 vectors have correct indexes and values inside i suppose). The last problems is that i get an error on a free call, the log is long so i post an attachment with it. the call is:
mkoctfile -v fl_spcompose.cc A = sprand(10, 10, 0.5); B = sprand(10, 10, 0.5); tic; C=fl_spcompose(A,B); toc Hope that someone can help me solving my problem (i'll attach the full sourcecode too, the part i'm talking about is from 319 to 371)
#include <octave/oct.h>
#include <octave/parse.h>
#include <pthread.h>
#define HELP " -- Loadable Function: C = fl_spcompose (A, B)\n -- Loadable Function: C = fl_spcompose (A, B, LOCK)\n -- Loadable Function: C = fl_spcompose (A, B, T)\n -- Loadable Function: C = fl_spcompose (A, B, T, S)\n -- Loadable Function: C = fl_spcompose (A, B, LOCK, T)\n\nReturns the T-Norm / S-Norm composition as basic inference mechanism of Fuzzy\nLogic. By default, it calculates the max-min composition.\n\nA and B must be matrices with conformant dimensions as in matrix product.\n\nWhen true, the boolean LOCK option forces to calculate the diagonal results\nonly and returns it as a column vector.\n\nThe arguments T and S allows you to specify a custom T-Norm and S-Norm function\nrespectively. They can be:\n - 'min': use the minimum function (default for T-Norm);\n - 'pro': use the product function;\n - 'max': use the maximum function (default for S-Norm);\n - 'sum': use the sum function;\n - function_handle: a user-defined function (at most 2 arguments).\n\nNote that only the predefined functions are calculated rapidly and in\nmultithread mode. Using a user-defined function as T-Norm and/or S-Norm\nwill result in a long time calculation. \n\n\n\n\n"
#define MIN_ARG 2
#define MAX_ARG 4
#define STR_MINNORM "min"
#define STR_PRONORM "pro"
#define STR_MAXNORM "max"
#define STR_SUMNORM "sum"
#define ERR_INVALIDFUNC "fl_spcompose: the specified T-Norm or S-Norm function is invalid!\nPlease read help for details"
#define ERR_ARG "fl_spcompose: invalid arguments number!\n"
#define ERR_MATRIXSIZE "fl_spcompose: incompatible matrices dimensions\n"
#define ERR_MATRIXTYPE "fl_spcompose: the first two argument must be matrices\n"
#define ERR_LOCKOPTION "fl_spcompose: the lock_option must be 0 or 1\n"
// Structure for thread arguments. Each thread will perform the computation between the start_index and end_index row.
struct threadArg
{
int start_index;
int end_index;
};
// Functions prototype declaration
double getElem(double vec[], int row, int col,int numCols);
int checkFunction(octave_function *func);
void *thread_function(void *arg);
int getCpuCore();
void compose();
int assignFunction(octave_value arg, int tnorm);
double (*calc_tnorm)(double,double);
double (*calc_snorm)(double,double);
double func_min(double first, double second);
double func_prod(double first, double second);
double func_max(double first, double second);
double func_sum(double first, double second);
double func_custom_tnorm(double first, double second);
double func_custom_snorm(double first, double second);
// T-Norm and S-Norm function pointers
octave_function *tnorm;
octave_function *snorm;
// Structure for the thread arguments
struct threadArg* thread_args = NULL;
// Input and output matrices
SparseMatrix a,b,c;
// Matrix component vectors
int *cidxB, *ridxB, *cidxC, *ridxC;
double *dataB, *dataC;
double *copyA;
// Matrices dimensions
////////////////////////////////////////////////////////
int rowsA,rowsB,colsA,colsB,nnzA,nnzB;
// Lock option. 1 = calculation executed only for the diagonal of the matrix
int lock_option = 0;
// The increment
int col_index_increment;
// Number of threads that will be created
int num_threads;
/* Main function - Check the arguments and call the compose method */
DEFUN_DLD (fl_spcompose, args, nargout, HELP)
{
int numargs = args.length(); // Arguments number
int specified_lockoption = 0; // 1 = the lock_option was specified in the arguments
// Set the num_thread to the actual available CPU cores
num_threads = getCpuCore();
// Set the lock_option to default value (0)
lock_option = 0;
// Set the T-Norm and S-Norm function pointer to default value (Minimum and Maximum)
calc_tnorm = func_min;
calc_snorm = func_max;
// Check if the argument number is correct
if(numargs < MIN_ARG || numargs > MAX_ARG)
{
error(ERR_ARG);
return octave_value_list();
}
// Check if the first two arguments are matrices
if(!args(0).is_matrix_type() || !args(1).is_matrix_type())
{
error(ERR_MATRIXTYPE);
return octave_value_list();
}
// Analyze the third argument (can be the lock_option, a custom T-Norm or a default string for T-Norm)
if(numargs > 2)
{
// Check if the third argument is the lock_option.
if(args(2).is_bool_scalar())
{
// Get the lock_option value
lock_option = args(2).int_value();
specified_lockoption = 1;
}
// Check if the third argument is a function
else if(args(2).is_function_handle())
{
// Check if the custom T-Norm function has at most two arguments
if(!checkFunction(args(2).function_value()))
{
error(ERR_INVALIDFUNC);
return octave_value_list();
}
else
{
// Set the custom T-Norm and force single thread mode
tnorm = args(2).function_value();
calc_tnorm = func_custom_tnorm;
num_threads = 1;
}
}
// Check if the third argument is a string and try to assign a default T-Norm
else if(!((args(2).is_sq_string() || args(2).is_dq_string()) && assignFunction(args(2),1)))
{
error(ERR_INVALIDFUNC);
return octave_value_list();
}
}
// Analyze the fourth argument (can be a custom T-Norm, a default string for T-Norm or a custom S-Norm)
if(numargs > 3)
{
// Check if the fourth argument is a function
if(args(3).is_function_handle())
{
if(specified_lockoption)
{
// Check if the custom T-Norm function has at most two arguments
if(!checkFunction(args(3).function_value()))
{
error(ERR_INVALIDFUNC);
return octave_value_list();
}
else
{
// Set the custom T-Norm and force single thread mode
tnorm = args(3).function_value();
calc_tnorm = func_custom_tnorm;
num_threads = 1;
}
}
else
{
// Check if the custom S-Norm function has at most two arguments
if(!checkFunction(args(3).function_value()))
{
error(ERR_INVALIDFUNC);
return octave_value_list();
}
else
{
// Set the custom S-Norm and force single thread mode
snorm = args(3).function_value();
calc_snorm = func_custom_snorm;
num_threads = 1;
}
}
}
// Check if the fourth argument is a string (could be a default T-Norm)
else if (args(3).is_sq_string() || args(3).is_dq_string())
{
if(specified_lockoption)
{
if (!assignFunction(args(3),1))
{
error(ERR_INVALIDFUNC);
return octave_value_list();
}
}
else
{
if (!assignFunction(args(3),0))
{
error(ERR_INVALIDFUNC);
return octave_value_list();
}
}
}
else
{
error(ERR_INVALIDFUNC);
return octave_value_list();
}
}
// Get the input matrices
a = args(0).sparse_matrix_value();
b = args(1).sparse_matrix_value();
// Get the dimensions of matrices
rowsA = a.dims()(0);
colsA = a.dims()(1);
nnzA = a.nnz();
rowsB = b.dims()(0);
colsB = b.dims()(1);
nnzB = b.nnz();
const SparseMatrix tmpA (a);
copyA = new double[rowsA*colsA];
for (int i=0; i<rowsA; i++)
for(int j=0; j<colsA;j++)
copyA[i*colsA+j] = tmpA(i,j);
int* tmp;
int lenght;
tmp = b.cidx();
lenght = sizeof(tmp)/sizeof(*tmp);
cidxB = new int[colsB+1];
octave_stdout << "colsB= " << colsB << "\n";
for (int i = 0; i <= colsB+1; i++) {
octave_stdout << "b.cidx(" << i << ")= " << b.cidx(i) << "\n";
cidxB[i] = b.cidx(i);
}
tmp = b.ridx();
lenght = sizeof(tmp)/sizeof(*tmp);
ridxB = new int[nnzB];
octave_stdout << "nnzB= " << nnzB << "\n";
for (int i = 0; i < nnzB; i++) {
octave_stdout << "b.ridx(" << i << ")= " << b.ridx(i) << "\n";
ridxB[i] = b.ridx(i);
}
double* tmpData;
/*
tmpData = a.data();
lenght = sizeof(tmpData)/sizeof(*tmpData);
dataA = new double[lenght];
for (int i = 0; i < lenght; i++)
dataA[i] = a.data(i);
*/
tmpData = b.data();
lenght = sizeof(tmpData)/sizeof(*tmpData);
dataB = new double[nnzB];
for (int i = 0; i < nnzB; i++)
dataB[i] = b.data(i);
// Check if the dimensions are compatible
if(colsA != rowsB)
{
error(ERR_MATRIXSIZE);
return octave_value_list();
}
// If the lock_option is true, the function will return a column vector.
// Sets proper matrix dimension and column index increment.
if(lock_option)
{
// A column index increment = colsB let the two FOR cycle to calculate only the first element
col_index_increment = colsB;
//c = SparseMatrix(1, rowsA, rowsA);
cidxC = new int[1];
ridxC = new int[rowsA];
dataC = new double[rowsA];
}
else
{
col_index_increment = 1;
//c = SparseMatrix(colsB, rowsA, colsB*rowsA);
cidxC = new int[colsB];
ridxC = new int[colsB*rowsA];
dataC = new double[colsB*rowsA];
}
// Start the matrix composition
compose();
printf("res\n");
if(lock_option)
{
ColumnVector cv_dataC (rowsA);
octave_stdout << "rowsA= " << rowsA << "\n";
for (int i = 0; i < rowsA; i++) {
octave_stdout << "dataC{" << i << "}= " << dataC[i] << "\n";
cv_dataC(i) = dataC[i];
}
ColumnVector cv_ridxC (rowsA);
for (int i = 0; i < rowsA; i++)
cv_ridxC(i) = ridxC[i];
ColumnVector cv_cidxC (rowsA+1);
for (int i = 0; i < rowsA+1; i++)
cv_cidxC(i) = cidxC[i];
c = SparseMatrix(cv_dataC, cv_ridxC, cv_cidxC, 1, rowsA);
}
else
{
ColumnVector cv_cidxC (colsB+1);
for (int i = 0; i < colsB+1; i++) {
octave_stdout << "cidxC{" << i << "}= " << cidxC[i] << "\n";
cv_cidxC(i) = cidxC[i];
}
for (int i = colsB+1; i < colsB+1; i++) {
cv_cidxC(i) = 0;
}
ColumnVector cv_dataC (cidxC[colsB]);
octave_stdout << "colsB*rowsA= " << cidxC[colsB] << "\n";
for (int i = 0; i < cidxC[colsB]; i++) {
octave_stdout << "dataC{" << i << "}= " << dataC[i] << "\n";
cv_dataC(i) = dataC[i];
}
ColumnVector cv_ridxC (cidxC[colsB]);
for (int i = 0; i < cidxC[colsB]; i++) {
octave_stdout << "ridxC{" << i << "}= " << ridxC[i] << "\n";
cv_ridxC(i) = ridxC[i];
}
c = SparseMatrix(cv_dataC, cv_ridxC, cv_cidxC, colsB, rowsA, false);
}
// Return the resulting matrix
c.maybe_compress();
return octave_value(c.transpose());
}
// Calculate the S-Norm/T-Norm composition
void compose()
{
// Create the thread args array
thread_args = new threadArg[num_threads];
// Define an array of threads
pthread_t th[num_threads];
// Define the number interval of rows for each thread
int interval = rowsA / num_threads;
printf("thread");
// Define the threads
for (int i=0; i<num_threads; i++)
{
// Set the proper row start_index and end_index in the thread argument
thread_args[i].start_index = i*interval;
thread_args[i].end_index = thread_args[i].start_index + interval;
if(i == num_threads - 1)
thread_args[i].end_index = rowsA;
// Start the thread
pthread_create(&th[i],NULL,thread_function, (void *) &thread_args[i]);
}
void *ans1;
// Wait the results from each thread
for(int i=0; i<num_threads; i++)
pthread_join(th[i],&ans1);
}
/* Function executed by each thread */
void *thread_function(void *arg)
{
// Get the structure from the thread arguments
struct threadArg *thread_args;
thread_args = (struct threadArg *) arg;
// Create constat versions of the input matrices to prevent them to be filled by zeros on reading
const SparseMatrix tmpA (a);
const SparseMatrix tmpB (b);
// Declare variables for the T-Norm and S-Norm values
double snorm_val;
double tnorm_val;
// Get the row start_index and end_index
int start_index = thread_args->start_index;
int end_index = thread_args->end_index;
// Initialize the number of nonzero elemnts in sparse matrix c
int nel = 0;
cidxC[0] = 0;
printf("for\n");
// Calculate the composition for the specified rows (between start_index and end_index)
for (int i = start_index; i < end_index; i++)
{
for(int j = lock_option*i; j < colsB; j = j+col_index_increment)
{
snorm_val = 0;
printf("nelforj\n");
// From b.cidx(j) to b.cidx(j+1)-1 there are all the nonzero elements of the column j of b (tmpB actually)
for (int k = cidxB[j]; k <= (cidxB[j+1]-1); k++)
{
printf("nelfork\n");
octave_stdout << "k= " << k << " cidxB[j+1]= " << cidxB[j+1] << "\n";
// b.ridx(k) is the row index of the actual nonzero element in b (tmpB actually)
// so it is also the column index of the element in a (tmpA actually)
tnorm_val = calc_tnorm(getElem(copyA, i, ridxB[k], colsA), dataB[k]);
snorm_val = calc_snorm(snorm_val, tnorm_val);
}
if (snorm_val != 0)
{
printf("aggiunto\n");
ridxC[nel] = j;
dataC[nel++] = snorm_val;
}
}
cidxC[i+1] = nel;
}
return((void *)0);
}
/* Returns 1 if an octave_function FUNC has at most 2 arguments. */
int checkFunction(octave_function *func)
{
octave_value_list testArgs;
testArgs(0) = 1;
testArgs(1) = 2;
feval(func,testArgs);
if(error_state)
return 0;
else
return 1;
}
/* Assign a default function or T-Norm or S-Norm.
If tnorm = 0 the function will be assigned as S-Norm, else the function will be assigned as T-Norm.
Returns 0 if it fails*/
int assignFunction(octave_value arg, int tnorm)
{
int res = 1;
if(arg.string_value() == STR_PRONORM)
if(tnorm)
calc_tnorm = func_prod;
else
calc_snorm = func_prod;
else if (arg.string_value() == STR_MAXNORM)
if(tnorm)
calc_tnorm = func_max;
else
calc_snorm = func_max;
else if (arg.string_value() == STR_SUMNORM)
if(tnorm)
calc_tnorm = func_sum;
else
calc_snorm = func_sum;
else if(!(arg.string_value() == STR_MINNORM))
res = 0;
return res;
}
/* Calculate the minimum between the two values */
double func_min(double first, double second)
{
if(first < second)
return first;
else
return second;
}
/* Calculate the product between the two values */
double func_prod(double first, double second)
{
return first*second;
}
/* Calculate the Maximum between the two values */
double func_max(double first, double second)
{
if(first > second)
return first;
else
return second;
}
/* Calculate the sum between the two values */
double func_sum(double first, double second)
{
return first+second;
}
/* Calculate a custom T-Norm between the two values */
double func_custom_tnorm(double first, double second)
{
octave_value_list fargs;
fargs(0) = octave_value(first);
fargs(1) = octave_value(second);
return feval(tnorm,fargs)(0).double_value();
}
/* Calculate a custom S-Norm between the two values */
double func_custom_snorm(double first, double second)
{
octave_value_list fargs;
fargs(0) = octave_value(first);
fargs(1) = octave_value(second);
return feval(snorm,fargs)(0).double_value();
}
/* Get the current cpu cores number */
int getCpuCore()
{
#ifdef WIN32
SYSTEM_INFO sysinfo;
GetSystemInfo( &sysinfo );
return sysinfo.dwNumberOfProcessors;
#else
return sysconf( _SC_NPROCESSORS_ONLN);
#endif
}
/* Get the (i,j)-th element from the vector vec. The column number of the original matrix (numCols) is required */
double getElem(double vec[], int row, int col,int numCols)
{
return vec[row*numCols+col];
}
log
Description: Binary data
------------------------------------------------------------------------------ Learn how Oracle Real Application Clusters (RAC) One Node allows customers to consolidate database storage, standardize their database environment, and, should the need arise, upgrade to a full multi-node Oracle RAC database without downtime or disruption http://p.sf.net/sfu/oracle-sfdevnl
_______________________________________________ Octave-dev mailing list Octave-dev@lists.sourceforge.net https://lists.sourceforge.net/lists/listinfo/octave-dev
