Dear PETSc Developers,
I am having trouble interfacing SuperLU_Dist as a direct solver for certain
problems in PETSc. The problem is that when interfacing with SuperLU_Dist,
you need your matrix to be of Type MPISEQAIJ when running MPI with one
processor. PETSc has long allowed the use of Matrix type MPIAIJ for all MPI
runs, including MPI with a single processor and that is still the case for
all of PETSc's native solvers. This however has been broken for the
SuperLU_Dist option. The following code snippet (in PETSc and not
SuperLU_Dist) is responsible for this restriction and I do not know if it
is by design or accident :
In file petsc-3.12.4/src/mat/impls/aij/mpi/superlu_dist/superlu_dist.c line
257 onwards :
ierr = MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);CHKERRQ(ierr);
if (size == 1) {
aa = (Mat_SeqAIJ*)A->data;
rstart = 0;
nz = aa->nz;
} else {
Mat_MPIAIJ *mat = (Mat_MPIAIJ*)A->data;
aa = (Mat_SeqAIJ*)(mat->A)->data;
bb = (Mat_SeqAIJ*)(mat->B)->data;
ai = aa->i; aj = aa->j;
bi = bb->i; bj = bb->j;
The code seems to check for number of processors and if it is = 1 conclude
that the matrix is a Mat_SeqAIJ and perform some operations. Only if
number-of-procs > 1 then it assumes that matrix is of type Mat_MPIAIJ. I
think this is unwieldy and lacks generality. One would like the same piece
of code to run in MPI mode for all processors with type Mat_MPIAIJ. Also
this restriction has suddenly appeared in a recent version. The issue was
not there until at least 3.9.4. So my question is from now (e.g. v12.4) on,
should we always use matrix type Mat_SeqAIJ when running on 1 processor
even with MPI enabled and use Mat_MPIAIJ for more than 1 processor. That is
use the number of processors in use as a criterion to set the matrix type ?
A an illustration, I have attached a minor modification of KSP example 12,
that used to work with all PETSc versions until at least 3.9.4 but now
throws a segmentation fault. It was compiled with MPI and run with mpiexec
-n 1 ./ex12
If I remove the "ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr);" it is okay.
I hope you can clarify my confusion.
Regards,
Rochan
static char help[] = "Solves a linear system in parallel with KSP.\n\
Input parameters include:\n\
-m <mesh_x> : number of mesh points in x-direction\n\
-n <mesh_n> : number of mesh points in y-direction\n\n";
/*T
Concepts: KSP^solving a system of linear equations
Concepts: KSP^Laplacian, 2d
Concepts: PC^registering preconditioners
Processors: n
T*/
/*
Demonstrates registering a new preconditioner (PC) type.
To register a PC type whose code is linked into the executable,
use PCRegister(). To register a PC type in a dynamic library use PCRegister()
Also provide the prototype for your PCCreate_XXX() function. In
this example we use the PETSc implementation of the Jacobi method,
PCCreate_Jacobi() just as an example.
See the file src/ksp/pc/impls/jacobi/jacobi.c for details on how to
write a new PC component.
See the manual page PCRegister() for details on how to register a method.
*/
/*
Include "petscksp.h" so that we can use KSP solvers. Note that this file
automatically includes:
petscsys.h - base PETSc routines petscvec.h - vectors
petscmat.h - matrices
petscis.h - index sets petscksp.h - Krylov subspace methods
petscviewer.h - viewers petscpc.h - preconditioners
*/
#include <petscksp.h>
PETSC_EXTERN PetscErrorCode PCCreate_Jacobi(PC);
int main(int argc,char **args)
{
Vec x,b,u; /* approx solution, RHS, exact solution */
Mat A; /* linear system matrix */
KSP ksp; /* linear solver context */
PetscReal norm; /* norm of solution error */
PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,its;
PetscErrorCode ierr;
PetscScalar v,one = 1.0;
PC pc; /* preconditioner context */
ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
ierr = PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the matrix and right-hand-side vector that define
the linear system, Ax = b.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create parallel matrix, specifying only its global dimensions.
When using MatCreate(), the matrix format can be specified at
runtime. Also, the parallel partitioning of the matrix can be
determined by PETSc at runtime.
*/
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
ierr = MatSetType(A,MATMPIAIJ);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
/*
Currently, all PETSc parallel matrix formats are partitioned by
contiguous chunks of rows across the processors. Determine which
rows of the matrix are locally owned.
*/
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
/*
Set matrix elements for the 2-D, five-point stencil in parallel.
- Each processor needs to insert only elements that it owns
locally (but any non-local elements will be sent to the
appropriate processor during matrix assembly).
- Always specify global rows and columns of matrix entries.
*/
for (Ii=Istart; Ii<Iend; Ii++) {
v = -1.0; i = Ii/n; j = Ii - i*n;
if (i>0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr);
}
/*
Assemble matrix, using the 2-step process:
MatAssemblyBegin(), MatAssemblyEnd()
Computations can be done while messages are in transition
by placing code between these two statements.
*/
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/*
Create parallel vectors.
- When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
we specify only the vector's global
dimension; the parallel partitioning is determined at runtime.
- When solving a linear system, the vectors and matrices MUST
be partitioned accordingly. PETSc automatically generates
appropriately partitioned matrices and vectors when MatCreate()
and VecCreate() are used with the same communicator.
- Note: We form 1 vector from scratch and then duplicate as needed.
*/
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
/*
Set exact solution; then compute right-hand-side vector.
Use an exact solution of a vector with all elements of 1.0;
*/
ierr = VecSet(u,one);CHKERRQ(ierr);
ierr = MatMult(A,u,b);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the linear solver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create linear solver context
*/
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
/*
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
/*
First register a new PC type with the command PCRegister()
*/
//ierr = PCRegister("ourjacobi",PCCreate_Jacobi);CHKERRQ(ierr);
/*
Set the PC type to be the new method
*/
ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
ierr = PCSetType(pc,PCLU);CHKERRQ(ierr);
ierr = PCFactorSetMatSolverType(pc, MATSOLVERSUPERLU_DIST);CHKERRQ(ierr);
/*
Set runtime options, e.g.,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
These options will override those specified above as long as
KSPSetFromOptions() is called _after_ any other customization
routines.
*/
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the linear system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Check solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Check the error
*/
ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
/*
Print convergence information. PetscPrintf() produces a single
print statement from all processes that share a communicator.
*/
ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr);
/*
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
*/
ierr = KSPDestroy(&ksp);CHKERRQ(ierr);
ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr);
/*
Always call PetscFinalize() before exiting a program. This routine
- finalizes the PETSc libraries as well as MPI
- provides summary and diagnostic information if certain runtime
options are chosen (e.g., -log_view).
*/
ierr = PetscFinalize();
return ierr;
}
