www.pudn.com > AdRBF.rar > PBCG.cpp

#include "PBCG.h" 
 
PBCG::PBCG() 
{ 
	sa = NULL; 
	ija = NULL; 
} 
 
PBCG::~PBCG() 
{ 
	if(sa != NULL) 
		delete sa; 
	if(ija != NULL) 
		delete ija; 
} 
 
void PBCG::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=dvector(1,n); 
	pp=dvector(1,n); 
	r=dvector(1,n); 
	rr=dvector(1,n); 
	z=dvector(1,n); 
	zz=dvector(1,n); 
 
	*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 nrerror("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("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b\biter=%4d err=%12.6f",*iter,*err/diagonal); 
          if (*err <= tol){ 
            printf("\n"); 
            break; 
          } 
	} 
 
	free_dvector(p,1,n); 
	free_dvector(pp,1,n); 
	free_dvector(r,1,n); 
	free_dvector(rr,1,n); 
	free_dvector(z,1,n); 
	free_dvector(zz,1,n); 
} 
 
void PBCG::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 PBCG::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 PBCG::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 PBCG::dsprsax(double sa[], unsigned long ija[], double x[], double b[], unsigned long n) 
{ 
	unsigned long i,k; 
 
	if (ija[1] != n+2) nrerror("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 PBCG::dsprstx(double sa[], unsigned long ija[], double x[], double b[], unsigned long n) 
{ 
	unsigned long i,j,k; 
	if (ija[1] != n+2) nrerror("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]; 
		} 
	} 
}