2011/1/7 David Bateman <dbate...@dbateman.org>:
> w4nderlust wrote:
>> Il giorno 06/gen/2011, alle ore 22.27, David Bateman ha scritto:
>>
>>
>>> w4nderlust wrote:
>>>
>>>> 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.
>>>>
>>> If you do something like
>>>
>>> const SparseMatrix a = args(0).sparse_matrix_value ();
>>> const *double adata = a.data();
>>> const *octave_idx_type acidx = a.cidx();
>>> const *octave_idx_type aridx = a.ridx();
>>>
>>> you can clean up a lot of the issues I believe you were having. As the
>>> rest of you code will then just work on these pointers there is no
>>> reason that the code shouldn't be thread-safe. However
>>>
>>
>> I'm trying right now, but when i was using the multithread on the nonsparse
>> version of this project (so using the octve matrix objects) it crashed when
>> repeating 1000 executions continuously, but when i modifidies to use simble
>> double * as matrix class replacement it worked fine. That's why i'm doing
>> all this :)
>>
> I think you missed my point. The above gives you double* and
> octave_idx_type* pointer to the underlying data structures of the Octave
> sparse matrix that you can manipulate directly. No need to copy from the
> matrix to a *double pointer.
Yes it's clear, i was explaining way i choosed to work this way insted
of using the class.
>
>> I think you are talking about this:
>> http://www.network-theory.co.uk/docs/octave3/octave_204.html 20.1.1
>> I was talking about the oct-file part (
>> http://www.gnu.org/software/octave/doc/interpreter/Sparse-Matrices-in-Oct_002dFiles.html
>> ).
>
> The two sections were written together. If you don't know how sparse
> matrices are used in the octave interpreter its perhaps permature to
> read the section on using them in oct-files. So section 21 of the manual
> should be read before the appendix.
>> I think is lacking in at least 2 points, i write them here not to contrast
>> your opinion, but simply because maybe someone would agree and improve the
>> documentation (maybe when i finish this project i could do something for
>> what i think i understand well
>>
> I wrote these sections, and although these were written rapidly, if they
> aren't read independently, and other basic concepts of oct-files are
> understood they should be understandable. Yes they might be improved,
> but I think you need to use sparse matrices a bit more in Octave first
> as there seems to be some basic points you're still missing.
I had to work in C++ from the principle so i really didn't care, now i
red it all.
>
>> - the example in 20.1.1 is confusing with the example in
>> http://www.gnu.org/software/octave/doc/interpreter/Creating-Sparse-Matrices-in-Oct_002dFiles.html#Creating-Sparse-Matrices-in-Oct_002dFiles
>> (the first one) as the cidx vector is completely different 8and i think it
>> is wrong here).
>>
> What is proposed in 20.1.2 (21.1.2 in the manual of 3.3.x and soon to be
> 3.4.0) rather than 20.1.1. The data structure cidx, ridx of a sparse
> matrix are not exposed in the interpreter of Octave and so working with
> them directly is not an option. What is proposed in section 2.1.2 is to
> create 3 vectors representing the triplets (row, col, val) that are then
> internally converted to cidx, ridx, data in an fairly efficient manner.
> In oct-files you have access to cidx, ridx, data and working directly
> with them is more efficient and trying to work on a sparse matrix in an
> oct-file directly with (row, col, val) triplets is going to be extremely
> slow.. This explains the difference between the two.
>
> If you want to do the same thing that the example in 21.1.2 suggests in
> a oct-file then do something like
>
> OCTAVE_LOCAL_BUFFER (double, val, nnz)
> OCTAVE_LOCAL_BUFFER (octave_idx_type, rows, nnz)
> OCTAVE_LOCAL_BUFFER (octave_idx_type, cols, nnz)
>
> for (octave_idx_type i = 0; i < nnz; i++)
> {
> val[i] = ....;
> rows[i] = ....;
> cols[i] = ....;
> }
>
> SparseMatix c (rows, cols, val)
>
Ok now it's really cleared, thankyou (and that also clarifies why they
have to be the same lenght). Probably if you add this clear
description to the doc, it would help someone else.
> instead. However working directly with cidx, ridx and data will save you
> a copying the matrix in memory and will be more efficient.
>> - the usage page is really poor, the only method shown
>> sparse_bool_matrix_value () .
> I have no idea what you mean by "usage page".. I presume you mean
> section A1.6.3. This section as it says is how to use the sparse matrix
> type in oct-files themselves and concerns only the interface from the
> octave_value type and back again so that you can deal with the data.
> Section A1.6.2 had already discussed how to deal with the data in the
> sparse matrices themselves.. Section A1.6.3 has the example
>
> octave_value_list retval;
> SparseMatrix sm = args(0).sparse_matrix_value ();
> SparseComplexMatrix scm =
> args(1).sparse_complex_matrix_value ();
> SparseBoolMatrix sbm = args(2).sparse_bool_matrix_value ();
> ...
> retval(2) = sbm;
> retval(1) = scm;
> retval(0) = sm;
>
> and all three types of sparse matrices are treated.. I can accept
> criticism but at least make it valid.
The point is that this doesn't explain how to use the class methods,
That would really be needed, at least i think so, you are free to
disagree :)
>
>> All he other methods are present in the octave-forge documentation but
>> there's no description for each one and some parameters are really obscure
>> to me (for example the boolean value at the end of the constructor with the
>> 3 vectors)
>>
>>
> I presume by the octave-forge documentation you mean the doxygen pages at
>
> http://octave.sourceforge.net/doxygen/html/class_sparse_matrix.html
>
> This is created from the octave source code directly by doxygen and has
> no more information than the source code itself. By boolean value at the
> end of the constructor I presume you mean the constructor given by
>
> http://octave.sourceforge.net/doxygen/html/class_sparse_matrix.html#3d38fab78c54c461739f46ae8e7196ff
>
> where the boolean value is called "sum_terms". Nothing stops you
> (row,col,val) triplets from duplicating (row,col) values. In this case
> you have two choices. Keep only the second value or sum the values with
> the same (row,col) values which is the default.. Its kind of implicit in
> the variable name "sum_terms" I would have thought.. Or you might have
> figured it out by reading "help sparse" in the section
>
> <quote>
> Note: if multiple values are specified with the same I, J indices, the
> corresponding values in S will
> be added.
>
> The following are all equivalent:
>
> s = sparse (i, j, s, m, n)
> s = sparse (i, j, s, m, n, "summation")
> s = sparse (i, j, s, m, n, "sum")
>
> Given the option "unique". if more than two values are specified for the
> same I, J indices, the last specified value will be used.
> </quote>
>
> Sure it might be clearer, but no one has complained about this point
> yet. The octave code is pretty clear and so you shouldn't hesitate to
> jump into it directly for more info.
It was just an example to show that not really everything is crystal
clear, that was really obscure to me. Probably it's my fault.
>
>>
>>> I see two solutions.
>>>
>>> 1) Precalculate the number of non-zero elements per-column. That is, in
>>> the code block I sent you for spmaxmin_new.cc calculate the output cidx
>>> at the same time you calculate nel, rather than in the second loop. If
>>> you first save the number of non-zero elements to ccidx you can even
>>> thread this. For example (with code that might be threaded, without the
>>> threading itself)
>>>
>>> c.xcidx (0) = 0;
>>> for (octave_idx_type i = 0; i < rowsA; i++)
>>> {
>>> // From here on this code could be threaded
>>> octave_idx_type kk = 0
>>>
>>> for(octave_idx_type j = 0; j < colsB; j++)
>>> {
>>> if (b.cidx(j) != b.cidx(j+1))
>>> kk++;
>>> else
>>> {
>>> for (octave_idx_type k = b.cidx(j);
>>> k <= (b.cidx(j+1) - 1); k++)
>>> {
>>> if (a(i, b.ridx(k)) != 0)
>>> {
>>> kk++;
>>> break;
>>> }
>>> }
>>> }
>>> }
>>> c.xcidx(i+1) = kk;
>>> }
>>>
>>> // Now that we've calculated the number of elements per column of
>>> the output matrix in code
>>> // that might be threaded, correct the indexing (Can't be threaded)
>>> for (octave_idx_type i = 0; i < rowsA; i++)
>>> c.xcidx(i+1) += c.xcidx(i);
>>>
>>> octave_idx_type nel = c.xcidx(rowsA);
>>>
>>> It might not be profitable to thread the above and it could be
>>> simplified to
>>>
>>> octave_idx_type nel = 0;
>>>
>>> c.xcidx(0) = 0;
>>> for (octave_idx_type i = 0; i < rowsA; i++)
>>> {
>>> for(octave_idx_type j = 0; j < colsB; j++)
>>> {
>>> if (b.cidx(j) != b.cidx(j+1))
>>> nel++;
>>> else
>>> {
>>> for (octave_idx_type k = b.cidx(j);
>>> k <= (b.cidx(j+1) - 1); k++)
>>> {
>>> if (a(i, b.ridx(k)) != 0)
>>> {
>>> nel++;
>>> break;
>>> }
>>> }
>>> }
>>> }
>>> c.xcidx(i+1) = nel;
>>> }
>>>
>>> You can then uses c.xcidx in your thread_function to index correctly
>>> into c.xridx and c.xdata.
>>>
>>>
>>> 2) You could create N different sparse matrices, one per thread, and
>>> then concatenate them together with the concat method. This will
>>> probably be slower and more memory hungry than the above so I won't
>>> develop this idea further.
>>>
>>>
>>
>> I think the first is the wa i will follow. Mabe i'll try to execute my test
>> adding directly to c(i,j) for each thread (this makes the code really
>> cleaner for the threads) and not working directly on the c index vectors. If
>> this still gives me creashes on multiple executions, i'll try to wark the
>> way you described on point 1.
>>
>
> I reaIly wouldn't suggest working directly with c(i,j) as writing to
> c(i,j) will essentially move the data in the matrix to to make way for
> the new data. If you want to work with c(i,j) then the loops should be
> in this order
>
> for (octave_idx_type i = 0; i < m; i++)
> for (octave_idx_type j = 0; j < n; j++)
> c(i,j) = ....;
>
> and not the reverse or any out of order row,col calculation as this will
> allow the new elements to be written directly to the end of the ridx and
> data arrays. However, to just simply search the matrix to ensure that
> such a copy isn't needed and modifying all the elements of cidx implied
> by a write to c(i,j) is a O(n) operation, whereas working directly with
> an element using cidx,ridx,data should be O(1). So I know which way I'd
> prefer to write my code.
>
> D.
>
>
That's the way i'm working right now.
Really not right now, because i discovered a bug in the single thread
version and i'm correcting it at this time. Can you help me about it?
The problem is this: if the result matrix is a vector 1xn or nx1 the
final c.transpose() gives a segmentation fault. I really can't
understand why.
You can easly reproduce the problem addint a true as third parameter
to the function (it calculates only the diagonal values and returns it
in a colum vector).
Like:
A = sprand(400, 400, 0.02);
B = sprand(400, 400, 0.02);
Z= fl_spcompose(A, B, true);
I attach the code
Thanks anyway ;)
P.S. My criticism is just a description of a problem i encountered
hoping that someone alse would not encounter them, not a critisism
against you or the octave project
#include <octave/oct.h>
#include <octave/parse.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"
// Functions prototype declaration
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);
int checkFunction(octave_function *func);
int assignFunction(octave_value arg, int tnorm);
// T-Norm and S-Norm function pointers
octave_function *tnorm;
octave_function *snorm;
// Input and output matrices
SparseMatrix a,b,c;
// 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;
/* 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 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;
}
}
// 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;
}
}
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;
}
}
}
// 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();
// Create constat versions of the input matrices to prevent them to be filled by zeros on reading
const SparseMatrix constA (a);
const SparseMatrix constB (b);
// 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();
// Check if the dimensions are compatible
if(colsA != rowsB)
{
error(ERR_MATRIXSIZE);
return octave_value_list();
}
// Count the number of non-zero elements first. This should be
// relatively low cost
/*
octave_idx_type nnzC = 0;
for (octave_idx_type i = 0; i < rowsA; i++)
{
for(octave_idx_type j = 0; j < colsB; j++)
{
if (b.cidx(j) != b.cidx(j+1))
nnzC++;
else
{
for (octave_idx_type k = b.cidx(j);
k <= (b.cidx(j+1) - 1); k++)
{
if (a(i, b.ridx(k)) != 0)
{
nnzC++;
break;
}
}
}
}
}
*/
// 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);
}
else
{
col_index_increment = 1;
c = SparseMatrix(colsB, rowsA, colsB*rowsA);
}
// Start the matrix composition
// Declare variables for the T-Norm and S-Norm values
double snorm_val;
double tnorm_val;
// Initialize the number of nonzero elemnts in sparse matrix c
int nel = 0;
c.xcidx(0) = 0;
// Calculate the composition for the specified rows (between 0 and rowsA)
for (int i = 0; i < rowsA; i++)
{
for(int j = lock_option*i; j < colsB; j = j+col_index_increment)
{
snorm_val = calc_tnorm(constA(i,0), constB(0,j));
// From b.cidx(j) to b.cidx(j+1)-1 there are all the nonzero elements of the column j of b (constB actually)
for (int k = constB.cidx(j); k <= (constB.cidx(j+1)-1); k++)
{
// b.ridx(k) is the row index of the actual nonzero element in b (constB actually)
// so it is also the column index of the element in a (constA actually)
tnorm_val = calc_tnorm(constA(i, constB.ridx(k)), constB.data(k));
snorm_val = calc_snorm(snorm_val, tnorm_val);
}
if (snorm_val != 0)
{
// Equivalent to c(i, j) = snorm_val;
c.xridx(nel) = j;
c.xdata(nel++) = snorm_val;
}
}
c.xcidx(i+1) = nel;
}
// Compress the result sparse matrix because it is initialized with a number of nonzero element probably greater than the real one
c.maybe_compress();
octave_stdout << "OK\n";
// Return the result spars matrix (it's trasposed because we have been working on the rows of its transpose
return octave_value(c.transpose());
}
/* 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();
}
/* 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;
}
------------------------------------------------------------------------------
Gaining the trust of online customers is vital for the success of any company
that requires sensitive data to be transmitted over the Web. Learn how to
best implement a security strategy that keeps consumers' information secure
and instills the confidence they need to proceed with transactions.
http://p.sf.net/sfu/oracle-sfdevnl
_______________________________________________
Octave-dev mailing list
Octave-dev@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/octave-dev
