[petsc-users] Calculate only the diagonal elements for matrix matrix multiplications

[Barry Smith]

   Say you want to compute the diagonal of A * B where A and B are PETSc SeqAIJ matrices. If we use "MATLAB" notation A(i,:) to represent the ith row of A and B(:,j) to represent the nth column of B then you are computing  for each i A(i,:) * B(:,i) that is the ith row of A with the ith column of B.

   PETSc SeqAIJ matrices are stored "by row" and for each row is stored the column indices of the nonzeros (sorted) and the numerical values for those columns. If the nonzero structure of B is not related to the non-zero structure of A then computing the diagonal is a set of "sparse" vector vector products. Generic sparse vector vector products are never super efficient and to make matters worse with SeqAIJ we don't have stored directly the columns of B (only the rows) 

  If you look at MatMatMult_SeqAIJ_SeqAIJ() you will see we have SIX implementations using different approaches to try to get reasonable performance of the sparse matrix-matrix product. So if this computation of the diagonal is an important kernel in your code (and takes significant time) then you have some work cut out for you to write the appropriate code using the SeqAIJ data structure directly. We don't provide anything out of the box for this computation. I suggest looking at the various approaches for sparse matrix matrix product and see what one looks reasonable to specialize to compute just the diagonal portion.

  Barry


> On Jun 22, 2016, at 3:41 PM, ehsan sadrfaridpour <it.sadr at gmail.com> wrote:
> 
> Now, they are sparse (AIJ) and sequential. But I will upgrade the code to parallel later.
> 
> 
> On Wed, Jun 22, 2016 at 4:34 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
>    Matrix A and B, sparse or dense, parallel or sequential?
> 
> 
> > On Jun 22, 2016, at 3:27 PM, ehsan sadrfaridpour <it.sadr at gmail.com> wrote:
> >
> > Hi,
> > I need the diagonal elements from result of multiplying matrix A, B together.
> > I can get by multiplying the 2 matrices and then call MatGetDiagonal.
> > However, I only need the values on the diagonal and the rest of the elements are useless for me.
> > As the size of matrices increase, I am afraid it affect the performance.
> >
> > So, I am looking for another method that only calculates the diagonal of multiplication and not the rest of elements.
> > If there is not such a method for matrix, how can I use the vector's operations to reach the same results?
> >
> >
> > Best,
> > Ehsan
> 
>