[petsc-users] MatPtAP for involving MPIDENSE and MPIAIJ matrices

[Bikash Kanungo]

HI,

I'm trying to evaluate B = P^T * A * P, where P is an MPIDENSE matrix and A
is an MPIAIJ matrix. Neither MatTransposeMatMult nor MatPtAP support such
an operation for the given type of matrices. I tried to look into
MATELEMENTAL class and didn't find any MatTransposeMatMult function for a
pair of Elemental matrices.

What would be a possible and efficient way to perform such the above
operation?

Thanks,
Bikash

[Matthew Knepley]

Post by Bikash Kanungo
HI,
I'm trying to evaluate B = P^T * A * P, where P is an MPIDENSE matrix and
A is an MPIAIJ matrix. Neither MatTransposeMatMult nor MatPtAP support such
an operation for the given type of matrices. I tried to look into
MATELEMENTAL class and didn't find any MatTransposeMatMult function for a
pair of Elemental matrices.
What would be a possible and efficient way to perform such the above
operation?

Jed, is TAIJ the way to do this?

Matt

Post by Bikash Kanungo
Thanks,
Bikash
--
Bikash S. Kanungo
PhD Student
Computational Materials Physics Group
Mechanical Engineering
University of Michigan

[Barry Smith]

HI,
I'm trying to evaluate B = P^T * A * P, where P is an MPIDENSE matrix and A is an MPIAIJ matrix. Neither MatTransposeMatMult nor MatPtAP support such an operation for the given type of matrices. I tried to look into MATELEMENTAL class and didn't find any MatTransposeMatMult function for a pair of Elemental matrices.
What would be a possible and efficient way to perform such the above operation?

Hong will see about providing the missing piece. What are you sizes of the matrix P? Presumably it has many rows but few columns?

Barry

Jed, is TAIJ the way to do this?
Matt
Thanks,
Bikash
--
Bikash S. Kanungo
PhD Student
Computational Materials Physics Group
Mechanical Engineering
University of Michigan
--
What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
-- Norbert Wiener

[Jed Brown]

Post by Matthew Knepley
Jed, is TAIJ the way to do this?

No, TAIJ is unrelated.

[Matthew Knepley]

Post by Jed Brown

Post by Matthew Knepley
Jed, is TAIJ the way to do this?

No, TAIJ is unrelated.

I do not understand "unrelated". My understanding was that TAIJ could be
used (with T = I)
to get the action of A on a set of vectors, which I think would be A P. Why
specifically would
you not use it?

Thanks,

Matt

[Jed Brown]

Post by Matthew Knepley
I do not understand "unrelated". My understanding was that TAIJ could be
used (with T = I)
to get the action of A on a set of vectors, which I think would be A P. Why
specifically would
you not use it?

MAIJ is enough for that; it would involve packing P into a Vec, but
Bikash might as well use MatMatMult to perform the same operation in a
more direct way (MatMatMult_MPIAIJ_MPIDense exists). The problem is
that MatTransposeMatMult_MPIDense_MPIDense needs to be implemented.

[Hong]

I plan to implement MatTransposeMatMult_MPIDense_MPIDense via

1. add MatTransposeMatMult_elemental_elemental()
2. C_dense = P_dense^T * B_dense
via MatConvert_dense_elemental() and MatConvert_elemental_dense()

Let me know if you have better suggestions.

Hong

Post by Jed Brown

Post by Matthew Knepley
I do not understand "unrelated". My understanding was that TAIJ could be
used (with T = I)
to get the action of A on a set of vectors, which I think would be A P.

Why

Post by Matthew Knepley
specifically would
you not use it?

MAIJ is enough for that; it would involve packing P into a Vec, but
Bikash might as well use MatMatMult to perform the same operation in a
more direct way (MatMatMult_MPIAIJ_MPIDense exists). The problem is
that MatTransposeMatMult_MPIDense_MPIDense needs to be implemented.

[Jed Brown]

r<#secure method=pgpmime mode=sign>

Post by Hong
I plan to implement MatTransposeMatMult_MPIDense_MPIDense via
1. add MatTransposeMatMult_elemental_elemental()
2. C_dense = P_dense^T * B_dense
via MatConvert_dense_elemental() and MatConvert_elemental_dense()

The above involves a ton of data movement and MPIDense is a logical
distribution for matrices with a modest number of columns. I think I
would just do the local GEMM and then MPI_Reduce_scatter it.

[Hong]

Post by Jed Brown

Post by Hong
I plan to implement MatTransposeMatMult_MPIDense_MPIDense via
1. add MatTransposeMatMult_elemental_elemental()
2. C_dense = P_dense^T * B_dense
via MatConvert_dense_elemental() and MatConvert_elemental_dense()

The above involves a ton of data movement and MPIDense is a logical
distribution for matrices with a modest number of columns. I think I
would just do the local GEMM and then MPI_Reduce_scatter it.

This would work.

However, I recall that you did a smart ordering which allows
MatConvert_mpidense_elemental() uses same physical matrix storage for petsc
and elemental, but logically in the layout of elemental. As an example,
petsc/src/mat/examples/tests/ex103.c:
mpiexec -n 6 ./ex103
Outplace MatConvert, A_elemental:
Mat Object: 6 MPI processes
type: elemental
Elemental matrix (cyclic ordering)
0 0 0 0 0
1 1 1 1 1
2 2 2 2 2
3 3 3 3 3
4 4 4 4 4
5 5 5 5 5
0 0 0 0 0
1 1 1 1 1
2 2 2 2 2
3 3 3 3 3

Elemental matrix (explicit ordering)
Mat Object: 6 MPI processes
type: mpidense
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
0.0000000000000000e+00 0.0000000000000000e+00
1.0000000000000000e+00 1.0000000000000000e+00 1.0000000000000000e+00
1.0000000000000000e+00 1.0000000000000000e+00
1.0000000000000000e+00 1.0000000000000000e+00 1.0000000000000000e+00
1.0000000000000000e+00 1.0000000000000000e+00
2.0000000000000000e+00 2.0000000000000000e+00 2.0000000000000000e+00
2.0000000000000000e+00 2.0000000000000000e+00
2.0000000000000000e+00 2.0000000000000000e+00 2.0000000000000000e+00
2.0000000000000000e+00 2.0000000000000000e+00
3.0000000000000000e+00 3.0000000000000000e+00 3.0000000000000000e+00
3.0000000000000000e+00 3.0000000000000000e+00
3.0000000000000000e+00 3.0000000000000000e+00 3.0000000000000000e+00
3.0000000000000000e+00 3.0000000000000000e+00
4.0000000000000000e+00 4.0000000000000000e+00 4.0000000000000000e+00
4.0000000000000000e+00 4.0000000000000000e+00
5.0000000000000000e+00 5.0000000000000000e+00 5.0000000000000000e+00
5.0000000000000000e+00 5.0000000000000000e+00

i.e., elemental and petsc dense matrices have same ownership.
If there is no data movement for MatConvert(), then it would be easier to
use elemental.

Hong

[Jed Brown]

Post by Hong

Post by Jed Brown

Post by Hong
I plan to implement MatTransposeMatMult_MPIDense_MPIDense via
1. add MatTransposeMatMult_elemental_elemental()
2. C_dense = P_dense^T * B_dense
via MatConvert_dense_elemental() and MatConvert_elemental_dense()

The above involves a ton of data movement and MPIDense is a logical
distribution for matrices with a modest number of columns. I think I
would just do the local GEMM and then MPI_Reduce_scatter it.

This would work.
However, I recall that you did a smart ordering which allows
MatConvert_mpidense_elemental() uses same physical matrix storage for petsc
and elemental,

Same storage for vectors. This is your code, but I think you'll find
that it moves the matrix entries around. (Note that Elemental [MC,MR]
is a 2D distribution while MPIDense is 1D.) Also, I think it would be
better if this particular operation did not depend on Elemental.

You could write a conversion to Elemental [VC,*], which would then match
the MPIDense distribution and thus not need to move any matrix entries.

[Hong]

Bikash,
I implemented MatTransposeMatMult_MPIDense_MPIDense()
in the branch hzhang/mattransmatmult_dense.
Once it passes our nightly tests, I'll merge it to petsc master branch.

See petsc/src/mat/examples/tests/ex104.c, in which I added a simple
example doing
B = P^T * A * P, where P is an MPIDENSE matrix and A is an MPIAIJ matrix.

Let us know if you see any bug or performance issues.

Hong

Post by Jed Brown

Post by Hong

Post by Jed Brown

Post by Hong
I plan to implement MatTransposeMatMult_MPIDense_MPIDense via
1. add MatTransposeMatMult_elemental_elemental()
2. C_dense = P_dense^T * B_dense
via MatConvert_dense_elemental() and MatConvert_elemental_dense()

The above involves a ton of data movement and MPIDense is a logical
distribution for matrices with a modest number of columns. I think I
would just do the local GEMM and then MPI_Reduce_scatter it.

This would work.
However, I recall that you did a smart ordering which allows
MatConvert_mpidense_elemental() uses same physical matrix storage for

petsc

Post by Hong
and elemental,

Same storage for vectors. This is your code, but I think you'll find
that it moves the matrix entries around. (Note that Elemental [MC,MR]
is a 2D distribution while MPIDense is 1D.) Also, I think it would be
better if this particular operation did not depend on Elemental.
You could write a conversion to Elemental [VC,*], which would then match
the MPIDense distribution and thus not need to move any matrix entries.

[Bikash Kanungo]

Thanks a ton, Hong. You guys are amazing! I'll report any bug or
performance issue if I find any.

Post by Hong
Bikash,
I implemented MatTransposeMatMult_MPIDense_MPIDense()
in the branch hzhang/mattransmatmult_dense.
Once it passes our nightly tests, I'll merge it to petsc master branch.
See petsc/src/mat/examples/tests/ex104.c, in which I added a simple
example doing
B = P^T * A * P, where P is an MPIDENSE matrix and A is an MPIAIJ matrix.
Let us know if you see any bug or performance issues.
Hong

Post by Jed Brown

Post by Hong

Post by Jed Brown

Post by Hong
I plan to implement MatTransposeMatMult_MPIDense_MPIDense via
1. add MatTransposeMatMult_elemental_elemental()
2. C_dense = P_dense^T * B_dense
via MatConvert_dense_elemental() and MatConvert_elemental_dense()

The above involves a ton of data movement and MPIDense is a logical
distribution for matrices with a modest number of columns. I think I
would just do the local GEMM and then MPI_Reduce_scatter it.

This would work.
However, I recall that you did a smart ordering which allows
MatConvert_mpidense_elemental() uses same physical matrix storage for

petsc

Post by Hong
and elemental,

Same storage for vectors. This is your code, but I think you'll find
that it moves the matrix entries around. (Note that Elemental [MC,MR]
is a 2D distribution while MPIDense is 1D.) Also, I think it would be
better if this particular operation did not depend on Elemental.
You could write a conversion to Elemental [VC,*], which would then match
the MPIDense distribution and thus not need to move any matrix entries.