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

[S V N Vishwanathan]

> 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.  

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.

If we set A to be MPIDENSE then PETSc cannot compute A^T A. So we
thought about using Elemental. But now the issue we are running into is
that we cannot form A^T T1 if A is of type MATELEMENTAL. 

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. 

thanks

vishy