並行求解三對角矩陣 CR方法 v2 從0開始，指針傳入

並行求解三對角矩陣 CR方法 v2 從0開始，指針傳入

根據文獻 

A comparison of sequential and parallel elimination methods for tridiagonal matrices

D.J. Evans

int pTM20(double *a, double *b, double *c, double *r,  double *x,int n){
	//原矩陣形式 要求 對角佔優 a[0]=0 c[n]=0
	//b0 c0
	//a1 b1 c1
	//   a2 b2 c2
	//    ···
	//         an bn

	double *f,*g,*h;
	f=new double[n+1];
	h=new double[n+1];
	g=new double[n+1];
	int p=n/2;
#pragma omp parallel sections
{
#pragma omp section
{
	g[0]=c[0]/b[0];
	h[0]=r[0]/b[0];
	for(int i=1;i<=p;i++){
		g[i]=c[i]/(b[i]-a[i]*g[i-1]);
		h[i]=(r[i]-a[i]*h[i-1])/(b[i]-a[i]*g[i-1]);
	}
}
#pragma omp section
{
	f[n]=a[n]/b[n];
	h[n]=r[n]/b[n];
	for(int j=n-1;j>=(p+1);j--){
		f[j]=a[j]/(b[j]-c[j]*f[j+1]);
		h[j]=(r[j]-c[j]*h[j+1])/(b[j]-c[j]*f[j+1]);
	}
}
}
#pragma omp parallel sections
{
#pragma omp section
{
	x[p]=(h[p]-g[p]*h[p+1])/(1-g[p]*f[p+1]);
	for(int i=p-1;i>=0;i--){
		x[i]=h[i]-g[i]*x[i+1];
	}
}
#pragma omp section
{
	x[p+1]=(h[p+1]-f[p+1]*h[p])/(1-g[p]*f[p+1]);
	for(int j=p+2;j<=n;j++){
		x[j]=h[j]-f[j]*x[j-1];
	}
}
}
	delete[] f;
	delete[] g;
	delete[] h;
	return 0;
}