[petsc-users] MatMatMult takes much time

[Alexander Grayver]

On 12.03.2011 15:54, Barry Smith wrote:
>     A*B where A is sparse and B dense can be implemented efficiently since it is a bunch of sparse matrix vector products (one for each column of B)
>
>     S*G where S is dense is unlikley to ever be particularly fast: Here is how it is currently implemented
>
> #undef __FUNCT__
> #define __FUNCT__ "MatMatMultNumeric_MPIDense_MPIAIJ"
> /*
>      Computes (B'*A')' since computing B*A directly is untenable

Ok, it explains why the code:
               call MatTranspose(G,MAT_INITIAL_MATRIX,GT,ierr);CHKERRQ(ierr)
               call MatTranspose(S,MAT_INITIAL_MATRIX,ST,ierr);CHKERRQ(ierr)
               call 
MatMatMult(GT,ST,MAT_INITIAL_MATRIX,PETSC_DEFAULT_DOUBLE_PRECISION,JT,ierr);CHKERRQ(ierr)

works faster than just
               call 
MatMatMult(S,G,MAT_INITIAL_MATRIX,PETSC_DEFAULT_DOUBLE_PRECISION,J,ierr);CHKERRQ(ierr)


>                 n                       p                          p
>          (              )       (              )         (                  )
>        m (      A       )  *  n (       B      )   =   m (         C        )
>          (              )       (              )         (                  )
>
> */
> PetscErrorCode MatMatMultNumeric_MPIDense_MPIAIJ(Mat A,Mat B,Mat C)
> {
>    PetscErrorCode     ierr;
>    Mat                At,Bt,Ct;
>
>    PetscFunctionBegin;
>    ierr = MatTranspose(A,MAT_INITIAL_MATRIX,&At);CHKERRQ(ierr);
>    ierr = MatTranspose(B,MAT_INITIAL_MATRIX,&Bt);CHKERRQ(ierr);
>    ierr = MatMatMult(Bt,At,MAT_INITIAL_MATRIX,1.0,&Ct);CHKERRQ(ierr);
>    ierr = MatDestroy(At);CHKERRQ(ierr);
>    ierr = MatDestroy(Bt);CHKERRQ(ierr);
>    ierr = MatTranspose(Ct,MAT_REUSE_MATRIX,&C);CHKERRQ(ierr);
>    ierr = MatDestroy(Ct);CHKERRQ(ierr);
>    PetscFunctionReturn(0);
> }
>
> If you can suggest a way to implement it much more efficiently please let us know.
>
> Are you sure you need this product and not A*B?

Unfortunately I can't suggest a better way. Well actually I need A*B, 
but my initial matrices have to be transposed for that.
But probably there is a way to avoid forming dense matrix at all, this 
way requires more code, that is why I tried another way with full dense 
matrix first.

Thanks Barry.

>
>    Barry
>
>
> On Mar 12, 2011, at 8:17 AM, Alexander Grayver wrote:
>
>> Hello,
>>
>> call MatMatMult(S,G,MAT_INITIAL_MATRIX,PETSC_DEFAULT_DOUBLE_PRECISION,J,ierr);CHKERRQ(ierr)
>>
>> I use this operation to compute matrix-matrix product. It's correct, but quite slow.
>>
>> For instance, my matrix S is a dense matrix size of<6x762048 complex double>  and matrix G is a sparse matrix with 2 nonzeros per row and size of<762048x254016 complex double>
>> Time spent on S*G was ~22 secs with 32 processes over 4 nodes. There is infiniband between nodes.
>>
>> Is it a real life or I'm doing something wrong?
>>
>> Regards,
>> Alexander