Bookkeeping

Important note: The Numerical Recipes routines for linear systems are all ``one based'' (all array and matrix elements begin with one, not zero). The simplest way to work with this is to build all of our linear algebra objects (such as G) also as one based and to keep our physics grid objects (such as phi) as zero based. This mixed strategy allows us to use the Numerical Recipes routines with minimal modification and requires no changes to our preexisting code.

The new algorithm requires that we keep track of the charge densities Rho and their residuals g from all previous iterations. Accordingly, we shall require two dimensional arrays (dmatrix's) for their storage in the main program.

To do this, modify the main program to store the Rho and g as it goes along, producing the new program ``test_it''. An appropriate code fragment is

#define Itmx 100
            .
            .
  Rho=dmatrix(1,Itmx+1,0,Nmx);
  g=dmatrix(1,Itmx+1,0,Nmx);
            .
            .
  for (it=1; it<=Itmx; it++) {
    Etot=getg(g[it],Rho[it],Z,lmax,nmax,nmaxmax,F,r,dr,N);
    printf("Iter:%4d  Etot:%20.12f  Error:%20.12f\n",
           it,Etot,fabs(Etot-(-74.473076803203738)));

    for (k=0; k<=N; k++)
      Rho[it+1][k]=Rho[it][k]+alpha*g[it][k];
  }
            .
            .
  free_dmatrix(Rho,1,Itmx+1,0,Nmx);
  free_dmatrix(g,1,Itmx+1,0,Nmx);

Note how we begin the matrices Rho and g with it=1 so that we can conform to the one-based routines from Numerical Recipes. Also, because the final iteration (Itmx) needs space for Rho[it+1], we declare space for 1..Itmx+1 in the dmatrix statements.

• Debugging