[petsc-users] [dev at libelemental.org] Re: The product of two MPIDENSE matrices

[Jed Brown]

S V N Vishwanathan <vishy at stat.purdue.edu> writes:

>> You can use MatTransposeMatMult with an MPIAIJ and a MPIDENSE.  (This
>> is the standard case mentioned earlier, and I don't understand what
>> Elemental is offering you in this situation.  Do you have other dense
>> matrices running around?)
>
> Perhaps I should give a bit of background on what we are trying to
> achieve and the experts in the lists can give us some advice.
>
> We are trying to find a k rank approximation to a n x m x p tensor
> T. Here, n, m, p are of the order of millions but k is very small,
> typically somewhere between 10 and 100. Also T is extremely sparse (nnz
> is typically around 500 million or so).
>
> Let A (size n x k), B(m x k) and C (p x k) be the factors. 
>
> We are flattening the tensor along each of the three dimensions and
> representing them as three matrices T1 (size n x mp), T2 (size m x np)
> and T3 (size p x mn), all in MPIAIJ format.  

Thanks for this description.

> One of the intermediate computations that we perform is form the k x k
> matrix A^T A. Another computation that we perform is to compute A^T T1.

I'm in favor of implementing the (very simple)
MatTransposeMatMult_MPIDense_MPIDense.  You would usually want the
result distributed redundantly, but PETSc doesn't have a dedicated
matrix format for redundant matrices.  We can add that (it would also be
useful in other places), but MatTransposeMatMult doesn't currently state
the desired output distribution and I want to make sure the bases are
covered if we do that.  In the mean time, you can use MatDenseGetArray,
dgemm, and MPI_Allreduce because that's what
MatTransposeMatMult_MPIDense_MPIDense would be doing when the result is
to be redundantly distributed.

> Also somewhere in our code, we form a k x k matrix, collect it in one of
> the processors and invert it and broadcast it back to all the other
> processors. We need to use Cholesky or LU factorization for this
> calculation. 

You might be better off forming the matrix redundantly (allreduce if
needed) and factoring it redundantly.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 835 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20131001/32365c2c/attachment.pgp>