並行求解三對角矩陣 Wang方法，分裂算法

本函數求解三對角矩陣 的並行算法 

傳入的數組不是從0開始，是從1開始，數組最大索引爲n

此算法解出的結果有偏差，我實在是找不出有什麼問題，徹底根據Wang的文章編寫

A Parallel Method for Tridiagonal Equations, H. H. WANG

ACM Transactions on Mathematical Software, Vol 7, No. 2, June 1981, Pages 170-183.

int pTM(double a[], double b[], double c[], double r[],  double x[],int n){
	//原矩陣形式 要求 對角佔優 c(1)=0 b(n)=0
	//a1 b1
	//c2 a2 b2
	//   c3 a3 b3
	//    ···
	//         cn an
	//k>2 
	int p,k;
	p=4;
	k=4;
	double *f, *g,*t;
	f=new double[n+1];
	g=new double[n+1];
	t=new double[k+2];
	#pragma omp for
	for(int i=1;i<=p-1;i++){//並行循環
		f[i*k+1]=c[i*k+1];
	}
	//#pragma omp for
	for(int i=0;i<=p-1;i++){//並行循環
		for(int j=2;j<=k;j++)
		{
			c[i*k+j]=c[i*k+j]/a[i*k+j-1];
			a[i*k+j]=a[i*k+j]-c[i*k+j]*b[i*k+j-1];
			r[i*k+j]=r[i*k+j]-c[i*k+j]*r[i*k+j-1];
			if(i!=0){
				f[i*k+j]=-c[i*k+j]*f[i*k+j-1];
			}
		}
	}
	//#pragma omp for
	for(int i=1;i<=p;i++){//並行循環
		g[i*k-1]=b[i*k-1];
	}
	//#pragma omp for
	for(int i=0;i<=p-1;i++){//並行循環
		for(int j=k-1;j>=2;j--)
		{
			b[i*k+j-1]=b[i*k+j-1]/a[i*k+j];
			g[i*k+j-1]=-b[i*k+j-1]*g[i*k+j];
			r[i*k+j-1]=r[i*k+j-1]-b[i*k+j-1]*r[i*k+j];
			if(i!=0){
				f[i*k+j-1]=f[i*k+j-1]-b[i*k+j-1]*f[i*k+j];
			}
		}
	}
	//#pragma omp for
	for(int i=1;i<=p-1;i++){//並行循環
		b[i*k]=b[i*k]/a[i*k+1];
		a[i*k]=a[i*k]-b[i*k]*f[i*k+1];
		g[i*k]=-b[i*k]*g[i*k+1];
		r[i*k]=r[i*k]-b[i*k]*r[i*k+1];
	}
	//#pragma omp for
	for(int i=1;i<=p-1;i++){//並行循環
		for(int j=1;j<=k;j++){
			t[j]=g[i*k+j-1];
		}
		t[k+1]=a[(i+1)*k];
		for(int j=1;j<=k;j++){
			f[i*k+j]=f[i*k+j]/a[i*k];
			t[j+1]=t[j+1]-f[i*k+j]*t[1];
			r[i*k+j]=r[i*k+j]-f[i*k+j]*r[i*k];
		}
		a[(i+1)*k]=t[k+1];
	}
	//#pragma omp for
	for(int i=p-1;i>=1;i--){//並行循環
		for(int j=1;j<=k;j++){
			g[i*k+j-1]=g[i*k+j-1]/a[(i+1)*k];
			r[i*k+j-1]=r[i*k+j-1]-g[i*k+j-1]*r[(i+1)*k];
			if(j!=k){
				g[j]=g[j]/a[k];
				r[j]=r[j]-g[j]*r[k];
			}
		}
	}
	//#pragma omp for
	for(int i=1;i<=n;i++){//並行循環
		x[i]=r[i]/a[i];
	}
	delete[] f;
	delete[] g;
	delete[] t;
	return 0;
}