Numerical Methods: Mathematica 3.0 Notebooks
(c) John H. Mathews, 1998
To accompany the text:
Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Simon & Schuster,  One Lake Street
Upper Saddle River, New Jersey 07458  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5
www.Prenhall.com
This free software is compliments of the author.
E-mail address:      in%"mathews@fullerton.edu"

CONTENTS

Chapter 3.   The Solution of Linear Systems  AX = B
	Algorithm 3.1  Back Substitution.
	Algorithm 3.2  Upper-Triangularization Followed by Back Substitution.
	Algorithm 3.3  PA = LU Factorization with Pivoting.
	Algorithm 3.4  Jacobi Iteration.
	Algorithm 3.5  Gauss-Seidel Iteration.

Clear["Global`*"];
Off[General::"spell1"];
Clear[a, a0, A, b, b0, B, BackSub, c, GaussJordan, GaussSeidel, i, j, 
        Jacobi, k, lu0, LUfactor, m, ma, max, mna, n, n1, na, nb, 
        Nonsingular, p, P0, x, row, rowk, rowp, sum, Singular, SolveLU, Tol, 
        V, x];
On[General::"spell1"];

Algorithm 3.1 Back Substitution. To solve the upper-triangular system:
          (Subscript[a, 1, 1]) (Subscript[x, 1]) + (Subscript[a, 1, 2]) (Subscript[x, 2]) + (Subscript[a, 1, 3]) (Subscript[x, 3]) + 
          ⋯ + (Subscript[a, 1, n - 1]) (Subscript[x, n - 1]) + 
          (Subscript[a, 1, n]) (Subscript[x, n])  =  (Subscript[b, 1])                 
      (Subscript[a, 2, 2]) (Subscript[x, 2]) + (Subscript[a, 2, 3]) (Subscript[x, 3]) + ⋯ + 
          (Subscript[a, 2, n - 1]) (Subscript[x, n - 1]) + (Subscript[a, 2, n]) (Subscript[x, n])  =  
        (Subscript[b, 2])                                  
      (Subscript[a, 3, 3]) (Subscript[x, 3]) + ⋯ + 
          (Subscript[a, 3, n - 1]) (Subscript[x, n - 1]) + (Subscript[a, 3, n]) (Subscript[x, n])  =  
        (Subscript[b, 3])\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
      ⋮\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
        ⋮\ \ \ \ \ \ \ \ \ ⋮\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ 
      \(a\_\(n - 1, n - 1\)\) x\_\(n - 1\) + \(a\_\(n - 1, n\)\) x\_n\  = \ 
        b\_\(n - 1\)\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(a\_\(n, n\)\) x\_n\  = \ b\_n\)
Proceed with the method only if all the diagonal elements are non-zero.   First compute  (Subscript[x, n]) = (((Subscript[b, n]))/((Subscript[a, n, n])))  and then use the rule
(Subscript[x, k]) = (((Subscript[b, k]) - (Subscript[a, k, k + 1]) (Subscript[x, k + 1]) - ⋯ - 
              (Subscript[a, k, n - 1]) (Subscript[x, n - 1]) - (Subscript[a, k, n]) (Subscript[x, n]))/((Subscript[a, k, k])))  for  k = n-1, n-2, ⋯ , 1.
Remark.  Additional steps have been placed in the algorithm which perform the side calculation of  det(A),  the determinant of  A. 
Section  3.3   Upper-Triangular Linear Systems   Page  145

Execute this cell to activate Back Substitution.
BackSub[a0_,b0_]
a0_ is the matrix A,
b0_ is the column vector B.

Off[General::"spell1"]; Clear[BackSub, a0, b0]; \(BackSub[a0_, b0_] := 
      Module[{a, b, x, k, sum, mna, nb, ma, na}, 
\ \ \ a = a0; 
\ \ \ 
        b = b0; 
\ \ \ mna = Dimensions[a]; 
\ \ \ nb = Length[b]; 
\ \ \ 
        ma = mna[[1]]; 
\ \ \ 
        na = mna[[2]]; 
\ \ \ det = 1; 
        
\ \ \ For[k = 1, k ≤ ma, \(k++\), 
          \(det = det\ a[[k, k]]; \)]; 

        \ \ \ \ \ \ 
        If[ma == na && det ≠ 0, 
\ \ \ \ \ \ 
          Module[{i, j}, 
\ \ \ \ \ \ \ \ \ x = Table[0, {i, ma}]; 

            \ \ \ \ \ \ \ \ \ 
            x[[ma]] = 
              b[[ma]]\/a
                  [[ma, ma]]; 

            \ \ \ \ \ \ \ \ \ Clear[i, j, k]; 
\ \ \ \ \ \ \ \ \ 
            For[i = ma - 1, 1 ≤ i, \(i--\), 
\ \ \ \ \ \ \ \ \ 
              sum = 0; 
\ \ \ \ \ \ \ \ \ 
              For[j = i + 1, j ≤ ma, \(j++\), 
                sum = sum + 
                    a[[i, j]]\ 
                      x[[j]]]; 

              \ \ \ \ \ \ \ \ \ 
              x[[i]] = 
                \(b[[i]] - sum\)\/a
                    [[i, i]]; ]; ], 

          \ \ \ \ \ \ \ \ \ x = MatrixIsSingular]; 
\ \ \ Return[x]; ]; 

    \(On[General::"spell1"]; \)\)

Chapter 3, Example 3.12, Page 144.
Use back substitution to solve the upper-triangular linear system  AX = B:

A = (({{3, -1, 2, 3}, {0, -2, 7, -4}, {0, 0, 6, 5}, {0, 0, 0, 3}}));
B = (({{20}, {-7}, {4}, {6}}));
Xs = (({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}}));
X = BackSub[A, B];
Print["Solve ", MatrixForm[A], MatrixForm[Xs], " = ", 
        MatrixForm[B]];
Print["X = ", MatrixForm[X]];

Verify by matrix multiplication to see if  AX  =  B.

Print[MatrixForm[A], MatrixForm[X], " = ", MatrixForm[A.X], 
      ",  B = ", MatrixForm[B]];

Compare with Mathematica's built in routine.

\(Z = LinearSolve[A, B]; 
Print["Z = ", MatrixForm[Z]]; 

    Print[MatrixForm[A], MatrixForm[Z], " = ", MatrixForm[A.Z], 
      ",  B = ", MatrixForm[B]]; \)

Algorithm 3.2 Upper-Triangularization Followed by Back Substitution.
To construct the solution to  AX=B,  by first reducing the augmentedmatrix  [A,B]  to upper triangular form and then performing back substitution.
Section  3.4   Gaussian Elimination and Pivoting   Page  156

Let the coefficients of  A  and the constants of  B  be stored in the augmented 
matrix  [A, B] = ((Subscript[a, i, j]))  which has n rows and  n+1  columns; that is, the column 
vector  B = ((Subscript[a, i, n + 1]))  is stored in column  n+1  of the augmented matrix  [A, B]. 
The system can be written
        
      (Subscript[a, 1, 1]) (Subscript[x, 1])   +     (Subscript[a, 1, 2]) (Subscript[x, 2])   + 
          ⋯ +      (Subscript[a, 1, n - 1]) (Subscript[x, n - 1])  + 
             (Subscript[a, 1, n]) (Subscript[x, n])    =  (Subscript[a, 1, n + 1])   
      (Subscript[a, 2, 1]) (Subscript[x, 1])   +     (Subscript[a, 2, 2]) (Subscript[x, 2 ])  + 
          ⋯ +      (Subscript[a, 2, n - 1]) (Subscript[x, n - 1])  + 
             (Subscript[a, 2, n]) (Subscript[x, n])    =  (Subscript[a, 2, n + 1])(Subscript[a, n - 1, 1]) (Subscript[x, 1])   +  (Subscript[a, n - 1, 2]) (Subscript[x, 2])  + 
          ⋯ +  (Subscript[a, n - 1, n - 1]) (Subscript[x, n - 1])   +  
          (Subscript[a, n - 1, n]) (Subscript[x, n])  =  (Subscript[a, n - 1, n + 1])   
      (Subscript[a, n, 1]) (Subscript[x, 1])   +    (Subscript[a, n, 2]) (Subscript[x, 2])    + 
          ⋯ +     (Subscript[a, n, n - 1]) (Subscript[x, n - 1])   + 
             (Subscript[a, n, n]) (Subscript[x, n])   =  (Subscript[a, n, n + 1])
Row operations will be used to eliminate  (Subscript[x, p])  in column  p.
REMARKS. Rather than actually interchange rows, row interchanges are recorded in the pointer vector Row(p).  If either the matrix A or vector B is to be used later, then a copy must be stored someplace else besides the augmented matrix.

Execute this cell to activate Upper-Triangularization Followed by Back Substitution.
GaussJordan[a0_]
a0_ is the augmented matrix [A,B].

Off[General::"spell1"]; Clear[GaussJordan, a0]; \(GaussJordan[a0_] := 
      Module[{a, i, j, mna, n, n1, p, row, x, cond, Nonsingular, Singular}, 

        		a = a0; 
		mna = Dimensions[a]; 
		
        n = mna[[1]]; 
		
        n1 = mna[[2]]; 
		
        x = Table[0, {i, n}]; 
		cond = Nonsingular; 
		Clear[row]; 

        		For[j = 1, j ≤ n, \(j++\), row[j] = j]; 
		
        For[p = 1, p ≤ n - 1, \(p++\), 
			
          \(Block[{k}, 
				
            For[k = p + 1, k ≤ n, \(k++\), 
					
              \(If[Abs[a[[row[k], p]]] > 
                  Abs[a[[row[p], p]]], 
                t = row[p]; row[p] = row[k]; row[k] = t]; \)]; 
				
            If[a[[row[p], p]] == 0, 

              					
              Print["A(", row[p], ",", p, ") = ", 
                a[[row[p], p]]]; 

              					cond = Singular; 
					
              Print["The matrix is singular."]; 
					
              Goto[getout]; ]; 
				Clear[k]; 
				
            For[k = p + 1, k ≤ n, \(k++\), 
					
              \(Module[{c}, 
                m = a[[row[k], p]]\/a
                      [[row[p], p]]; 

                						
                For[c = p + 1, c ≤ n + 1, \(c++\), 
							
                  a[[row[k], c]] = 
                    a[[row[k], c]] - 
                      m\ a[[row[p], c
                          ]]]; ]; \)]; ]; \)]; 
		
        If[a[[row[n], n]] == 0, 
			
          Print["A(", row[n], ",", n, ") = ", 
            a[[row[n], n]]]; 
			
          cond = Singular; 
			Print["The matrix is singular."]; 

          			Goto[getout]; ]; 
		
        x[[n]] = 
          a[[row[n], n + 1]]\/a
              [[row[n], n]]; 
		
        Clear[k]; 
		
        For[k = n - 1, 1 ≤ k, \(k--\), 
			sum := 0; 
			
          For[c = k + 1, c ≤ n, \(c++\), 
				
            sum = sum + 
                a[[row[k], c]]\ 
                  x[[c]]; 
				
            x[[k]] = 
              \(a[[row[k], n + 1]] - sum
                  \)\/a[[row[k], k
                  ]]; ]; ]; 
		Label[getout]; 
		
        Return[x]; ]; 
\(On[General::"spell1"]; \)\)

Chapter 3, Exercise 8, Page 158.  Use upper-triangularization followed by 
back-substitution to solve the linear system  AX = B:

A = (({{4, 8, 4, 0}, {1, 5, 4, -3}, {1, 4, 7, 2}, {1, 3, 0, -2}}));
B = (({{8}, {-4}, {10}, {-4}}));
Xs = (({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}, {(Subscript[x, 4])}}));
AB = (({{4, 8, 4, 0, 8}, {1, 5, 4, -3, -4}, {1, 4, 7, 2, 10}, {1, 3, 0, -2, -4}}));
X = GaussJordan[AB];
Print["Solve ", MatrixForm[A], MatrixForm[Xs], " = ", 
        MatrixForm[B]];
Print["X = ", MatrixForm[X]];

Verify by matrix multiplication to see if  AX  =  B.

Print[MatrixForm[A], MatrixForm[X], " = ", MatrixForm[A.X], 
      ",  B = ", MatrixForm[B]];

Compare with Mathematica's built in routine.

\(Z = LinearSolve[A, B]; 
Print["Z = ", MatrixForm[Z]]; 

    Print[MatrixForm[A], MatrixForm[Z], " = ", MatrixForm[A.Z], 
      ",  B = ", MatrixForm[B]]; \)

Algorithm 3.3 PA = LU Factorization with Pivoting.
To construct the solution to the linear system  AX=B  by performing the steps:
1.  Find a permutation matrix  P , lower-triangular matrix  L   and 
     upper-triangular matrix  U  which satisfy  PA=LU.
2.  Compute  PB  and form the equivalent linear system  LUX = PB.
3.  Solve the lower-triangular system  LY = PB  for  Y.
4.  Solve the upper-triangular system  UX = Y  for  X.
Remarks: Since the diagonal elements of L are all 1's these values do not need to be stored.  The coefficients of  L  below the main diagonal and the nonzero coefficients of  U  overwrite the matrix  A. The permutation matrix  P  is not actually needed.  The pointer vector  row serves in its place to keep track of row interchanges.  If matrix  A  is to be used later, then a copy must be stored someplace else.  The value  det  is the determinant of A.
Section  3.6   Triangular Factorizaton   Page  175

Execute this cell to activate PA = LU Factorization with Pivoting.
LUfactor[a0_]
SolveLU[lu0_,b0_]
a0_ is the matrix A,
b0_ is the column vector B,
lu0_ is the factored matrix.

Off[General::"spell1"]; Clear[LUfactor, a0, row]; \(\(LUfactor[a0_] := 
      Module[{a, i, j, p, mna, n, m, rowk, rowp}, 
		a = a0; 
		
        mna = Dimensions[a]; 
		
        n = mna[[1]]; 
		
        m = mna[[2]]; 
		det = 1; 

        		For[j = 1, j ≤ n, \(j++\), row[j] = j]; 
		
        For[p = 1, p ≤ n - 1, \(p++\), 
			
          \(Block[{c, k}, 
				
            For[k = p + 1, k ≤ n, \(k++\), 
					
              \(If[Abs[a[[row[k], p]]] > 
                  Abs[a[[row[p], p]]], 
                t = row[p]; 
						row[p] = row[k]; 
						
                row[k] = t; 
						det = \(-det\)]; \)]; 
				
            det = det\ a[[row[p], p]]; 

            				If[det == 0, Print["The matrix is singular."]; 

              					Goto[getout]; ]; 
				Clear[k]; 
				
            For[k = p + 1, k ≤ n, \(k++\), 
					
              rowk = row[k]; 
					rowp = row[p]; 
					
              a[[rowk, p]] = 
                a[[rowk, p]]\/a
                    [[rowp, p]]; 

              					
              For[c = p + 1, c ≤ n, \(c++\), 
						
                \(a[[rowk, c]] = 
                  a[[rowk, c]] - 
                    a[[rowk, p]]\ 
                      a[[rowp, c]]; 
                \)]; ]; ]; \)]; 
		Label[getout]; 
		
        det = det\ a[[row[n], n]]; 
		
        Return[a]; ]; \)\)

Clear[SolveLU, lu0, b0]; \(SolveLU[lu0_, b0_] := 
	
      Module[{a, c, i, j, k, n, p, x, rowk, sum}, 
		a = lu0; 
		
        b = b0; 
		n = Length[a]; 
		x = Table[0, {i, n}]; 
		
        For[k = 1, k ≤ n, \(k++\), 
			
          \(If[a[[row[k], k]] == 0, 
            Goto[getout]]; \)]; 
		
        x[[1]] = 
          b[[row[1]]]; 
		Clear[k]; 

        		For[k = 2, k ≤ n, \(k++\), 
			sum = 0; 
			
          rowk = row[k]; 
			
          For[c = 1, c ≤ k - 1, \(c++\), 
				
            sum = sum + 
                a[[rowk, c]]\ 
                  x[[c]]]; 
			
          x[[k]] = 
            b[[rowk]] - sum; ]; 
		
        x[[n]] = 
          x[[n]]\/a
              [[row[n], n]]; 
		
        Clear[k]; 
		
        For[k = n - 1, 1 ≤ k, \(k--\), 
			sum = 0; 
			
          rowk = row[k]; 
			
          For[c = k + 1, c ≤ n, \(c++\), 
				
            sum = sum + 
                a[[rowk, c]]\ 
                  x[[c]]]; 
			
          x[[k]] = 
            \(x[[k]] - sum\)\/a
                [[rowk, k]]; ]; 
		
        Label[getout]; 
		Return[x]; ]; 
\(On[General::"spell1"]; 
    \)\)

Chapter 3, Exercise 7, Page 178.
Use LU factorization with pivoting to solve the linear system  AX = B:

A = (({{2, 4, -4, 0}, {1, 5, -5, -3}, {2, 3, 1, 3}, {1, 4, -2, 2}}));
LU = LUfactor[A];
Print["LU = ", MatrixForm[LU]];

B = (({{12}, {18}, {8}, {8}}));
X = SolveLU[LU, B];
Print["Solve ", MatrixForm[A], MatrixForm[Xs], " = ", 
        MatrixForm[B]];
Print["X = ", MatrixForm[X]];

Verify by matrix multiplication to see if  AX  =  B.

Print[MatrixForm[A], MatrixForm[X], " = ", MatrixForm[A.X], 
      ",  B = ", MatrixForm[B]];

Compare with Mathematica's built in routine.

\(Z = LinearSolve[A, B]; 
Print["Z = ", MatrixForm[Z]]; 

    Print[MatrixForm[A], MatrixForm[Z], " = ", MatrixForm[A.Z], 
      ",  B = ", MatrixForm[B]]; \)

Algorithm 3.4 Jacobi Iteration.  To solve the linear system  AX = B  by starting with  (Subscript[P, 0])=0,  and generating a sequence  {(Subscript[P, k])}  which converges to the solution  P,  i.e.  AP = B.  A sufficient condition for the method to be applicable is that  A  is diagonally  dominant. 
Section  3.7   Iterative Methods for Linear Systems   Page  186

Execute this cell to activate Jacobi Iteration.
Jacobi[A_,B_,P0_,Tol_,max_]
A_ is the matrix A,
B_ is the column vector B,
P0_ is the starting vector,
Tol_ is the tolerance,
max_ is the maximum number of iterations.

Off[General::"spell1"]; Clear[Jacobi, A, B, V, P0, Tol, max]; \(Jacobi[A_, B_, P0_, Tol_, max_] := 
	
      Block[{r, c, k, n, cond, stat, Xold, row, DiagonallyDominant, 
          NotDominant}, 
		
        norm[V_] := N[\@ReplaceAll[V\^2, List -> Plus]]; 
		sep = 1; 
        
		k = 0; 
		cond = 1; 
		stat = DiagonallyDominant; 
		
        X = N[P0]; 
		Xold = P0; 
		n = Length[X]; 
		
        For[r = 1, r ≤ n, \(r++\), 
			row = 0; 
			
          For[c = 1, c ≤ n, \(c++\), 
				
            row = N[row + 
                  Abs[A[[r, c]]]]]; 

          			If[
            row ≥ 
              2\ Abs[A[[r, r]]], 
            cond = 0]]; 
		
        If[cond == 0, stat = NotDominant; Goto[getout]]; 
		
        While[k < max && sep > Tol, 
			
          For[r = 1, r ≤ n, \(r++\), 
            sum = N[B[[r]]]; 
				
            For[c = 1, c ≤ n, \(c++\), 
					
              \(If[c ≠ r, 
                sum = N[
                    sum - A[[r, c]]\ 
                        Xold[[c]]]]; 
              \)]; 
				
            X[[r]] = 
              N[sum\/A[[r, r]]]]; 

          			sep = norm[X - Xold]; 
			
          Print[MatrixForm[Transpose[N[X, 10]]]]; 
			Xold = X; 
          k = k + 1]; 
		Label[getout]; ]; 
On[General::"spell1"]; 
    \)

Chapter 3, Example 3.25, Page 180.
Use Jacobi iteration to solve the linear system  AX = B:

Clear[A, B, P0, X];
A =(({{4, -1, 1}, {4, -8, 1}, {-2, 1, 5}}));
B = (({{7}, {-21}, {15}}));
P0 = (({{1}, {2}, {2}}));
Xs = (({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}}));
Jacobi[A, B, P0, 0.00000001, 100];
Print["Solve ", MatrixForm[A], MatrixForm[Xs], " = ", 
        MatrixForm[B]];
Print["X = ", MatrixForm[X]];

Verify by matrix multiplication to see if  AX  =  B.

Print[MatrixForm[A], MatrixForm[X], " = ", MatrixForm[A.X], 
      ",  B = ", MatrixForm[B]];

Compare with Mathematica's built in routine.

\(Z = LinearSolve[A, B]; 
Print["Z = ", MatrixForm[Z]]; 

    Print[MatrixForm[A], MatrixForm[Z], " = ", MatrixForm[A.Z], 
      ",  B = ", MatrixForm[B]]; \)

Algorithm 3.5 Gauss-Seidel Iteration.  To solve the linear system  AX = B  by starting with  (Subscript[P, 0])=0,  and generating a sequence  {(Subscript[P, k])}  which converges to the solution  P,  i.e.  AP = B.  A sufficient condition for the method to be applicable is that  A  is diagonally dominant. 
Section  3.7   Iterative Methods for Linear Systems   Page  187

Execute this cell to activate Gauss-Seidel Iteration.
GaussSeidel[A_,B_,P0_,Tol_,max_]
A_ is the matrix A,
B_ is the column vector B,
P0_ is the starting vector,
Tol_ is the tolerance,
max_ is the maximum number of iterations.

Off[General::"spell1"]; Clear[GaussSeidel, A, B, V, P0, Tol, max]; \(GaussSeidel[A_, B_, P0_, Tol_, max_] := 
	
      Block[{r, c, k, n, Xold, cond, row, stat, DiagonallyDominant, 
          NotDominant}, 
		
        norm[V_] := N[\@ReplaceAll[V\^2, List -> Plus]]; 
		sep = 1; 
        
		k = 0; 
		cond = 1; 
		stat = DiagonallyDominant; 
		
        X = N[P0]; 
		Xold = P0; 
		n = Length[X]; 
		
        For[r = 1, r ≤ n, \(r++\), 
			row = 0; 
			
          For[c = 1, c ≤ n, \(c++\), 
				
            row = N[row + 
                  Abs[A[[r, c]]]]]; 

          			If[
            row ≥ 
              2\ Abs[A[[r, r]]], 
            cond = 0]]; 
		
        If[cond == 0, stat = NotDominant; Goto[getout]]; 
		
        While[k < max && sep > Tol, 
			
          For[r = 1, r ≤ n, \(r++\), 
            sum = N[B[[r]]]; 
				
            For[c = 1, c ≤ n, \(c++\), 
					
              \(If[c ≠ r, 
                sum = N[
                    sum - A[[r, c]]\ 
                        X[[c]]]]; \)]; 

            				X[[r]] = 
              N[sum\/A[[r, r]]]]; 

          			sep = norm[X - Xold]; 
			
          Print[MatrixForm[Transpose[N[X, 10]]]]; 
			Xold = X; 
          k = k + 1]; 
		Label[getout]; ]; 
On[General::"spell1"]; 
    \)

Chapter 3, Example 3.25, Page 180.
Use Gauss-Seidel iteration to solve the linear system:

Clear[A, B, P0, X];
A =(({{4, -1, 1}, {4, -8, 1}, {-2, 1, 5}}));
B = (({{7}, {-21}, {15}}));
P0 = (({{1}, {2}, {2}}));
Xs = (({{(Subscript[x, 1])}, {(Subscript[x, 2])}, {(Subscript[x, 3])}}));
GaussSeidel[A, B, P0, 0.00000001, 100]\(Print["Solve ", MatrixForm[A], MatrixForm[Xs], " = ", 
      MatrixForm[B]]; 
Print["X = ", MatrixForm[X]]; \)

Verify by matrix multiplication to see if  AX  =  B.

Print[MatrixForm[A], MatrixForm[X], " = ", MatrixForm[A.X], 
      ",  B = ", MatrixForm[B]];

Compare with Mathematica's built in routine.

\(Z = LinearSolve[A, B]; 
Print["Z = ", MatrixForm[Z]]; 

    Print[MatrixForm[A], MatrixForm[Z], " = ", MatrixForm[A.Z], 
      ",  B = ", MatrixForm[B]]; \)