Boost logo

Ublas :

Date view Thread view Subject view Author view

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)

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, 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
>                 needs to "randomly" access to sub-matrices. 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]
> 
> 
> 
>
Date view Thread view Subject view Author view