Question from the NR

This is a discussion on Question from the NR within the C Programming forums, part of the General Programming Boards category; I'm trying to solve Ax=b using linbcg.c (p.86 NR in c, 2nd version), A is sparse. All my bs are ...

  1. #1
    Registered User
    Join Date
    Apr 2011
    Posts
    12

    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;
    }
    Last edited by smartfish; 06-12-2011 at 03:45 PM. Reason: indent

  2. #2
    ATH0 quzah's Avatar
    Join Date
    Oct 2001
    Posts
    14,826
    Press edit. Go indent that. No one wants to read a hundred lines of un-indented crap.


    Quzah.
    Hope is the first step on the road to disappointment.

  3. #3
    Registered User
    Join Date
    Apr 2011
    Posts
    12
    sorry, the code was directly copied from the book, didn't notice that.

  4. #4
    Registered User
    Join Date
    Apr 2011
    Posts
    12
    I got the code figured out myself. case closed

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Replies: 1
    Last Post: 03-23-2011, 09:00 AM
  2. *szString = things question/memory question
    By Jang in forum C Programming
    Replies: 3
    Last Post: 01-20-2011, 03:59 AM
  3. Newbish Question file reading question....
    By kas2002 in forum C Programming
    Replies: 23
    Last Post: 05-17-2007, 12:06 PM
  4. Self regiserting DLLs question and libraries question.
    By ApocalypticTime in forum Windows Programming
    Replies: 2
    Last Post: 03-22-2003, 01:02 PM
  5. A question of color...(not a racial question)
    By Sebastiani in forum Windows Programming
    Replies: 7
    Last Post: 01-15-2003, 07:05 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21