Go to the source code of this file.
Functions | |
| void | linbcg (int n, double b[], double x[], int itol, double tol, int itmax) |
|
||||||||||||||||||||||||||||
|
Definition at line 40 of file linbcg.h. References asolve(), atimes(), EPS, err, p, pp, printf(), r, rr, snrm(), z, and zz. Referenced by poisson3D(). 00041 { 00042 int j; 00043 double ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zminrm,znrm; 00044 int iter; 00045 double err=1.e12; 00046 00047 // reset the allocated arrays 00048 for(j=1;j<=n;j++){ 00049 p[j]=0.0; 00050 r[j]=0.0; 00051 z[j]=0.0; 00052 pp[j]=0.0; 00053 rr[j]=0.0; 00054 zz[j]=0.0; 00055 } 00056 00057 // {int i; for(i=1;i<=Ng;i++) printf("b[%d]=%g\n",i,b[i]);} 00058 // calculate the initial residual 00059 iter=0; 00060 atimes(n,x,r,0); 00061 for(j=1;j<=n;j++){ 00062 r[j]=b[j]-r[j]; 00063 rr[j]=r[j]; 00064 // printf("j = %d r = %g b = %g rr = %g\n",j,r[j],b[j],rr[j]); 00065 } 00066 00067 // uncomment this line to get the "minimum residual" variant of the algorithm 00068 // atimes(n,r,rr,0); 00069 00070 if(itol==1){ 00071 bnrm=snrm(n,b,itol); 00072 asolve(n,r,z,0); 00073 } 00074 else if(itol==2){ 00075 asolve(n,b,z,0); 00076 bnrm=snrm(n,z,itol); 00077 asolve(n,r,z,0); 00078 } 00079 else if(itol==3 || itol==4){ 00080 asolve(n,b,z,0); 00081 bnrm=snrm(n,z,itol); 00082 asolve(n,r,z,0); 00083 znrm=snrm(n,z,itol); 00084 } 00085 else{ 00086 printf("linbcg error : illegal itol\n"); 00087 exit(0); 00088 } 00089 00090 // printf("bnrm=%g\n",bnrm); 00091 while(iter<=itmax){ // <--- main loop 00092 // printf("iter = %d itmax = %d\n",iter,itmax); 00093 ++(iter); 00094 asolve(n,rr,zz,1); // final 1 indicates use of transpose matrix A tilde 00095 for(bknum=0.0,j=1;j<=n;j++) bknum+=z[j]*rr[j]; 00096 // printf("bknum=%g\n",bknum); 00097 // calculate coefficient bk and direction vectors p and pp 00098 if(iter==1){ 00099 for(j=1;j<=n;j++){ 00100 p[j]=z[j]; 00101 pp[j]=zz[j]; 00102 } 00103 } 00104 else{ 00105 bk=bknum/bkden; 00106 for(j=1;j<=n;j++){ 00107 p[j]=bk*p[j]+z[j]; 00108 pp[j]=bk*pp[j]+zz[j]; 00109 } 00110 } 00111 00112 // for(j=1;j<=n;j++) printf("p[%d]=%g\n",j,p[j]); 00113 // calculate coefficient ak, new iterate x, and new residuals r and rr 00114 bkden=bknum; 00115 atimes(n,p,z,0); 00116 for(akden=0.0,j=1;j<=n;j++) akden+=z[j]*pp[j]; 00117 // printf("akden=%g\n",akden); 00118 ak=bknum/akden; 00119 // printf("ak=%g\n",ak); 00120 atimes(n,pp,zz,1); 00121 for(j=1;j<=n;j++){ 00122 x[j] += ak*p[j]; 00123 r[j] -= ak*z[j]; 00124 rr[j]-= ak*zz[j]; 00125 } 00126 00127 // solve A tilde . z = r and check stopping criterion 00128 asolve(n,r,z,0); 00129 if(itol==1) err=snrm(n,r,itol)/bnrm; 00130 else if(itol==2) err=snrm(n,z,itol)/bnrm; 00131 else if(itol==3 || itol==4){ 00132 zminrm=znrm; 00133 znrm=snrm(n,z,itol); 00134 if(fabs(zminrm-znrm)>EPS*znrm){ 00135 dxnrm=fabs(ak)*snrm(n,p,itol); 00136 err=znrm/fabs(zminrm-znrm)*dxnrm; 00137 } 00138 else{ 00139 err=znrm/bnrm; // error may not be accurate, so loop again 00140 // printf("iter=%4d err=%12.6f\n",iter,err); 00141 continue; 00142 } 00143 xnrm=snrm(n,x,itol); 00144 if(err<=0.5*xnrm) err/=xnrm; 00145 else{ 00146 err=znrm/bnrm; // error may not be accurate, so loop again 00147 // printf("iter=%4d err=%12.6f\n",iter,err); 00148 continue; 00149 } 00150 } 00151 // printf("*** err=%g ",err); 00152 if(err<=tol){ 00153 printf("\nerr = %g < tol = %g\n",err,tol); 00154 break; 00155 } 00156 } 00157 }
|