プログラマ専用SNS ミクプラ

CG法での計算例があるのですが、ICCG法の解き方が分かりません。

ご指導のほどよろしくお願いします。

/* CG 法でのポアソン方程式の計算 */
for(i=1;i<=mxz;i++){
v_=w_=0.0;
p_=d_=cnm;
}
for(i=1;i<=mxz;i++){
w_[1]=mA[2][1]*p_[1]+ mA[3][1]*p_[2] mA[4][1]*p_[mx1+1];
for(i=2;i<=mx1;i++)
w_=mA[1]*p_[i-1]+ mA[2]*p_ mA[3]*p_[i+1]+mA[4][i]*p_[i+mx1];
for(i=mx1+1;i<=mxz-mx1;i++)
w_[i]=mA[0][i]*p_[i-mx1]+mA[1][i]*p_[i-1]+mA[2][i]*p_[i]+mA[3][i]*p_[i+1]
38
+mA[4][i]*p_[i+mx1];
for(i=mxz-mx1+1;i<=mxz-1;i++)
w_[i]=mA[0][i]*p_[i-mx1]+mA[1][i]*p_[i-1]+mA[2][i]*p_[i]+mA[3][i]*p_[i+1];
w_[mxz]=mA[0][mxz]*p_[mxz-mx1]+mA[1][mxz]*p_[mxz-1]+mA[2][mxz]*p_[mxz];
eu_=dj_=0.0;
for(j=1;j<=mxz;i++){
dj_+=d_[j]*d_[j];
eu_+=p_[j]*w_[j];
}
if(eu_=0.0)
return;
e_=dj_/eu_;
q_=0.0;
for(j=1;j<=mxz;j++){
v_[j]+=e_*p_[j];
d_[j]-=e_*w_[j];
q_+=d[j]*d[j];
}
q_/=dj_;
for(j=1;j<=mxz;j++)
p_[j]=d_[j]+q_*p_[j];
}
}