the idea that communication could occur in assembly could, for example, be an exercise based on vecmatksp.c

e.g. constantine suggests: put rank == 0 conditional around vecsetvalues

The problem is that I have not yet found any online MINPACK source that contains it. Looking for "incidence-degree ordering" (IDO) and "smallest-last ordering" (SLO) as they were included into PETSc source.

Note j[4] is initialized already. It is much clearer code if I do

double v[4][4] = {{1.0, 2.0, 3.0, 0.0},
                  {2.0, 1.0, -2.0, -3.0},
...

note I have already moved the VecDuplicate() which allocates x to the KSP part of the code

Chapter 3: red can be gray; compare a figure in Chapter 5

are their others?

Barry says:

ch1/e.c instead of using PetscPrintf() with PETSC_COMM_SELF and getting scrambled output why not use PetscSynchronizedPrintf()? It only requires one more like the PetscSynchronizedFlush() call. Why have crazy output since you don't need to.

So I need to make my vaguely-conceived point about no-a-priori-message-order more clear.

right now it uses KSP and KSPSetOperators()

Currently it is convert.c; this should be re-done in python. The other C code (new name femfish.c; see #19) reads a petsc binary.

Constantine reports plap.c does not compile in clang. Problem is in these lines in plap.c:

static double zq[3][3] = { {0.0,NAN,NAN},
                           {-1.0/sqrt(3.0),1.0/sqrt(3.0),NAN},
                           {-sqrt(3.0/5.0),0.0,sqrt(3.0/5.0)} },
              wq[3][3] = { {2.0,NAN,NAN},
                           {1.0,1.0,NAN},
                           {5.0/9.0,8.0/9.0,5.0/9.0} };

Two changes here, and in similar: remove sqrt() function call, and (minor) irrelevant static, which I do not want to explain anyway.

triangle will use zero-based if inputs are zero-based

fix

How are you doing it? It looks like you use \inputwhole{} which I have never seen. I am using lstlisting. Is there a reason you did not use this? I would like to standardize, and am willing to listen to argument, but I hate to keep changing.

e.g. in chapter 9 was

test: runadvect_1 runadvect_2 runadvect_3 runad3_1 runad3_2

now

test_advect: runadvect_1 runadvect_2 runadvect_3
test_ad3: runad3_1 runad3_2
test: test_advect test_ad3

From issue #11.

c/ch2/tri.c is an example. Make this uniform so that

$ ./program -help |head          # determine that "prg_" is prefix
$ ./program -help |grep prg_     # list the options with this prefix

is part of the "help system".

Then document this feature in Chapters 2 and 3 where it will first be used.

It needs to work with different gmsh formats. Currently it expects format version 2.2, but that is at least at 4.1 now.

The script already reads a section at a time. Note that from the Gmsh format spec we have

To represent a simple mesh, the minimal sections that should be present in the file are $MeshFormat, $Nodes and $Elements. Nodes are assumed to be defined before elements.

But in Chapter 10 we also need $PhysicalNames. All other sections (e.g. $Entities) are already ignored by the design of msh2petsc.py.

The Gmsh doc site has a page for the legacy format which seems to be 2.2.

Plans:

• Change check_mesh_format() to get_mesh_format() which also gets the version number. (Major,Minor? Or a legacy flag?)
• Pass the version number into `read_SECTION()' functions, and have them act accordingly.

Intro

• PETSc strategy, p6.
  "compose Newton's method and mesh topology". Not clear what this means (to me).
  "choice of preconditioners". ?
• an introduction to the problems that are covered in the book would
  be helpful to get an idea of where the book is going, and to see if I
  should jump to a chapter that interests me (though seasoned people may
  not need this...)

Chapter 1.
p11. "on this machine there no existing MPI installation" (missing "is")

About "facts of life". It wasn't clear to me that they will be coming
up in every chapter, but this may be useful to me beyond PETSc so I
would have liked to see a little more about it.

The unifying theme about facts of life was interesting, but I had to
read it a few times to see what it meant. I think the point is that
they affect stability of simulations.

For now, Facts I -- V are scattered in the book, which is fine. But might as well make them macros which re-appear in a one-page Appendix.

Motivated by #35.

Barry says

I recommend starting an index now. In the PETSc users manual we have a subject index with the latex macro \sindex and a function index with \findex. I think it is better to do it now and make sure you introduce all relevant terms as you along rather than trying to do it all at the end.

Wrong in Chapter 2: "These preconditioners do not communicate any information between rows at all, and thus they exhibit no dependence."

The right thing is along the lines of: "Preconditioners are fast if they avoid communication, i.e. needing info in rows they don't own. Some preconditioners (e.g. jacobi) only use info from single rows, so they are insensitive to the parallel decomposition. Better preconditioners, such as block-Jacobi and ASM, act on multiple rows owned by the processor. Different parallel decompositions therefore determine different actual preconditioners."

Many notes and corrections for the Dec 23 draft
2872_001.pdf

Here are some comments (using revision 58a575b). (You said "I would welcome any kind of feedback", right?)

• tex/preface.tex:8

 and probably some in \Matlab or python

The name of the programming language Python should be capitalized (judging by the way it is spelled on www.python.org).

• tex/preface.tex:12

 Use \texttt{git} to get the example C codes

"Git" should be capitalized also.

• tex/preface.tex:62

 passing arguments by value and reference, 

The C language does not have the concept of a reference: in C we can pass arguments by value or by pointer.

The Wikipedia page on the C language states that Pass-by-reference is simulated in C by explicitly passing pointer values.

• tex/preface.tex:76

 Finally, this book assumes a \emph{bash} shell

According to its official webpage Bash should be capitalized.

• tex/getstarted.tex:16

   We will return to this ``trace-back'' mechanism later in the Chapter.

Add a FIXME or something like \ref{FIXME} (this way LaTeX will keep reminding you about a broken reference).

• tex/getstarted.tex:32

   If it is not already installed on your machine, go to
 \begin{quote}
 \href{http://www.mcs.anl.gov/petsc/download/index.html}{\texttt{www.mcs.anl.gov/petsc/download/}}
 \end{quote}
 to download the source code

Consider creating a \newcommand to make it easier to update URLs such as this one.

• tex/getstarted.tex:124

 A traceback is a list of the nested methods, in reverse order

Technically these are "nested function calls" or "nested stack frames" (there are no methods in C).

• tex/getstarted.tex:131

  on  $N$ processors

Very, very minor: consider using "processor" only when referring to physical processors and "process" otherwise.

It is fine for there to be a separate code to read triangle file and write petsc binary. But the PDE-solving code should be one code.

This is a Barry comment, on Chapter 4, in issue #11.

I have read chapter 5 and I noticed that TS object use the SNES object by default. SNES is a nonlinear solver which calls KSP solver many times in every time step. Is it more expensive than using KSP solver directly?

Uniform code structure and explanations argues for this. Note that the objective evaluation has this outermost loop meaning already. As a result, there will be new if statements inside the loops in FormFunctionLocal().

do it

Barry's comments on this:

ch1/e.c you have chosen to use PetscReal, but not PetscInt, PetscMPIInt etc however in the MPI calls you use MPI_DOUBLE MPI_SUM instead of MPIU_REAL, MPIU_SUM which are actually needed for quad precision; so the benefit of PetscReal doesn't really exist. I suggest actually just using double (and no special PETSc types like PetscReal) and just expand your discussion of complex numbers to say that PETSc can be built to use single, double, quad precisions and 32 or 64 bit integers but for simplicity sake all your codes are just double with 32 bit integers.

Although this gives me "permission" to use double, which would indeed improve the look of my codes, in fact I am not sure I want to give up on having my codes also work in quad precision.

The reason I suggested not doing that is you need the ugly MPIU_stuff (see below) or quad doesn't work and the MPIU_stuff is confusing to new PETSc users. Note that if your code does NO MPI message passing or Allreducing then you don't need any MPIU_stuff.

So: Is there a "how to write a petsc code that is double-->quad portable" doc or man page somewhere?

Sadly, not really. I should add something in the users manual.

For example, I am not sure how switching to quad would affect addressing and int types, if at all, but (I think) you are hinting it does.

int size has nothing to do with quad and doesn't affect addressing either

And direct MPI calls have no flexibility about number types?

The issue is that MPI does not know anything about quad precision at all. Thus MPI_SUM, MPI_MAX etc don't work with quad precision. One needs to use the PETSc provided MPIU_SUM etc and one needs to use MPIU_REAL instead of MPI_DOUBLE; when using quad precision MPIU_REAL is a PETSc defined MPI data type. When using double precision MPIU_REAL is MPI_DOUBLE. One also needs to use a quad precision BLAS/LAPACK; the only one I know about is our own --download-f2cblaslapack.

Barry says

Page 24 you start talking about spectrum and condition number and norms of matrices here without introducing them. I suggest you insert at this point the definition of standard 1,2, infinity vector norms and matrix norms and you can also introduce how they may be computed in PETSc via VecNorm() , VecNormType, etc etc.

Fix.

These are Tables 4.2 and 4.3 in Chapter 4. Constantine suggested this, and it makes sense. Transpose of 4.3 will be especially beneficial.

From issue #11.

I used revision 58a575b.

• tex/linearsystem.tex:152

  Thus if $Ax=b$ for some $x$ then row $i$ of $A$ is on the rank $m$ processor if and if entry $i$ of $b$ is on the rank $m$ processor.

  Should this be "if and only if"?

• tex/linearsystem.tex:558

  Our $A$ is allocated as dense so the exact LU factorization requires no fill-in.

  Define fill-in or provide a reference.

• tex/linearsystem.tex:691

  to time some bigger solves

  Consider saying "some bigger systems". "Solve" is rarely a noun.

• tex/linearsystem.tex:706

  We can try a direct solve

  Should probably say "a direct solver".

• tex/linearsystem.tex:758

  X11 windows are installed.

  Consider saying "if PETSc was built with the X Window System support" (see https://en.wikipedia.org/wiki/X_Window_System).

• tex/linearsystem.tex:794

  For any nontrivial preconditioner, even one based on diagonal blocks of a tridiagonal matrix as here, different sets of rows ``communicate'' at the preconditioning stage depending on the processor count.

  This sentence should probably be re-phrased. (I can't parse it.)

• tex/linearsystem.tex:798

  significant spectral effect.

  Hmm. There must be a better way to say this.

• tex/structured.tex:36

  $c\rho \partial u/\partial t = - \Div\bq + f$

  Add \, (c\rho\, \partial...) to get better spacing.

• tex/structured.tex:243

  i.e.~pass-by-reference

  I.e. via an output argument.

• tex/structured.tex:42

  in the next Chapter we apply a finite element approach

  Reference the chapter. (It's not the next one currently.)

• tex/structured.tex:547

  and we explicitly ask for the \pVec objects to be assembled by calling \texttt{VecAssemblyBegin/End()}

  There is no need to do this if you access Vecs using DMDAVecGetArray() and DMDAVecRestoreArray().

  When VecSetValues() is used, VecAssemblyBegin/End() is needed, though.

Like the "Facts of Life" idea (issue #36), things like "As practical advice for the debugging stage" (Chapter 4) could be a scaffolding to help the reader implement a nontrivial (and not in book) example. Could be a bulleted list in an appendix that repeats (and points to) in-text paragraphs.

The strip.sh fails on my Apple. This patch should be portable I think

$ git diff
diff --git a/tex/strip.sh b/tex/strip.sh
index e094b7b..b982d6a 100755
--- a/tex/strip.sh
+++ b/tex/strip.sh
@@ -43,9 +43,9 @@ cd cstrip/
 set -x

 for NAME in *.c; do
-    sed -i 's/ierr = //' $NAME
-    sed -i 's/CHKERRQ(ierr);//' $NAME
-    sed -i '/\/\/STRIP/d' $NAME
-    sed -i '/PetscErrorCode ierr/d' $NAME
+    sed -i "" 's/ierr = //g' $NAME
+    sed -i "" 's/CHKERRQ(ierr);//g' $NAME
+    sed -i "" '/\/\/STRIP/d' $NAME
+    sed -i "" '/PetscErrorCode ierr/d' $NAME
 done

• The make target rune_1 in c/ch1/ is an example.
• Put outputs in output/ subdirectory.
• No need to be super thorough, but there should be at least one such test for each code.
• No need to test solns/ codes.

Hello Ed,
I am unable to create the .dat files. Any idea what is wrong? I have the latest version of petsc
Thank you for the help.
PS: The book is wonderful !

debasish@debasish-HP-Z840:~/Desktop/softwares/p4pdes/c/ch7$ ./minimal -snes_fd_color -snes_grid_sequence 7 -ksp_view_mat binary::A.dat -ksp_view_rhs binary::b.dat
[0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Unable to open file
[0]PETSC ERROR: Cannot open file for WRITE: H�U��yH�E��
[0]PETSC ERROR: See https://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.15.0, unknown
[0]PETSC ERROR: ./minimal on a linux-c-dbg named debasish-HP-Z840 by debasish Sat Jun 12 12:18:17 2021
[0]PETSC ERROR: Configure options --with-debugging=1
[0]PETSC ERROR: #1 PetscBinaryOpen() at /home/debasish/Desktop/softwares/petsc/src/sys/fileio/sysio.c:505
[0]PETSC ERROR: #2 PetscViewerFileSetUp_BinarySTDIO() at /home/debasish/Desktop/softwares/petsc/src/sys/classes/viewer/impls/binary/binv.c:1460
[0]PETSC ERROR: #3 PetscViewerSetUp_Binary() at /home/debasish/Desktop/softwares/petsc/src/sys/classes/viewer/impls/binary/binv.c:1515
[0]PETSC ERROR: #4 PetscViewerSetUp() at /home/debasish/Desktop/softwares/petsc/src/sys/classes/viewer/interface/view.c:330
[0]PETSC ERROR: #5 PetscOptionsGetViewer() at /home/debasish/Desktop/softwares/petsc/src/sys/classes/viewer/interface/viewreg.c:354
[0]PETSC ERROR: #6 KSPSetFromOptions() at /home/debasish/Desktop/softwares/petsc/src/ksp/ksp/interface/itcl.c:568
[0]PETSC ERROR: #7 SNESSetFromOptions() at /home/debasish/Desktop/softwares/petsc/src/snes/interface/snes.c:1113
[0]PETSC ERROR: #8 main() at /home/debasish/Desktop/softwares/p4pdes/c/ch7/minimal.c:150
[0]PETSC ERROR: PETSc Option Table entries:
[0]PETSC ERROR: -ksp_view_mat binary::A.dat
[0]PETSC ERROR: -ksp_view_rhs binary::b.dat
[0]PETSC ERROR: -snes_fd_color
[0]PETSC ERROR: -snes_grid_sequence 7
[0]PETSC ERROR: ----------------End of Error Message -------send entire error message to [email protected]

MPI_ABORT was invoked on rank 0 in communicator MPI_COMM_WORLD
with errorcode 65.

NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.
You may or may not see output from other processes, depending on
exactly when Open MPI kills them.

This package has been (and truly still is) a very very helpful for learning Petsc for PDEs.
Although i am still exploring the library, I am impatient to find out how to generate an unstructured
simplicial meshes ( triangular and tetrahedral one using Triangle/Tetgen software).
As petsc may install those generators, I would like to know how to give the solver just a *.poly file
and then it will generate the desired mesh.

All the best,
zaq

and a bit of blue/red in text on implementing b.c.s

solid gray looks like black in reproductions

Really looking for this book. Do you have a timeline for the release of this book? :)

do it

Hi,

It seems be to info->my-1 instead of info->my in this line:

content in the new 7. Preconditioners will do the job, but for sure forward refs from linearsystem.tex and structured.tex need fixing

no reason for structuredpoisson.[h|c] to be separated out

Barry says

page 21: "Before PETSc computes the Mat-Vec product, PETSc communicates "scatters" the whole vector to each process". This is wrong in multiple ways it is really "While PETSc computes the Mat-Vec product, PETSc communicators (scatters) the parts of the vector needed by each process using MPI calls"

Fix.

This idea comes from Barry's issue #11 comment on page 20.