On Thu, Jul 9, 2015 at 11:02 PM, Marco Zocca <[email protected]> wrote:
> 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. > You can use MatSetValuesStencil() as is done here: http://www.mcs.anl.gov/petsc/petsc-current/src/ksp/ksp/examples/tutorials/ex50.c.html on line 162. > 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: > Actually, this is only common is toy example. Almost no real PDEs have a fixed stencil because you have material coefficients, or nonlinearities, or constraints, etc. Matt > 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 > -- 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
