Indeed - I was just using the default solver (GMRES with ILU).
Using just standard LU (direct solve with "-pc_type lu -ksp_type
preonly"), I find elemental to be extremely slow even for a 1000x1000
matrix. For MPIaij it's throwing me an error if I tried "-pc_type lu".
I'm attaching the code here, in case you'd like to have a look at what
I've been trying to do.
The two configurations of interest are,
$> mpirun -n 4 ./ksps -N 1000 -mat_type mpiaij
$> mpirun -n 4 ./ksps -N 1000 -mat_type elemental
(for the GMRES with ILU) and,
$> mpirun -n 4 ./ksps -N 1000 -mat_type mpiaij -pc_type lu -ksp_type
preonly
$> mpirun -n 4 ./ksps -N 1000 -mat_type elemental -pc_type lu
-ksp_type preonly
elemental seems to perform poorly in both cases.
Nidish
On 8/7/20 12:50 AM, Barry Smith wrote:
What is the output of -ksp_view for the two case?
It is not only the matrix format but also the matrix solver that
matters. For example if you are using an iterative solver the
elemental format won't be faster, you should use the PETSc MPIDENSE
format. The elemental format is really intended when you use a direct
LU solver for the matrix. For tiny matrices like this an iterative
solver could easily be faster than the direct solver, it depends on
the conditioning (eigenstructure) of the dense matrix. Also the
default PETSc solver uses block Jacobi with ILU on each process if
using a sparse format, ILU applied to a dense matrix is actually LU so
your solver is probably different also between the MPIAIJ and the
elemental.
Barry
On Aug 7, 2020, at 12:30 AM, Nidish <[email protected]
<mailto:[email protected]>> wrote:
Thank you for the response.
I've just been running some tests with matrices up to 2e4 dimensions
(dense). When I compared the solution times for "-mat_type elemental"
and "-mat_type mpiaij" running with 4 cores, I found the mpidense
versions running way faster than elemental. I have not been able to
make the elemental version finish up for 2e4 so far (my patience runs
out faster).
What's going on here? I thought elemental was supposed to be superior
for dense matrices.
I can share the code if that's appropriate for this forum (sorry, I'm
new here).
Nidish
On Aug 6, 2020, at 23:01, Barry Smith <[email protected]
<mailto:[email protected]>> wrote:
On Aug 6, 2020, at 7:32 PM, Nidish <[email protected]
<mailto:[email protected]>> wrote: I'm relatively new to PETSc,
and my applications involve (for the most part) dense matrix
solves. I read in the documentation that this is an area
PETSc does not specialize in but instead recommends external
libraries such as Elemental. I'm wondering if there are any
"best" practices in this regard. Some questions I'd like
answered are: 1. Can I just declare my dense matrix as a
sparse one and fill the whole matrix up? Do any of the others
go this route? What're possible pitfalls/unfavorable outcomes
for this? I understand the memory overhead probably shoots up.
No, this isn't practical, the performance will be terrible.
2. Are there any specific guidelines on when I can expect
elemental to perform better in parallel than in serial?
Because the computation to communication ratio for dense matrices is
higher than for sparse you will see better parallel performance for dense
problems of a given size than sparse problems of a similar size. In other words
parallelism can help for dense matrices for relatively small problems, of
course the specifics of your machine hardware and software also play a role.
Barry
Of course, I'm interesting in any other details that may be
important in this regard. Thank you, Nidish
--
Nidish
char help[] = "Testing out KSP\n"
"second line\n";
#include <iostream>
#include <petscksp.h>
#include <petscpc.h>
int main(int nargs, char *sargs[])
{
PetscMPIInt rank, size;
PetscInt ierr;
PetscInt N=10;
ierr = PetscInitialize(&nargs, &sargs, (char*)0, help); if (ierr) return ierr;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
MPI_Comm_size(PETSC_COMM_WORLD, &size);
PetscOptionsGetInt(NULL, NULL, "-N", &N, NULL);
Mat A;
Vec x,b;
PetscInt m, n;
MatCreate(PETSC_COMM_WORLD, &A);
MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, N, N);
// MatSetType(A, MATELEMENTAL);
MatSetFromOptions(A);
// MatMPIAIJSetPreallocation(A, N/size, NULL, N*(size-1)/size, NULL);
MatSetUp(A);
MatGetLocalSize(A, &m, &n);
// PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d]=(%d, %d)\n", rank, m, n);
// PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
MatType mtp;
MatGetType(A, &mtp);
if (!strcmp(mtp, MATMPIAIJ)) { // Allocate if mpiaij
MatMPIAIJSetPreallocation(A, m, NULL, N-m, NULL);
}
PetscRandom pr=NULL;
PetscRandomCreate(PETSC_COMM_WORLD, &pr);
PetscInt lo, hi;
PetscScalar tmdat=0.0;
MatGetOwnershipRange(A, &lo, &hi);
// MatSetRandom(A, pr);
for (int i=lo; i<hi; ++i) {
for (int j=0; j<N; ++j) {
// PetscRandomGetValue(pr, &tmdat);
MatSetValue(A, i, j, tmdat, ADD_VALUES);
}
MatSetValue(A, i, i, 10.0*pow(-1, i), ADD_VALUES);
}
// PetscRandomGetValue(pr, PetscScalar *);
MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
MatCreateVecs(A, &x, &b);
VecSetFromOptions(x);
VecSetFromOptions(b);
VecGetOwnershipRange(b, &lo, &hi);
for (int i=lo; i<hi; ++i)
VecSetValue(b, i, (PetscScalar)10.0, ADD_VALUES);
VecAssemblyBegin(b);
VecAssemblyEnd(b);
KSP ksp;
KSPCreate(PETSC_COMM_WORLD, &ksp);
KSPSetOperators(ksp, A, A);
ierr = KSPSetTolerances(ksp,1.e-2/((N+1)*(N+1)),1.e-50,PETSC_DEFAULT,
PETSC_DEFAULT);CHKERRQ(ierr);
KSPSetFromOptions(ksp);
PC pc;
KSPGetPC(ksp, &pc);
PCSetFromOptions(pc);
// PCSetType(pc, PCLU);
// PCSetType(pc, PCCHOLESKY);
KSPSolve(ksp, b, x);
VecAssemblyBegin(x);
VecAssemblyEnd(x);
Vec err;
VecDuplicate(b, &err);
MatMult(A, x, err);
VecAXPY(err, -1.0, b);
PetscScalar errnorm;
VecNorm(err, NORM_2, &errnorm);
PetscPrintf(PETSC_COMM_WORLD, "err = %e\n", errnorm);
// VecView(b, PETSC_VIEWER_STDOUT_WORLD);
// PetscPrintf(PETSC_COMM_WORLD, "-------------------------\n");
// PetscPrintf(PETSC_COMM_WORLD, "Solution\n");
// VecView(x, PETSC_VIEWER_STDOUT_WORLD);
// PetscPrintf(PETSC_COMM_WORLD, "-------------------------\n");
// MatMult(A, x, b);
// PetscPrintf(PETSC_COMM_WORLD, "Ax=\n");
// VecView(b, PETSC_VIEWER_STDOUT_WORLD);
// PetscPrintf(PETSC_COMM_WORLD, "-------------------------\n");
// PetscPrintf(PETSC_COMM_WORLD, "Done\n");
MatDestroy(&A);
VecDestroy(&x);
VecDestroy(&b);
KSPDestroy(&ksp);
PetscFinalize();
return 0;
}
