[petsc-users] human-readable examples and PETSc idioms
Marco Zocca
zocca.marco at gmail.com
Thu Jul 9 23:02:35 CDT 2015
Hello,
I'm looking for some non-optimized PETSc code; namely, I struggle a
bit with generalizing the provided examples. On one hand inlined
elementary operations help keep track of the FLOPs but make for
hard-to read code.
For instance, KSP example 2, (
http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex2.c.html
), isn't there a better abstraction to defining the finite-difference
grid Laplacian than the following?
for (Ii=Istart; Ii<Iend; Ii++) {
v = -1.0; i = Ii/n; j = Ii - i*n;
if (i>0) {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
if (j>0) {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
}
I see this approach as very error prone, yet this is one of the
simplest examples.
Another problematic trend I see is hardcoding constants and
assumptions on the functions, e.g. in the following (taken from
http://www.mcs.anl.gov/petsc/petsc-current/src/snes/examples/tutorials/ex2.c.html
, removed comments for space) we are interested in forming a matrix A
from the values contained in a vector x; i.e. some 'map' operation
that sweeps x and prepares rows of A at a time according to some fixed
stencil, with separate treatment of boundary values .
Isn't there an idiom to abstract out this functionality? It's one of
the most common operations in numerical PDE codes:
VecGetArrayRead(x,&xx);
VecGetSize(x,&n);
d = (PetscReal)(n - 1); d = d*d;
for (i=1; i<n-1; i++) {
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];
MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES);
}
i = 0; A[0] = 1.0;
MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
i = n-1; A[0] = 1.0;
MatSetValues(jac,1,&i,1,&i,A,INSERT_VALUES);
VecRestoreArrayRead(x,&xx);
I think it would be extremely helpful if the developers of the
examples spent a little time in adding theory annotations to the code,
and I would be very thankful if anybody could point me to relevant
tutorials and literature.
Thank you in advance and kind regards,
Marco Zocca
More information about the petsc-users
mailing list