
Ublas : 
Subject: Re: [ublas] New uBLAS maintainer
From: Riccardo Rossi (rrossi_at_[hidden])
Date: 20100727 03:16:36
Dear David,
i was digging in our code but i found only a specialized version
of the matmat routine to a case we use often (A*diag*B)
in any case this two routines, taken from scipy.sparse, do the matrix
matrix product in a very efficient way.
The idea is that in a first pass they compute the graph and then they
fill it (unordered if i remember correctly)
hope this can be helpful
ciao
Riccardo
/*
* Pass 1 computes CSR row pointer for the matrix product C = A * B
*
*/
template <class I>
void csr_matmat_pass1(const I n_row,
const I n_col,
const I Ap[],
const I Aj[],
const I Bp[],
const I Bj[],
I Cp[])
{
// method that uses O(n) temp storage
std::vector<I> mask(n_col, 1);
Cp[0] = 0;
I nnz = 0;
for(I i = 0; i < n_row; i++){
for(I jj = Ap[i]; jj < Ap[i+1]; jj++){
I j = Aj[jj];
for(I kk = Bp[j]; kk < Bp[j+1]; kk++){
I k = Bj[kk];
if(mask[k] != i){
mask[k] = i;
nnz++;
}
}
}
Cp[i+1] = nnz;
}
}
/*
* Pass 2 computes CSR entries for matrix C = A*B using the
* row pointer Cp[] computed in Pass 1.
*
*/
template <class I, class T>
void csr_matmat_pass2(const I n_row,
const I n_col,
const I Ap[],
const I Aj[],
const T Ax[],
const I Bp[],
const I Bj[],
const T Bx[],
I Cp[],
I Cj[],
T Cx[])
{
std::vector<I> next(n_col,1);
std::vector<T> sums(n_col, 0);
I nnz = 0;
Cp[0] = 0;
for(I i = 0; i < n_row; i++){
I head = 2;
I length = 0;
I jj_start = Ap[i];
I jj_end = Ap[i+1];
for(I jj = jj_start; jj < jj_end; jj++){
I j = Aj[jj];
T v = Ax[jj];
I kk_start = Bp[j];
I kk_end = Bp[j+1];
for(I kk = kk_start; kk < kk_end; kk++){
I k = Bj[kk];
sums[k] += v*Bx[kk];
if(next[k] == 1){
next[k] = head;
head = k;
length++;
}
}
}
for(I jj = 0; jj < length; jj++){
if(sums[head] != 0){
Cj[nnz] = head;
Cx[nnz] = sums[head];
nnz++;
}
I temp = head;
head = next[head];
next[temp] = 1; //clear arrays
sums[temp] = 0;
}
Cp[i+1] = nnz;
}
}
 On Fri, 20100723 at 14:44 +0200, David Bellot wrote: > Hi Riccardo, > > sorry for answer so late to you email. As it's now time for uBLAS to > get a new lifting, thanks to the new coming version of Boost, I was > reading old emails to gather all improvements. > I just realized you were saying that you have an efficient > sparse_matrix product. Well, of course we are interested. More than > you can imagine. > If you can contribute code, that will be most welcome. > > The current stage of uBLAS is the following: version 1.44 will come > with a new (beta) documentation in Doxygen, that's the big part of it > and a new improvements. It will be released next week but meanwhile it > is accessible from ublas.sourceforge.net for the documentation and > svn.boost.org for everything else (the documentation is also included > in the source code). > > Don't hesitate to send me code and contributions. > Many thanks and regards, > David Bellot > >  > David Bellot, PhD > david.bellot_at_[hidden] > http://david.bellot.free.fr > > > On Mon, Mar 15, 2010 at 15:59, <rrossi_at_[hidden]> wrote: > Hi David, > > i just wanted to say that if you find it useful we > have a sparse_matrix * sparse_matrix version that should be > quite efficient. > > > > unfortunately it is not really in ublas format, as the ublas > does not allow constructing a CSR matrix using the vectors of > indices and the vector of data. Let me know if you need it. > > > > We also have some CUDA solver, which may be of interest. Once > again it uses directly the three vectors of CSR, even if in > our code we internally use ublas ... > > > > let me know if some of this could be of help > > greetings > > Riccardo > > > > > > On Sun 14/03/10 9:44 PM , David Bellot > <david.bellot_at_[hidden]> wrote: > > Dear uBLAS users, > > recently Gunter Winkler asked for someone to take over > the maintenance of the uBLAS library. As a fervent > user of uBLAS and a strong believer in Open Source and > Free Software, I decided to propose myself as the new > duty man. I will therefore have the honor to be in > charge of uBLAS and will do my best to make uBLAS > reach its goals of versatility, performance and ease > of use. > > Let me quickly introduce myself: my name is David > Bellot, I hold a PhD in Computer Science and do > research (currently in the finance world) on Machine > Learning and probabilistic artificial intelligence, > hence my strong interest in uBLAS. > > First of all, I want to say a big Thank You to Gunter > for all the great work he did on uBLAS and a personal > Thank You to help me getting onboard and starting my > new duty as uBLAS maintainer. > > I would also like to talk about potential projects for > the next versions of uBLAS. For that purpose, I hope > everybody will be interested in bringing new ideas, > wishes and even new pieces of code: > > (1) as you can imagine, in machine learning, one often > needs to "randomly" access to submatrices. A good > framework is already in place for matrix_view, I would > like to extend it so that to make it as versatile as > it is in other libraries or even Matlab. > > (2) after reading last week emails, I think we could > provide basic implementations of a few standard > algorithms like inversion, solvers, etc... > > (3) bindings are a hot topic. Let's be pragmatic: it's > not supposed to be part of uBLAS but having a standard > interface would add a strong value to uBLAS. And, I am > like others, I want to play with my brand new nVidia > card. > > (4) another hot topic which is a recurrent complain > about uBLAS: the product of 2 matrices. Do we want > prod(A,B) or A*B. Let's think about it because other > libraries implemented A*B in a very efficient manner > too. > > (5) Bindings for big libraries are also important and > subject to discussion. I think we have to work more on > the interface between all standard libraries when it > is needed because, at the end of the day, people also > want to use uBLAS to make computations with existing > standard and not just write brand new algorithms. > > (6) I will join Gunter in his effort to provide new > documentation, covering more topics, with tutorial and > advanced topics. uBLAS is a great library and a good > documentation is of primary interest. That is one of > the most important topic for me (yes, way more than > prod(A,B) versus A*B) > > Please everybody contribute with your own ideas and > desiderata. > Let's work and make uBLAS simply the best. > > With my Best Regards, > David Bellot > >  > David Bellot, PhD > david.bellot_at_[hidden] > http://david.bellot.free.fr > > > > > > > _______________________________________________ > ublas mailing list > ublas_at_[hidden] > http://lists.boost.org/mailman/listinfo.cgi/ublas > Sent to: david.bellot_at_[hidden] > > > >