Routine Name: gaussSeidelSolver
Author: Brandon Furman
Language: C++
Description/Purpose: This function solves the square linear system of equations Ax = b using the Gauss-Seidel algorithm. The Gauss-Seidel algorithm is an iterative algorithm that is guaranteed to converge if the coefficient matrix is diagonally dominant or symmetric positive-definite.
Input: This function requires the following 5 items as inputs:
Output: This function returns an array1D object that is the solution to the linear system of equations. If the algorithm is unable to converge to a solution within the specified number of iterations, an exception is thrown.
Usage/Example: Usage of this function is straightforward. Consider the problem of solving the linear system with A = [[7,3,1],[-3,10,2],[1,7,-15]] and b = [3,4,2]. The following code will solve the given system:
int m = 3;
int maxIter = 100;
double tol = 0.0001;
//Create and populate the coefficient matrix.
array2D A;
A.allocateMem(m, m);
A(0, 0) = 7; A(0, 1) = 3; A(0, 2) = 1;
A(1, 0) = -3; A(1, 1) = 10; A(1, 2) = 2;
A(2, 0) = 1; A(2, 1) = 7; A(2, 2) = -15;
//Create and populate the vector of constant terms.
array1D b;
b.allocateMem(m);
b(0) = 3; b(1) = 4; b(2) = 2;
//Create and populate the initial guess.
array1D x;
x.allocateMem(m);
x(0) = 1.5; x(1) = -1.5; x(2) = 1;
//Use Jacobi Iteration to solve the system Ax = b.
x = gaussSeidelSolver(A, b, x, tol, maxIter);
//Print the results.
for (int i = 0; i < 3; i++) {
std::cout << x(i) << " ";
}
This code outputs the following to the console:
0.223222 0.448796 0.0910068
which is the correct solution to the linear system. This function can also be used to solve must larger systems of equations. The following code generates a random coefficient 1000 x 1000 matrix and solves it. The answer is then verified to be correct by checking the L2 norm of the residual.
int m = 1000;
int maxIter = 10000;
double tol = 0.0001;
//Create and populate a random
//diagonally dominant 1000 x 1000 matrix.
array2D A;
A = randDiagDomMat(m);
//Create and populate a random vector
//of length 1000
array1D b;
b = randVec(m);
//Create and populate an initial guess.
array1D x;
x = randVec(m);
//Solve the system of equations Ax = b for x.
x = gaussSeidelSolver(A, b, x, tol, maxIter);
//Calculate the residue: r = b - Ax.
array1D res;
res = multMatVec(A, x);
res = subVec(b, res);
//Calculate the L2 norm of the residue.
double error;
error = twoNormVec(res);
//Print the error in the solution.
std::cout << error << std::endl;
which outputs the following to the console:
0.00395319
which shows that the function did indeed find a correct solution. If the function is unable to converge to an answer within the number of iterations specified by maxIter, then an exception will be thrown. This exception can be caught with a try-catch block and results in the following message:
gaussSeidelSolver: Failed to converge
Implementation/Code: The gaussSeidelSolver() function is implemented as follows:
array1D gaussSeidelSolver(array2D& A, array1D& b, array1D x0, double tol, int maxIter) {
int m = A.getRows();
int n = A.getCols();
int p = b.getLength();
if (m != n || m != p) {
throw "gaussSeidelSolver: Incompatible Matrix Sizes";
}
array1D x1;
x1.allocateMem(p);
x1 = x0;
int cntr = 0;
double error = 10.0*tol;
double sum = 0.0;
while (error > tol && cntr < maxIter) {
cntr += 1;
cntr2 = cntr2 + 1;
for (int i = 0; i < m; i++) {
sum = 0.0;
for (int j = 0; j < n; j++) {
if (i != j) {
sum = sum + A(i, j)*x1(j);
}
}
x1(i) = (1.0 / A(i, i))*(b(i) - sum);
}
error = absErrVecTwoNorm(x1, x0);
x0 = x1;
}
return x1;
}
Last Modified: April/2019