Tridiagonal matrix algorithm - Revision history

Cedar101: /* Variants */ $, {{mvar}}$

2025-01-14T05:10:08Z

Variants: <math display="block">, {{mvar}}

← Previous revision		Revision as of 05:10, 14 January 2025
Line 176:		Line 176:
	In some situations, particularly those involving [[periodic boundary conditions]], a slightly perturbed form of the tridiagonal system may need to be solved:		In some situations, particularly those involving [[periodic boundary conditions]], a slightly perturbed form of the tridiagonal system may need to be solved:

	:<math>		<math display="block">
	\begin{alignat}{4}		\begin{alignat}{4}
	& a_1 x_{n} && + b_1 x_1 && + c_1 x_2 && = d_1 \\		& a_1 x_{n} && + b_1 x_1 && + c_1 x_2 && = d_1 \\
Line 242:		Line 242:
	.		.
	</math>		</math>
	Where <math> \gamma\in\mathbb{R} </math> is a parameter to be chosen. The matrix ~~<math>~~ A ~~</math>~~ can be reconstructed as <math> A = B + u v^\mathsf{T}</math>. The solution is then obtained in the following way:<ref>{{Cite journal \|last1=Batista \|first1=Milan \|last2=Ibrahim Karawia \|first2=Abdel Rahman A. \|date=2009 \|title=The use of the Sherman–Morrison–Woodbury formula to solve cyclic block tri-diagonal and cyclic block penta-diagonal linear systems of equations \|url=http://dx.doi.org/10.1016/j.amc.2009.01.003 \|journal=Applied Mathematics and Computation \|volume=210 \|issue=2 \|pages=558–563 \|doi=10.1016/j.amc.2009.01.003 \|issn=0096-3003}}</ref> first we solve two tridiagonal systems of equations applying the Thomas algorithm:		Where <math> \gamma\in\mathbb{R} </math> is a parameter to be chosen. The matrix {{mvar\|A}} can be reconstructed as <math> A = B + u v^\mathsf{T}</math>. The solution is then obtained in the following way:<ref>{{Cite journal \|last1=Batista \|first1=Milan \|last2=Ibrahim Karawia \|first2=Abdel Rahman A. \|date=2009 \|title=The use of the Sherman–Morrison–Woodbury formula to solve cyclic block tri-diagonal and cyclic block penta-diagonal linear systems of equations \|url=http://dx.doi.org/10.1016/j.amc.2009.01.003 \|journal=Applied Mathematics and Computation \|volume=210 \|issue=2 \|pages=558–563 \|doi=10.1016/j.amc.2009.01.003 \|issn=0096-3003}}</ref> first we solve two tridiagonal systems of equations applying the Thomas algorithm:
	<math display="block">By=d \qquad \qquad Bq=u </math>		<math display="block">By=d \qquad \qquad Bq=u </math>

	Then we reconstruct the solution ~~<math>~~x ~~</math>~~ using the [[Sherman–Morrison formula\|Shermann-Morrison formula]]:		Then we reconstruct the solution {{mvar\|x}} using the [[Sherman–Morrison formula\|Shermann-Morrison formula]]:
	<math display="block">\begin{align}		<math display="block">\begin{align}
	x &=A^{-1}d		x &=A^{-1}d
Line 255:		Line 255:

	The implementation as a [[C (programming language)\|C]] function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused:		The implementation as a [[C (programming language)\|C]] function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused:
	<syntaxhighlight lang="c">		<syntaxhighlight lang="c" line>
	void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double u[restrict X]) {		void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double u[restrict X]) {
	/*		/*
Line 332:		Line 332:
	<math display="block">x_i = u_i + x_1 v_i \qquad i=1,2,\dots, n</math>		<math display="block">x_i = u_i + x_1 v_i \qquad i=1,2,\dots, n</math>

	Indeed, multiplying each equation of the second auxiliary system by <math>x_1</math>, adding with the corresponding equation of the first auxiliary system and using the representation <math>x_i = u_i + x_1 v_i</math>, we immediately see that equations number ~~<math>~~2~~</math>~~ through ~~<math>~~n~~</math>~~ of the original system are satisfied; it only remains to satisfy equation number ~~<math>~~1~~</math>~~. To do so, consider formula for <math>i=2</math> and <math>i=n</math> and substitute <math>x_2 = u_2 + x_1 v_2</math>and <math>x_n = u_n + x_1 v_n</math> into the first equation of the original system. This yields one scalar equation for <math>x_1</math>:		Indeed, multiplying each equation of the second auxiliary system by <math>x_1</math>, adding with the corresponding equation of the first auxiliary system and using the representation <math>x_i = u_i + x_1 v_i</math>, we immediately see that equations number {{val\|2}} through {{mvar\|n}} of the original system are satisfied; it only remains to satisfy equation number {{val\|1}}. To do so, consider formula for <math>i=2</math> and <math>i=n</math> and substitute <math>x_2 = u_2 + x_1 v_2</math>and <math>x_n = u_n + x_1 v_n</math> into the first equation of the original system. This yields one scalar equation for <math>x_1</math>:
	<math display="block"> b_1x_1+c_1(u_2+x_1v_2)+a_1(u_n+x_1v_n) = d_1</math>		<math display="block"> b_1x_1+c_1(u_2+x_1v_2)+a_1(u_n+x_1v_n) = d_1</math>

Line 339:		Line 339:

	The implementation as a [[C (programming language)\|C]] function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused:		The implementation as a [[C (programming language)\|C]] function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused:
	<syntaxhighlight lang="c">		<syntaxhighlight lang="c" line>
	void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double v[restrict X]) {		void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double v[restrict X]) {
	/* first solve a system of length X - 1 for two right hand sides, ignoring ix == 0 */		/* first solve a system of length X - 1 for two right hand sides, ignoring ix == 0 */
Line 374:		Line 374:
	</syntaxhighlight>		</syntaxhighlight>

	In both cases the auxiliary systems to be solved are genuinely tri-diagonal, so the overall computational complexity of solving system <math> Ax = d </math> remains linear with the respect to the dimension of the system ~~<math>~~ n ~~</math>~~, that is <math>O(n)</math> arithmetic operations.		In both cases the auxiliary systems to be solved are genuinely tri-diagonal, so the overall computational complexity of solving system <math> Ax = d </math> remains linear with the respect to the dimension of the system {{mvar\|n}}, that is <math>O(n)</math> arithmetic operations.

	In other situations, the system of equations may be '''block tridiagonal''' (see [[block matrix]]), with smaller submatrices arranged as the individual elements in the above matrix system (e.g., the 2D [[Poisson equation discretized into block tridiagonal\|Poisson problem]]). Simplified forms of Gaussian elimination have been developed for these situations.<ref name="Quarteroni2007">{{cite book\|author-link1= Alfio Quarteroni\| last1=Quarteroni\|first1=Alfio \|last2=Sacco\|first2=Riccardo \| last3=Saleri\|first3=Fausto \| year=2007\|title=Numerical Mathematics\|publisher= Springer, New York\|isbn= 978-3-540-34658-6 \| chapter=Section 3.8}}</ref>		In other situations, the system of equations may be '''block tridiagonal''' (see [[block matrix]]), with smaller submatrices arranged as the individual elements in the above matrix system (e.g., the 2D [[Poisson equation discretized into block tridiagonal\|Poisson problem]]). Simplified forms of Gaussian elimination have been developed for these situations.<ref name="Quarteroni2007">{{cite book\|author-link1= Alfio Quarteroni\| last1=Quarteroni\|first1=Alfio \|last2=Sacco\|first2=Riccardo \| last3=Saleri\|first3=Fausto \| year=2007\|title=Numerical Mathematics\|publisher= Springer, New York\|isbn= 978-3-540-34658-6 \| chapter=Section 3.8}}</ref>

46.1.62.187: /* Method */ Fixed typo in code comment: v should be d

2025-01-13T14:46:48Z

Method: Fixed typo in code comment: v should be d

← Previous revision		Revision as of 14:46, 13 January 2025
Line 78:		Line 78:
	solves Ax = d, where A is a tridiagonal matrix consisting of vectors a, b, c		solves Ax = d, where A is a tridiagonal matrix consisting of vectors a, b, c
	X = number of equations		X = number of equations
	x[] = initially contains the input v, and returns x. indexed from [0, ..., X - 1]		x[] = initially contains the input, d, and returns x. indexed from [0, ..., X - 1]
	a[] = subdiagonal, indexed from [1, ..., X - 1]		a[] = subdiagonal, indexed from [1, ..., X - 1]
	b[] = main diagonal, indexed from [0, ..., X - 1]		b[] = main diagonal, indexed from [0, ..., X - 1]

JAELloyd: In first implementation of cyclic_thomas, the computation of 'fact' used beta/gamma instead of the alpha/gamma.

2024-11-29T22:35:15Z

In first implementation of cyclic_thomas, the computation of 'fact' used beta/gamma instead of the alpha/gamma.

← Previous revision		Revision as of 22:35, 29 November 2024
Line 297:		Line 297:
	}		}

	const double fact = (x[0] + x[X - 1] * ~~beta~~ / gamma) / (1.0 + u[0] + u[X - 1] * ~~beta~~ / gamma);		const double fact = (x[0] + x[X - 1] * alpha / gamma) / (1.0 + u[0] + u[X - 1] * alpha / gamma);

	/* loop from 0 to X - 1 inclusive */		/* loop from 0 to X - 1 inclusive */

JAELloyd: Corrected the code for cyclic_thomas. cmod[0] and m (under "handle X - 1") were not being initialized correctly. This is based on my having just implemented and tested this code.

2024-11-25T06:53:43Z

Corrected the code for cyclic_thomas. cmod[0] and m (under "handle X - 1") were not being initialized correctly. This is based on my having just implemented and tested this code.

← Previous revision		Revision as of 06:53, 25 November 2024
Line 274:		Line 274:
	const double gamma = -b[0];		const double gamma = -b[0];

	cmod[0] = ~~alpha~~ / (b[0] - gamma);		cmod[0] = c[0] / (b[0] - gamma);
	u[0] = gamma / (b[0] - gamma);		u[0] = gamma / (b[0] - gamma);
	x[0] /= (b[0] - gamma);		x[0] /= (b[0] - gamma);
Line 287:		Line 287:

	/* handle X - 1 */		/* handle X - 1 */
	const double m = 1.0 / (b[X - 1] - alpha * beta / gamma - ~~beta~~ * cmod[X - 2]);		const double m = 1.0 / (b[X - 1] - alpha * beta / gamma - a[X - 1] * cmod[X - 2]);
	u[X - 1] = (alpha - a[X - 1] * u[X - 2]) * m;		u[X - 1] = (alpha - a[X - 1] * u[X - 2]) * m;
	x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;		x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;

169.231.25.141: The previous code wouldn't have run. You can see this because the loop is from ix = 1 to X-1 inclusive, and accesses c[ix], but c only has elements from 0 to X-2. This causes an error on the ix = X-1 iteration, however, note that scratch[X-1] is never used in the computation, so we can simply only execute the line when ix is less than X-1.

2024-04-08T21:00:54Z

The previous code wouldn't have run. You can see this because the loop is from ix = 1 to X-1 inclusive, and accesses c[ix], but c only has elements from 0 to X-2. This causes an error on the ix = X-1 iteration, however, note that scratch[X-1] is never used in the computation, so we can simply only execute the line when ix is less than X-1.

← Previous revision		Revision as of 21:00, 8 April 2024
Line 90:		Line 90:
	/* loop from 1 to X - 1 inclusive */		/* loop from 1 to X - 1 inclusive */
	for (int ix = 1; ix < X; ix++) {		for (int ix = 1; ix < X; ix++) {
			if (ix < X-1){
	scratch[ix] = c[ix] / (b[ix] - a[ix] * scratch[ix - 1]);		scratch[ix] = c[ix] / (b[ix] - a[ix] * scratch[ix - 1]);
			}
	x[ix] = (x[ix] - a[ix] * x[ix - 1]) / (b[ix] - a[ix] * scratch[ix - 1]);		x[ix] = (x[ix] - a[ix] * x[ix - 1]) / (b[ix] - a[ix] * scratch[ix - 1]);
	}		}

2602:61:71B4:1202:4CD9:9E58:ADE6:DE5A: /* Variants */ replace C-like pseudocode that counts from 1 with C that counts from zero, and does not divide in the back substitution

2024-01-19T23:30:59Z

Variants: replace C-like pseudocode that counts from 1 with C that counts from zero, and does not divide in the back substitution

@@ Line 336: / Line 336: @@
 <math display="block"> x_1 = \frac{d_1-a_1u_n-c_1u_2}{b_1+a_1v_n+c_1v_2}</math>
 <syntaxhighlight lang="c">
-	double u[n+1] = {}, v[n+1] = {};
+    cmod[1] = c[1] / b[1];
-	double w, F[n+1] = {};
+    v[1] = -a[1] / b[1];
-	F[2] = -A[2];
+    x[1] = x[1] / b[1];
-	for (i = 3; i < n + 1; i++) {
+        cmod[ix] = c[ix] * m;
-		w = A[i] / B[i - 1];
+        v[ix] = (0.0f  - a[ix] * v[ix - 1]) * m;
-		B[i] = B[i] - w * C[i - 1];
+        x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
-	v[n] = F[n] / B[n];
+    cmod[X - 1] = c[X - 1] * m;
-	for (i = n - 1; i > 1; i--) {
+    v[X - 1] = (-c[0]    - a[X - 1] * v[X - 2]) * m;
-		u[i] = (D[i] - C[i] * u[i + 1]) / B[i];
+    x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;
-	for (i = 2; i < n + 1; i++) {
+        x[ix] -= cmod[ix] * x[ix + 1];
+}
 </syntaxhighlight>

2602:61:71B4:1202:4CD9:9E58:ADE6:DE5A: replaced bitrotted comment

2024-01-19T21:42:08Z

replaced bitrotted comment

← Previous revision		Revision as of 21:42, 19 January 2024
Line 252:		Line 252:


	The implementation in [[C (programming language)\|C]] ~~without~~ ~~preserving~~ ~~the~~ ~~coefficient~~ ~~vectors~~:		The implementation as a [[C (programming language)\|C]] function, which uses scratch space to avoid modifying its inputs for a-c, allowing them to be reused:
	<syntaxhighlight lang="c">		<syntaxhighlight lang="c">
	void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double u[restrict X]) {		void cyclic_thomas(const int X, double x[restrict X], const double a[restrict X], const double b[restrict X], const double c[restrict X], double cmod[restrict X], double u[restrict X]) {

2602:61:71B4:1202:4CD9:9E58:ADE6:DE5A: /* Variants */ replace C-like pseudocode with actual C that indexes from zero, and does not divide in the back substitution

2024-01-19T21:40:04Z

Variants: replace C-like pseudocode with actual C that indexes from zero, and does not divide in the back substitution

@@ Line 254: / Line 254: @@
 The implementation in [[C (programming language)|C]] without preserving the coefficient vectors:
 <syntaxhighlight lang="c">
-	double u[n+1] = {}, v[n+1] = {};
+    const double beta = c[X - 1];
-	double w;
+    const double gamma = -b[0];
-	u[1] = gamma;
+    u[0] = gamma / (b[0] - gamma);
-	u[n] = C[n];
+    x[0] /= (b[0] - gamma);
-	B[1] = B[1] - gamma;
+        x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m;
-		w = A[i] / B[i - 1];
+    u[X - 1] = (alpha    - a[X - 1] * u[X - 2]) * m;
-		B[i] = B[i] - w * C[i - 1];
+    x[X - 1] = (x[X - 1] - a[X - 1] * x[X - 2]) * m;
-	y[n] = D[n] / B[n];
+        u[ix] -= cmod[ix] * u[ix + 1];
-	q[n] = u[n] / B[n];
+        x[ix] -= cmod[ix] * x[ix + 1];
-		y[i] = (D[i] - C[i] * y[i + 1]) / B[i];
+    const double fact = (x[0] + x[X - 1] * beta / gamma) / (1.0 + u[0] + u[X - 1] * beta / gamma);
+}
 </syntaxhighlight>

2602:61:71B4:1202:5524:CE74:52A9:84C: /* Method */ replace VB script example with a cleaned up C based example from the article history, without the expensive division in the back substitution

2024-01-19T19:11:49Z

Method: replace VB script example with a cleaned up C based example from the article history, without the expensive division in the back substitution

@@ Line 70: / Line 70: @@
 : <math>x_i = \cfrac{d_i - c_i x_{i+1}}{b_i} \quad \text{for } i = n - 1, n - 2, \dots, 1.</math>
-<syntaxhighlight lang="vbscript">
+<syntaxhighlight lang="c">
-        D(i) = D(i) - W * D(i - 1)
+     X = number of equations
-    For i = N - 1 To 1 Step -1
+     b[] = main diagonal, indexed from [0, ..., X - 1]
-        X(i) = (D(i) - C(i) * X(i + 1)) / B(i)
+     c[] = superdiagonal, indexed from [0, ..., X - 2]
 </syntaxhighlight>

2602:61:71B4:1202:5524:CE74:52A9:84C: /* Variants */ whitespace for readability

2024-01-19T16:29:25Z

Variants: whitespace for readability

← Previous revision		Revision as of 16:29, 19 January 2024
Line 243:		Line 243:
	//Apply Thomas Alg., unknowns x[1],...,x[n]		//Apply Thomas Alg., unknowns x[1],...,x[n]
	void ThomasAlg(int n, double x[n+1], double A[n+2], double B[n+2], double C[n+2], double D[n+2]) {		void ThomasAlg(int n, double x[n+1], double A[n+2], double B[n+2], double C[n+2], double D[n+2]) {
	double u[n+1]={},v[n+1]={};		double u[n+1] = {}, v[n+1] = {};
	double y[n+1]={},q[n+1]={};		double y[n+1] = {}, q[n+1] = {};
	double value=0;		double value = 0;
	double w;		double w;
	int i;		int i;

	u[1]=gamma;		u[1] = gamma;
	u[n]=C[n];		u[n] = C[n];

	v[1]=1;		v[1] = 1.0;
	v[n]=A[1]/gamma;		v[n] = A[1] / gamma;

	//create matrix B		//create matrix B
	A[1]=0;		A[1] = 0;
	B[1]=B[1]-gamma;		B[1] = B[1] - gamma;
	B[n]=B[n]-(C[n]*A[n])/gamma;		B[n] = B[n] - (C[n] * A[n]) / gamma;
	C[n]=0;		C[n] = 0;

	for(i=2;i<n+1;i++){		for (i = 2; i < n + 1; i++) {
	w=A[i]/B[i-1];		w = A[i] / B[i - 1];
	B[i]=B[i]-w*C[i-1];		B[i] = B[i] - w * C[i - 1];
	D[i]=D[i]-w*D[i-1];		D[i] = D[i] - w * D[i - 1];
	u[i]=u[i]-w*u[i-1];		u[i] = u[i] - w * u[i - 1];
	}		}
	y[n]=D[n]/B[n];		y[n] = D[n] / B[n];
	q[n]=u[n]/B[n];		q[n] = u[n] / B[n];

	for(i=n-1;i>0;i--){		for (i = n - 1; i > 0; i--) {
	y[i]=(D[i]-C[i]*y[i+1])/B[i];		y[i] = (D[i] - C[i] * y[i + 1]) / B[i];
	q[i]=(u[i]-C[i]*q[i+1])/B[i];		q[i] = (u[i] - C[i] * q[i + 1]) / B[i];
	}		}
	value=(v[1]y[1]+v[n]y[n])/(1+v[1]q[1]+v[n]q[n]);		value = (v[1] * y[1] + v[n] * y[n]) / (1.0 + v[1] * q[1] + v[n] * q[n]);

	for(i=1;i<n+1;i++){		for (i = 1; i < n + 1; i++) {
	x[i]=y[i]-q[i]*value;		x[i] = y[i] - q[i] * value;
	}		}
	}		}
Line 319:		Line 319:
	//Apply Thomas Alg., unknowns x[1],...,x[n]		//Apply Thomas Alg., unknowns x[1],...,x[n]
	void ThomasAlg(int n, double x[n+1], double A[n+2], double B[n+2], double C[n+2], double D[n+2]) {		void ThomasAlg(int n, double x[n+1], double A[n+2], double B[n+2], double C[n+2], double D[n+2]) {
	double u[n+1]={},v[n+1]={};		double u[n+1] = {}, v[n+1] = {};
	double w,F[n+1]={};		double w, F[n+1] = {};
	F[2]=-A[2];		F[2] = -A[2];
	F[n]=-C[n];		F[n] = -C[n];
	int i;		int i;
	u[1]=0;		u[1] = 0;
	v[1]=1;		v[1] = 1;
	for(i=3;i<n+1;i++){		for (i = 3; i < n + 1; i++) {
	w=A[i]/B[i-1];		w = A[i] / B[i - 1];
	B[i]=B[i]-w*C[i-1];		B[i] = B[i] - w * C[i - 1];
	D[i]=D[i]-w*D[i-1];		D[i] = D[i] - w * D[i - 1];
	F[i]=F[i]-w*F[i-1];		F[i] = F[i] - w * F[i - 1];
	}		}
	u[n]=D[n]/B[n];		u[n] = D[n] / B[n];
	v[n]=F[n]/B[n];		v[n] = F[n] / B[n];
	for(i=n-1;i>1;i--){		for (i = n - 1; i > 1; i--) {
	u[i]=(D[i]-C[i]*u[i+1])/B[i];		u[i] = (D[i] - C[i] * u[i + 1]) / B[i];
	v[i]=(F[i]-C[i]*v[i+1])/B[i];		v[i] = (F[i] - C[i] * v[i + 1]) / B[i];
	}		}
	x[1]=(D[1]-A[1]u[n]-C[1]u[2])/(B[1]+A[1]v[n]+C[1]v[2]);		x[1] = (D[1] - A[1] * u[n] - C[1] * u[2]) / (B[1] + A[1] * v[n] + C[1] * v[2]);

	for(i=2;i<n+1;i++){		for (i = 2; i < n + 1; i++) {
	x[i]=u[i]+x[1]*v[i];		x[i] = u[i] + x[1] * v[i];
	}		}
	}		}