Strongly Implicit Solver (SIS) 代碼

這兩個函數分別解5對角 和 9對角矩陣

根據文獻

a Strongly Implicit Solver for Two-Dimensional Elliptic Differential Equations

Lee, Shong-Leih

Numerical Heat Transfer, Part B: Fundamentals, vol. 16, issue 2, pp. 161-178

int SIS5solver(double *Ab, double *Ad, double *Ae, 
			  double *Af, double *Ah, 
			  double *x,double* r, int m,int n)
{
	double *La,*Lb,*Lc,*Ld,*Ue,*Uf,*Ug;
	double *ALF,*BAT,*H,*RES;
	double RESN;
	double SOR=1;
	double TOL=1E-6;
	bool COV=false;
	int MAXITER=1000;
	int ITER;
	La=new double [n+1];
	Lb=new double [n+1];
	Lc=new double [n+1];
	Ld=new double [n+1];
	Ue=new double [n+1];
	Uf=new double [n+1];
	Ug=new double [n+1];
	ALF=new double [n+1];
	BAT=new double [n+1];
	H=new double [n+1];
	RES=new double [n+1];
	for(int i=0;i<=n;i++){
		La[i]=0;
		Lb[i]=0;
		Lc[i]=0;
		Ld[i]=0;
		Uf[i]=0;
		Ug[i]=0;
		ALF[i]=0;
		BAT[i]=0;
		H[i]=0;
		RES[i]=0;
	}
	for(int i=0;i<=n;i++)
	{
		//La-----------------
		
		//Lb-----------------
		if(i>=m){
			Lb[i]=Ab[i] - La[i]*Ue[i-m-1];
		}else{
			Lb[i]=0;
		}
		//Lc-----------------
		if(i>=m-1){
			ALF[i]=Lb[i]*Ue[i-m];
		}else{
			ALF[i]= 0;
		}
		//Ld-----------------
		if(i>=1){
			Lc[i]=Ad[i] - La[i]*Uf[i-m-1];
		}else{
			Lc[i]=0;
		}
		//Le-----------------
		Ld[i]=Ae[i] - La[i]*Ug[i-m-1] - Lb[i]*Uf[i-m] - Lc[i]*Ue[i-1];
		//Uf-----------------
		if(i<=n-1){	
			Ue[i]=(Af[i] - Lb[i]*Ug[i-m])/Ld[i];
		}else{
			Ue[i]=0;
		}
		//Ug-----------------
		if(i<=n-m+1){
			BAT[i]=Lc[i]*Uf[i-1];
		}else{
			BAT[i]=0;
		}
		//Uh-----------------
		if(i<=n-m){
			Uf[i]=(Ah[i] - Lc[i]*Ug[i-1])/Ld[i];
		}else{
			Uf[i]=0;
		}
		//Us-----------------
		
	}
	ITER=0;
	do{
		for(int i=0;i<=n;i++)
		{
			H[i]=r[i]+ALF[i]*(i-m+1>=0?x[i-m+1]:0)+BAT[i]*(i+m-1<=n?x[i+m-1]:0);
		}

		for(int i=0;i<=n;i++)
		{
			H[i]=(H[i]-La[i]*(i-m-1>=0?H[i-m-1]:0)
					-Lb[i]*(i-m>=0?H[i-m]:0)
					-Lc[i]*(i-1>=0?H[i-1]:0))/Ld[i];
		}
		
		for(int i=n;i>=0;i--)
		{
			H[i]=(H[i]-Ue[i]*(i+1<=n?H[i+1]:0)
					-Uf[i]*(i+m<=n?H[i+m]:0)
					-Ug[i]*(i+m+1<=n?H[i+m+1]:0));
		}
		for(int i=0;i<=n;i++)
		{
			x[i]=x[i]+SOR*(H[i]-x[i]);
		}
		RESN=0;
		for(int i=0;i<=n;i++)
		{
			RES[i]=r[i]-(+Ab[i]*(i-m>=0?x[i-m]:0)
						+Ad[i]*(i-1>=0?x[i-1]:0)
						+Ae[i]*x[i]
						+Af[i]*(i+1<=n?x[i+1]:0)
						+Ah[i]*(i+m<=n?x[i+m]:0));
			RESN+=abs(RES[i]);
		}
		ITER++;
	}while(abs(RESN)>1E-10&&ITER<MAXITER);
	
	delete [] La;
	delete [] Lb;
	delete [] Lc;
	delete [] Ld;
	delete [] Ue;
	delete [] Uf;
	delete [] Ug;
	delete [] ALF;
	delete [] BAT;
	return ITER>MAXITER?-1:ITER;
}

int SIS9solver(double *Aa, double *Ab, double *Ac, double *Ad, double *Ae, 
			  double *Af, double *Ag, double *Ah, double *As, 
			  double *x,double* r, int m,int n)
{
	double *La,*Lb,*Lc,*Ld,*Ue,*Uf,*Ug;
	double *ALF,*BAT,*H,*RES;
	double RESN;
	double SOR=1;
	double TOL=1E-6;
	bool COV=false;
	int MAXITER=1000;
	int ITER;
	La=new double [n+1];
	Lb=new double [n+1];
	Lc=new double [n+1];
	Ld=new double [n+1];
	Ue=new double [n+1];
	Uf=new double [n+1];
	Ug=new double [n+1];
	ALF=new double [n+1];
	BAT=new double [n+1];
	H=new double [n+1];
	RES=new double [n+1];
	for(int i=0;i<=n;i++){
		La[i]=0;
		Lb[i]=0;
		Lc[i]=0;
		Ld[i]=0;
		Uf[i]=0;
		Ug[i]=0;
		ALF[i]=0;
		BAT[i]=0;
		H[i]=0;
		RES[i]=0;
	}
	for(int i=0;i<=n;i++)
	{
		//La-----------------
		if(i>=m+1){
			La[i]= Aa[i];
		}else{
			La[i]=0;
		}
		//Lb-----------------
		if(i>=m){
			Lb[i]=Ab[i] - La[i]*Ue[i-m-1];
		}else{
			Lb[i]=0;
		}
		//Lc-----------------
		if(i>=m-1){
			ALF[i]=Lb[i]*Ue[i-m] - Ac[i];
		}else{
			ALF[i]= 0;
		}
		//Ld-----------------
		if(i>=1){
			Lc[i]=Ad[i] - La[i]*Uf[i-m-1];
		}else{
			Lc[i]=0;
		}
		//Le-----------------
		Ld[i]=Ae[i] - La[i]*Ug[i-m-1] - Lb[i]*Uf[i-m] - Lc[i]*Ue[i-1];
		//Uf-----------------
		if(i<=n-1){	
			Ue[i]=(Af[i] - Lb[i]*Ug[i-m])/Ld[i];
		}else{
			Ue[i]=0;
		}
		//Ug-----------------
		if(i<=n-m+1){
			BAT[i]=Lc[i]*Uf[i-1] - Ag[i];
		}else{
			BAT[i]=0;
		}
		//Uh-----------------
		if(i<=n-m){
			Uf[i]=(Ah[i] - Lc[i]*Ug[i-1])/Ld[i];
		}else{
			Uf[i]=0;
		}
		//Us-----------------
		if(i<=n-m-1){
			Ug[i]=As[i]/Ld[i];
		}else{
			Ug[i]= 0;
		}
	}
	ITER=0;
	do{
		for(int i=0;i<=n;i++)
		{
			H[i]=r[i]+ALF[i]*(i-m+1>=0?x[i-m+1]:0)+BAT[i]*(i+m-1<=n?x[i+m-1]:0);
		}

		for(int i=0;i<=n;i++)
		{
			H[i]=(H[i]-La[i]*(i-m-1>=0?H[i-m-1]:0)
					-Lb[i]*(i-m>=0?H[i-m]:0)
					-Lc[i]*(i-1>=0?H[i-1]:0))/Ld[i];
		}
		
		for(int i=n;i>=0;i--)
		{
			H[i]=(H[i]-Ue[i]*(i+1<=n?H[i+1]:0)
					-Uf[i]*(i+m<=n?H[i+m]:0)
					-Ug[i]*(i+m+1<=n?H[i+m+1]:0));
		}
		for(int i=0;i<=n;i++)
		{
			x[i]=x[i]+SOR*(H[i]-x[i]);
		}
		RESN=0;
		for(int i=0;i<=n;i++)
		{
			RES[i]=r[i]-(Aa[i]*(i-m-1>=0?x[i-m-1]:0)
						+Ab[i]*(i-m>=0?x[i-m]:0)
						+Ac[i]*(i-m+1>=0?x[i-m+1]:0)
						+Ad[i]*(i-1>=0?x[i-1]:0)
						+Ae[i]*x[i]
						+Af[i]*(i+1<=n?x[i+1]:0)
						+Ag[i]*(i+m-1<=n?x[i+m-1]:0)
						+Ah[i]*(i+m<=n?x[i+m]:0)
						+As[i]*(i+m+1<=n?x[i+m+1]:0));
			RESN+=abs(RES[i]);
		}
		ITER++;
	}while(abs(RESN)>1E-10&&ITER<MAXITER);
	
	delete [] La;
	delete [] Lb;
	delete [] Lc;
	delete [] Ld;
	delete [] Ue;
	delete [] Uf;
	delete [] Ug;
	delete [] ALF;
	delete [] BAT;
	return ITER>MAXITER?-1:ITER;
}