#include<stdio.h>#include<math.h>#include"PCBCGSolver.h"#define EPS 1.0e-14void PCBCGSolver::dsprsax(double sa[], unsigned long ija[], double x[], double b[], unsigned long n){		unsigned long i,k;	if (ija[1] != n+2) printf("dsprsax: mismatched vector and matrix");	for (i=1;i<=n;i++) {		b[i]=sa[i]*x[i];		for (k=ija[i];k<=ija[i+1]-1;k++) b[i] += sa[k]*x[ija[k]];	}}void PCBCGSolver::dsprstx(double sa[], unsigned long ija[], double x[], double b[], unsigned long n){		unsigned long i,j,k;	if (ija[1] != n+2) printf("mismatched vector and matrix in dsprstx");	for (i=1;i<=n;i++) b[i]=sa[i]*x[i];	for (i=1;i<=n;i++) {		for (k=ija[i];k<=ija[i+1]-1;k++) {			j=ija[k];			b[j] += sa[k]*x[i];		}	}}void PCBCGSolver::asolve(unsigned long n, double b[], double x[], int itrnsp){	unsigned long i;	for(i=1;i<=n;i++) x[i]=(sa[i] != 0.0 ? b[i]/sa[i] : b[i]);}void PCBCGSolver::atimes(unsigned long n, double x[], double r[], int itrnsp){  if (itrnsp) dsprstx(sa,ija,x,r,n);  else dsprsax(sa,ija,x,r,n);}double PCBCGSolver::snrm(unsigned long n, double sx[], int itol){	unsigned long i,isamax;	double ans;	if (itol <= 3) {		ans = 0.0;		for (i=1;i<=n;i++) ans += sx[i]*sx[i];		return sqrt(ans);	} else {		isamax=1;		for (i=1;i<=n;i++) {			if (fabs(sx[i]) > fabs(sx[isamax])) isamax=i;		}		return fabs(sx[isamax]);	}}void PCBCGSolver::linbcg(unsigned long n, double b[], double x[], int itol, double tol,	int itmax, int *iter, double *err){					unsigned long j;	double ak,akden,bk,bkden,bknum,bnrm,dxnrm,xnrm,zm1nrm,znrm;	double *p,*pp,*r,*rr,*z,*zz;	p=new double[n+1];	pp=new double[n+1];	r=new double[n+1];	rr=new double[n+1];	z=new double[n+1];	zz=new double[n+1];	*iter=0;	atimes(n,x,r,0);	for (j=1;j<=n;j++) {		r[j]=b[j]-r[j];		rr[j]=r[j];	}	znrm=1.0;	if (itol == 1) bnrm=snrm(n,b,itol);	else if (itol == 2) {		asolve(n,b,z,0);		bnrm=snrm(n,z,itol);	}	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 printf("illegal itol in linbcg");	asolve(n,r,z,0);	while (*iter <= itmax) {		++(*iter);		zm1nrm=znrm;		asolve(n,rr,zz,1);		for (bknum=0.0,j=1;j<=n;j++) bknum += z[j]*rr[j];		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;		atimes(n,p,z,0);		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);		if (itol == 1 || itol == 2) {			znrm=1.0;			*err=snrm(n,r,itol)/bnrm;		} else if (itol == 3 || itol == 4) {			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;				continue;			}			xnrm=snrm(n,x,itol);			if (*err <= 0.5*xnrm) *err /= xnrm;			else {				*err=znrm/bnrm;				continue;			}		}		//printf("iter=%4d err=%12.6f\n",*iter,*err);	if (*err <= tol) break;	}	delete [] p;	delete [] pp;	delete [] r;	delete [] rr;	delete [] z;	delete [] zz;}#undef EPS