To whom it may concern, I would like to modify the code in SNES's ex2.c to run in parallel.
https://www.mcs.anl.gov/petsc/petsc-3.6.4/src/snes/examples/tutorials/ex2.c.html I modified my ex2.c, as attached here. Even though ex3.c is the parallel version, I would like my implementation to "not contain" DMDA. (Despite version 3.6.4, it works for version 3.15.0, which I use in my debugging purpose.) The "crucial" parts of ex2.c that I modified are (apart from deletion of comments, exact solution verification, etc): - (line 63-64) matrix preallocation - use MatGetOwnershipRange on a PETSc matrix jac (in FormJacobian) - use VecGetOwnershipRange on every PETSc vector in the program I try to run the program with mpiexec. For n = 1 it works fine for -pc_type lu -snes_type newtonls. However, when I increase n to be n = 2, the code does not iterate further. (It does not go into the FormJacobian for after the first iteration.) I try several different pc_type and snes_type, but still do not get the result as obtained with n = 1. I attached the terminal printout also. Since the manual (3.15.0) said that (p81, section 2.4, 1st paragraph) "Also, the SNES interface is identical for the uniprocess and parallel cases; the only difference in the parallel version is that each process typically forms only its local contribution to various matrices and vectors" Therefore, I think that by modified only the vectors and matrices part, I should be able to get parallelism. Regards, Suwun PS1: I compile using mpicc and use mpiexec found in the directory ${PETSC_DIR}/arch-mumps-opt/bin/ to run the code. PS2: I attached the install.sh (I install with MUMPS) that I use to install PETSc. Note that I also turn on --with-debuggin=yes PS3: my OS is Ubuntu 20.04.2. Intel CPU. My environment should have everything necessary (gfortran, gcc, g++, valgrind, cmake).
ssuwunnarat@ssu-u2004-laptop:~/Documents/Inverse_Design_Drive/test_code/test_makesnes_parallel$ ~/soft/petsc-3.15.0/arch-mumps-opt/bin/mpiexec -n 1 main_exec -pc_type lu -snes_type newtonls ---------------- no. of MPI comm size = 1 ----------------- (rstart_r, rend_r) = (0,10) (rstart_F, rend_F) = (0,10) ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,10) vecnorm(f) = 8.0664075e+00 iter = 0, SNES Function norm 8.06641 [0]PETSC ERROR: PETSc installed without X windows or Microsoft Graphics on this machine proceeding without graphics ----------------- FormJacobian Called -------------------- (rstart, rend) = (0,10) matnorm = 5.5703709e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,10) vecnorm(f) = 1.9286608e+00 iter = 1, SNES Function norm 1.92866 ----------------- FormJacobian Called -------------------- (rstart, rend) = (0,10) matnorm = 5.6044198e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,10) vecnorm(f) = 1.3201094e-02 iter = 2, SNES Function norm 0.0132011 ----------------- FormJacobian Called -------------------- (rstart, rend) = (0,10) matnorm = 5.6016400e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,10) vecnorm(f) = 7.6723813e-07 iter = 3, SNES Function norm 7.67238e-07 ----------------- FormJacobian Called -------------------- (rstart, rend) = (0,10) matnorm = 5.6016191e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,10) vecnorm(f) = 1.0771452e-14 iter = 4, SNES Function norm 1.07715e-14 number of SNES iterations = 4 ssuwunnarat@ssu-u2004-laptop:~/Documents/Inverse_Design_Drive/test_code/test_makesnes_parallel$ ~/soft/petsc-3.15.0/arch-mumps-opt/bin/mpiexec -n 2 main_exec -pc_type lu -snes_type newtonls ---------------- no. of MPI comm size = 2 ----------------- (rstart_r, rend_r) = (0,5) (rstart_F, rend_F) = (0,5) ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 iter = 0, SNES Function norm 104.44 [0]PETSC ERROR: PETSc installed without X windows or Microsoft Graphics on this machine proceeding without graphics [1]PETSC ERROR: PETSc installed without X windows or Microsoft Graphics on this machine proceeding without graphics ----------------- FormJacobian Called -------------------- (rstart, rend) = (0,5) matnorm = 5.5911552e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.5752076e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.1580581e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0740469e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0535850e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0474159e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0454111e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0447401e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0445133e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0444363e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0444102e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0444013e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443983e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443973e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443969e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443968e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443968e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 ----------------- FormFunction Called -------------------- (rstart_f, rend_f) = (0,5) vecnorm(f) = 1.0443967e+02 number of SNES iterations = 0 ssuwunnarat@ssu-u2004-laptop:~/Documents/Inverse_Design_Drive/test_code/test_makesnes_parallel$
#include <petscsnes.h>
#include <stdio.h>
#include <math.h>
#include <nlopt.h>
extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
extern PetscErrorCode FormInitialGuess(Vec);
extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void*);
extern PetscErrorCode PrintVecToFile_r_i(Vec,char*,char*,int);
// Solve x''(t) + x^2 (t) = f(t) with f(t) given as f(t) = 6*t + t^6
// expect our solution (in the x_e_real_vec.qxz file) as x(t) = t^3
typedef struct {
PetscViewer viewer;
} MonitorCtx;
int main(int argc,char **argv)
{
SNES snes;
KSP ksp;
PC pc;
Vec x,r,F;
Mat J;
MonitorCtx monP;
PetscErrorCode ierr;
PetscInt its,n = 10,i,maxit,maxf,rstart_F,rend_F,rstart_r,rend_r;
PetscMPIInt size;
PetscScalar h,xp,v,r_val;
PetscReal abstol,rtol,stol;
PetscInitialize(&argc,&argv,(char*)0,NULL);
ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);CHKERRQ(ierr);
h = 1.0/(n-1);
PetscPrintf(PETSC_COMM_WORLD,"\n ---------------- no. of MPI comm size = %i ----------------- \n", size);fflush(stdout);
ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr); // Parallel vector
ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);
ierr = VecSetFromOptions(x);CHKERRQ(ierr);
ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
ierr = VecDuplicate(x,&F);CHKERRQ(ierr);
ierr = VecGetOwnershipRange(r,&rstart_r,&rend_r);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"\n (rstart_r, rend_r) = (%i,%i) \n", rstart_r, rend_r);fflush(stdout);
r_val = 0.0;
for (i=rstart_r; i<rend_r; i++) {
ierr = VecSetValues(r,1,&i,&r_val,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = VecAssemblyBegin(r); CHKERRQ(ierr); ierr = VecAssemblyEnd(r); CHKERRQ(ierr);
ierr = SNESSetFunction(snes,r,FormFunction,(void*)F);CHKERRQ(ierr);
ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);
ierr = MatSetFromOptions(J);CHKERRQ(ierr);
if (size==1){ ierr = MatSeqAIJSetPreallocation(J, 3, PETSC_NULL);CHKERRQ(ierr); }
else { ierr = MatMPIAIJSetPreallocation(J, 3, PETSC_NULL, 3, PETSC_NULL);CHKERRQ(ierr); }
ierr = SNESSetJacobian(snes,J,J,FormJacobian,NULL);CHKERRQ(ierr);
ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,400,400,&monP.viewer);CHKERRQ(ierr);
ierr = SNESMonitorSet(snes,Monitor,&monP,0);CHKERRQ(ierr);
ierr = SNESGetKSP(snes, &ksp);CHKERRQ(ierr);
ierr = KSPSetOperators(ksp, J, J);CHKERRQ(ierr);
ierr = KSPGetPC(ksp, &pc);CHKERRQ(ierr);
ierr = KSPSetTolerances(ksp,1.e-4,PETSC_DEFAULT,PETSC_DEFAULT,20);CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
ierr = SNESGetTolerances(snes,&abstol,&rtol,&stol,&maxit,&maxf);CHKERRQ(ierr);
ierr = VecGetOwnershipRange(F,&rstart_F,&rend_F);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"\n (rstart_F, rend_F) = (%i,%i) \n", rstart_F, rend_F);fflush(stdout);
xp = 0.0;
for (i=rstart_F; i<rend_F; i++) {
v = 6.0*xp + PetscPowScalar(xp+1.e-12,6.0); /* +1.e-12 <--- 0 */
ierr = VecSetValues(F,1,&i,&v,INSERT_VALUES);CHKERRQ(ierr);
xp += h;
}
ierr = VecAssemblyBegin(F); CHKERRQ(ierr); ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
ierr = FormInitialGuess(x);CHKERRQ(ierr);
ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr);
ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"number of SNES iterations = %D\n\n",its);CHKERRQ(ierr);
char* printF_write_real = "x_e_real_vec.qxz";
char* printF_write_imag = "x_e_imag_vec.qxz";
ierr = PrintVecToFile_r_i(x, printF_write_real, printF_write_imag,n);CHKERRQ(ierr);
ierr = VecDestroy(&x);CHKERRQ(ierr);
ierr = VecDestroy(&r);CHKERRQ(ierr);
ierr = VecDestroy(&F);CHKERRQ(ierr);
ierr = MatDestroy(&J);CHKERRQ(ierr);
ierr = SNESDestroy(&snes);CHKERRQ(ierr);
ierr = PetscViewerDestroy(&monP.viewer);CHKERRQ(ierr);
ierr = PetscFinalize();
return ierr;
}
/* ------------------------------------------------------------------- */
PetscErrorCode FormInitialGuess(Vec x)
{
PetscErrorCode ierr;
PetscScalar pfive = .90;
PetscInt i, rstart_x, rend_x;
ierr = VecGetOwnershipRange(x,&rstart_x,&rend_x);CHKERRQ(ierr);
for (i=rstart_x; i<rend_x; i++) {
ierr = VecSetValues(x,1,&i,&pfive,INSERT_VALUES);CHKERRQ(ierr);
}
ierr = VecAssemblyBegin(x); CHKERRQ(ierr); ierr = VecAssemblyEnd(x); CHKERRQ(ierr);
return 0;
}
/* ------------------------------------------------------------------- */
PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)
{
Vec g = (Vec)ctx;
const PetscScalar *xx,*gg;
PetscScalar d, f_insval;
PetscErrorCode ierr;
PetscInt i,n,rstart_f,rend_f;
PetscReal normsize_f;
ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
ierr = VecGetArrayRead(g,&gg);CHKERRQ(ierr);
ierr = VecGetSize(x,&n);CHKERRQ(ierr);
d = (PetscReal)(n - 1); d = d*d;
ierr = VecGetOwnershipRange(f,&rstart_f,&rend_f);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"\n ----------------- FormFunction Called -------------------- \n");//fflush(stdout);
PetscPrintf(PETSC_COMM_WORLD,"\n (rstart_f, rend_f) = (%i,%i) \n", rstart_f, rend_f);//fflush(stdout);
for (i=rstart_f; i<rend_f; i++) {
//printf("i = %i \n", i);fflush(stdout);
if (i==0){
f_insval = xx[0] - 0.0;
ierr = VecSetValue(f, i, f_insval, INSERT_VALUES);CHKERRQ(ierr);
}
else if (i==n-1){
f_insval = xx[n-1] - 1.0;
ierr = VecSetValue(f, i, f_insval, INSERT_VALUES);CHKERRQ(ierr);
}
else
{
f_insval = d*(xx[i-1] - 2.0*xx[i] + xx[i+1]) + xx[i]*xx[i] - gg[i];
ierr = VecSetValue(f, i, f_insval, INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = VecAssemblyBegin(f);CHKERRQ(ierr); ierr = VecAssemblyEnd(f);CHKERRQ(ierr);
ierr = VecNorm(f, NORM_2, &normsize_f);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n vecnorm(f) = %.7e\n", normsize_f);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
ierr = VecRestoreArrayRead(g,&gg);CHKERRQ(ierr);
return 0;
}
/* ------------------------------------------------------------------- */
PetscErrorCode FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *dummy)
{
const PetscScalar *xx;
PetscScalar A[3],d;
PetscErrorCode ierr;
PetscInt i,n,j[3];
PetscReal normsize;
ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
ierr = VecGetSize(x,&n);CHKERRQ(ierr);
d = (PetscReal)(n - 1); d = d*d;
PetscInt rstart, rend;
ierr = MatGetOwnershipRange(jac,&rstart,&rend);CHKERRQ(ierr);
PetscPrintf(PETSC_COMM_WORLD,"\n ----------------- FormJacobian Called -------------------- \n");//fflush(stdout);
PetscPrintf(PETSC_COMM_WORLD,"\n (rstart, rend) = (%i,%i) \n", rstart, rend);//fflush(stdout);
for (i=rstart; i<rend; i++) {
//printf("i = %i \n", i);fflush(stdout);
if (i==0){
A[0] = 1.0;
ierr = MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);CHKERRQ(ierr);
}
else if (i==n-1){
A[0] = 1.0;
ierr = MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);CHKERRQ(ierr);
}
else
{
j[0] = i - 1; j[1] = i; j[2] = i + 1;
A[0] = A[2] = d; A[1] = -2.0*d + 2.0*xx[i];
ierr = MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);CHKERRQ(ierr);
}
}
ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
if (jac != B) {
ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
ierr = MatNorm(jac,NORM_FROBENIUS, &normsize);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n matnorm = %.7e\n", normsize);CHKERRQ(ierr);
return 0;
}
PetscErrorCode Monitor(SNES snes,PetscInt its,PetscReal fnorm,void *ctx)
{
PetscErrorCode ierr;
MonitorCtx *monP = (MonitorCtx*) ctx;
Vec x;
ierr = PetscPrintf(PETSC_COMM_WORLD,"iter = %D, SNES Function norm %g\n",its,(double)fnorm);CHKERRQ(ierr);
ierr = SNESGetSolution(snes,&x);CHKERRQ(ierr);
ierr = VecView(x,monP->viewer);CHKERRQ(ierr);
return 0;
}
PetscErrorCode PrintVecToFile_r_i(Vec v, char *filename_r, char *filename_i, int n){
PetscErrorCode ierr;
PetscInt size;
ierr = VecGetSize(v, &size);
PetscScalar *_v;
/*real part*/
FILE *fp;
fp = fopen(filename_r, "w");
ierr = VecGetArray(v,&_v);CHKERRQ(ierr);
for(int i = 0; i < n; i++){
ierr = PetscFPrintf(PETSC_COMM_WORLD,fp, "%.9e\n",creal(_v[i]));CHKERRQ(ierr);
}
ierr = VecRestoreArray(v,&_v);CHKERRQ(ierr);
fclose(fp);
/*imag part*/
FILE *fp_conj;
fp_conj = fopen(filename_i, "w");
ierr = VecGetArray(v,&_v);CHKERRQ(ierr);
for(int i = 0; i < n; i++){
ierr = PetscFPrintf(PETSC_COMM_WORLD,fp, "%.9e\n",cimag(_v[i]));CHKERRQ(ierr);
}
ierr = VecRestoreArray(v,&_v);CHKERRQ(ierr);
fclose(fp_conj);
return 0;
}
install.sh
Description: application/shellscript
