Thread: Question from the NR

I'm trying to solve Ax=b using linbcg.c (p.86 NR in c, 2nd version), A is sparse. All my bs are 0s. The final solution x are suppose to be all 0s, too. But this leads to the problem that when initial guessing is 0s. All z's and rr's are zero and bkden are zeros.This made the code crash. How do I overcome that? (below are directly copied from NR.

Code:

*iter=0;
atimes(n,x,r,0); //Input to atimes is x[1..n], output is r[1..n];
for (j=1;j<=n;j++) 
{
     r[j]=b[j]-r[j];
    rr[j]=r[j];
}
/* atimes(n,r,rr,0); */ 
//Uncomment this line to get the “minimum residual”
if (itol == 1)
{ //variant of the algorithm.
   bnrm=snrm(n,b,itol);
   asolve(n,r,z,0); //Input to asolve is r[1..n], output is z[1..n]; the final 0 indicates that the matrix A (not its transpose) is to be used.
}
else if (itol == 2)
 {
          asolve(n,b,z,0);
          bnrm=snrm(n,z,itol);
          asolve(n,r,z,0);
  }
else if (itol == 3 || itol == 4) 
{
          asolve(n,b,z,0);
          bnrm=snrm(n,z,itol);
          asolve(n,r,z,0);
          znrm=snrm(n,z,itol);
}
 else nrerror("illegal itol in linbcg");
while (*iter <= itmax) 
{
          // Main loop.
          ++(*iter);
          asolve(n,rr,zz,1);//  Final 1 indicates use of transpose matrix AT.
for (bknum=0.0,j=1;j<=n;j++) bknum += z[j]*rr[j];
//Calculate coefficient bk and direction vectors p and pp.
if (*iter == 1)
{
   for (j=1;j<=n;j++)
   {
        p[j]=z[j];
        pp[j]=zz[j];
   }
}
else
 {
       bk=bknum/bkden;
       for (j=1;j<=n;j++) 
            {
             p[j]=bk*p[j]+z[j];
            pp[j]=bk*pp[j]+zz[j];
             }
 }
bkden=bknum; //Calculate coefficient ak, newi terate x, and new
atimes(n,p,z,0); //residuals r and rr.
for (akden=0.0,j=1;j<=n;j++) akden += z[j]*pp[j];
ak=bknum/akden;
atimes(n,pp,zz,1);
for (j=1;j<=n;j++) 
{
     x[j] += ak*p[j];
     r[j] -= ak*z[j];
    rr[j] -= ak*zz[j];
}
asolve(n,r,z,0); //Solve A · z = r and check stopping criterion.
if (itol == 1)
    *err=snrm(n,r,itol)/bnrm;
else if (itol == 2)
    *err=snrm(n,z,itol)/bnrm;
else if (itol == 3 || itol == 4) 
{
          zm1nrm=znrm;
          znrm=snrm(n,z,itol);
          if (fabs(zm1nrm-znrm) > EPS*znrm)
         {
              dxnrm=fabs(ak)*snrm(n,p,itol);
             *err=znrm/fabs(zm1nrm-znrm)*dxnrm;
          }
        else 
         {
             *err=znrm/bnrm; Error may not be accurate, so loop again.
             continue;
         }
         xnrm=snrm(n,x,itol);
        if (*err <= 0.5*xnrm) *err /= xnrm;
        else 
        {
            *err=znrm/bnrm; Error may not be accurate, so loop again.
             continue;
         }
}
printf("iter=%4d err=%12.6f\n",*iter,*err);
if (*err <= tol) break;
}

Press edit. Go indent that. No one wants to read a hundred lines of un-indented crap.


Quzah.