# Ublas :

Subject: Re: [ublas] New uBLAS maintainer
From: Riccardo Rossi (rrossi_at_[hidden])
Date: 2010-07-27 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)

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
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];
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 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){
length++;
}
}
}

for(I jj = 0; jj < length; jj++){

nnz++;
}

next[temp] = -1; //clear arrays
sums[temp] = 0;
}

Cp[i+1] = nnz;
}
}

```--
On Fri, 2010-07-23 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 on-board 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
>                 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)
>
>                 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]
>
>
>
>
```