Direct Solutions to Linear Systems

Numerical Methods

David Mayerich

Scalable Tissue Imaging and Modeling (STIM) Laboratory

Department of Electrical and Computer Engineering

Cullen College of Engineering

University of Houston

David Mayerich

STIM Laboratory, University of Houston

Partial Pivoting

Roundoff in Gaussian Elimination

Partial Pivoting

Pivoting Strategies

David Mayerich

STIM Laboratory, University of Houston

Roundoff in Gaussian Elimination

  • What is the exact solution to the following linear system using 3 digits of precision?

David Mayerich

STIM Laboratory, University of Houston

\begin{array}{llllll} 1.23 \times 10^{-4} x &+&y&=&1\\ 1.00\times 10^0 x &+& y &=& 1 \end{array}
\begin{bmatrix} 0.000123 & 1.00\\ 1.00 & 1.00 \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} 1.00\\ 2.00 \end{bmatrix}
\left[ \begin{array}{cc|c} 0.000123 & 1.00 & 1.00 \\ 1.00 & 1.00 & 2.00 \end{array} \right]
r_2 \leftarrow r_2 - \frac{1}{0.000123}r_1
r_2 \leftarrow r_2 - 8130r_1
\left[ \begin{array}{cc|c} 0.000123 & 1.00 & 1.00 \\ 0.00 & -8130 & -8130 \end{array} \right]
-8130y=-8130
y=1.00
0.000123x+1.00(1.00)=1.00
0.000123x=0.00
x=0.00
1.00 - 8130 \rightarrow -8130
2.00 - 8130 \rightarrow -8130

Swapping Pivot Equations

  • What is the exact solution to the following linear system using 3 digits of precision?

David Mayerich

STIM Laboratory, University of Houston

\begin{array}{llllll} 1.23 \times 10^{-4} x &+&y&=&1\\ 1.00\times 10^0 x &+& y &=& 1 \end{array}
\begin{bmatrix} 0.000123 & 1.00\\ 1.00 & 1.00 \end{bmatrix} \begin{bmatrix} x\\y \end{bmatrix} = \begin{bmatrix} 1.00\\ 2.00 \end{bmatrix}
r_2 \leftarrow r_2 - {0.000123}r_1
1.00 y=1.00
y=1.00
1.00 x+1.00(1.00)=2.00
1.00x=1.00
x=1.00
\left[ \begin{array}{cc|c} 1.00 & 1.00 & 2.00\\ 0.00 & 1.00 & 1.00 \\ \end{array} \right]
\left[ \begin{array}{cc|c} 1.00 & 1.00 & 2.00\\ 0.000123 & 1.00 & 1.00 \\ \end{array} \right]
1.00 - 0.000123 \rightarrow 1.00
1.00 - 0.000123 \rightarrow 1.00

Partial Pivoting

  • The pivot row is scaled to eliminate leading coefficients in other rows

David Mayerich

STIM Laboratory, University of Houston

\begin{bmatrix} 0.10 & 300 & -2.00\\ 2.00 & 10.0 & -1.00\\ 3.00 & 6.00 & 12.0 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix} 7.00\\ 10.0\\ 15.0 \end{bmatrix}

  • Select the pivot equation with the largest leading coefficient and swap rows

  • Perform row reduction

r_2 \leftarrow r_2 - \frac{2.00}{3.00}r_1
r_2 \leftarrow r_2 - \frac{0.10}{3.00}r_1
\left[ \begin{array}{ccc|c} 3.00 & 6.00 & 12.00 & 15.0\\ 0.00 & 6.00 & -9.00 & 0.00\\ 0.00 & 300 & -2.40 & 6.50 \end{array} \right]
\left[ \begin{array}{ccc|c} 3.00 & 6.00 & 12.00 & 15.0\\ 2.00 & 10.0 & -1.00 & 10.0\\ 0.10 & 300 & -2.00 & 7.00 \end{array} \right]

Partial Pivoting (iteration 2)

  • The pivot row is scaled to eliminate leading coefficients in other rows

David Mayerich

STIM Laboratory, University of Houston

  • Select the pivot equation with the largest leading coefficient and swap rows

  • Perform row reduction

r_2 \leftarrow r_2 - \frac{6.00}{300}r_1
\left[ \begin{array}{ccc|c} 3.00 & 6.00 & 12.00 & 15.0\\ 0.00 & 6.00 & -9.00 & 0.00\\ 0.00 & 300 & -2.40 & 6.50 \end{array} \right]
\left[ \begin{array}{ccc|c} 3.00 & 6.00 & 12.00 & 15.0\\ 0.00 & 300 & -2.40 & 6.50\\ 0.00 & 0.00 & -8.95 & -0.13 \end{array} \right]
\left[ \begin{array}{ccc|c} 3.00 & 6.00 & 12.00 & 15.0\\ 0.00 & 300 & -2.40 & 6.50\\ 0.00 & 6.00 & -9.00 & 0.00 \end{array} \right]

Comparison

David Mayerich

STIM Laboratory, University of Houston

134 x_3 = 0.00
-8.95 x_3 = -0.13
-5990 x_2 - 41.0(0.00) = -130
300 x_2 - 2.40(0.0145)=6.50
0.10x_1 - 300 (0.0217) - 2.00(0.00)=7.00
3.00x_1 - 6.00(0.0217) + 12.0(0.0145)=15.0
x_3= 0.0145
x_2= 0.0218
x_1 = 4.90
x_3 = 0.00
x_2=0.0217
x_1 = 4.9
x_3 = 0.0145
x_2 = 0.0218
x_1 = 4.9
\epsilon_r=0.0\%
\epsilon_r=0.0\%
\epsilon_r=0.0\%
\epsilon_r=100\%
\epsilon_r=0.5\%
\epsilon_r=0.0\%
\left[ \begin{array}{ccc|c} 0.10 & 300 & -2.00 & 7.00\\ 0.00 & -5990 & -41.0 & -130\\ 0.00 & 0.00 & 134 & 0.00 \end{array} \right]

Gaussian Elimination

\left[ \begin{array}{ccc|c} 3.00 & 6.00 & 12.0 & 15.0\\ 0.00 & 300 & -2.40 & 6.50\\ 0.00 & 0.00 & -8.95 & -0.13 \end{array} \right]

Partial Pivoting

Implementation Strategies

  • Swapping is slow: introduces an \(O(n)\) memory copy every iteration:

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1\\ a_{21} & a_{22} & \cdots & a_{2n} & b_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \end{array} \right]

Implementation Strategies

  • Swapping is slow: introduces an \(O(n)\) memory copy every iteration:

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1\\ a_{21} & a_{22} & \cdots & a_{2n} & b_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \end{array} \right]
\left[ \begin{array}{cccc|c} a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \end{array} \right]

Implementation Strategies

  • Swapping is slow: introduces an \(O(n)\) memory copy every iteration:

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1\\ a_{21} & a_{22} & \cdots & a_{2n} & b_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \end{array} \right]
\left[ \begin{array}{cccc|c} a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \end{array} \right]

Implementation Strategies

  • Swapping is slow: introduces an \(O(n)\) memory copy every iteration:

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1\\ a_{n1} & a_{n2} & \cdots & a_{nn} & b_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \end{array} \right]
\left[ \begin{array}{cccc|c} a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \end{array} \right]

Implementation Strategies

  • Swapping is slow: introduces an \(O(n)\) memory copy every iteration:

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1\\ a_{n1} & a_{n2} & \cdots & a_{nn} & b_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \end{array} \right]

  • Use an \(n\)-dimensional index vector \(\ell \in \mathbb{N}^n\)

\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1\\ a_{21} & a_{22} & \cdots & a_{2n} & b_2\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \end{array} \right]
\ell = \begin{bmatrix} 1\\ 2\\ \vdots\\ n \end{bmatrix}
\begin{bmatrix} 1\\ n\\ \vdots\\ 2 \end{bmatrix}
int i = 2 // row index
int j = 1 // column index

// direct access
x = C[i * n + j];  // x = a_21

// indirect access
y = C[l[i] * n + j];  // y = a_n1

Partial Pivoting Algorithm

David Mayerich

STIM Laboratory, University of Houston

void gaussian(double* A, double* b, int N){
  for(int k=0; k<N-1; k++){ 			// for each pivot
    for(int i=k+1; i<N; i++){ 			// for each equation
      float m=A[i*N+k] / A[k*N+k];		// calculate scale
      for(int j=k; j<N; j++){ 			// for each coefficient
        A[i*N+j] -= m * A[k*N+j];		// subtract
      }
      b[i] -= m * b[k];					// forward elimination
    }
  }
}

Partial Pivoting Algorithm

David Mayerich

STIM Laboratory, University of Houston

void ppivot(double* A, double* b, int N){
  int* L;
  init(L);								// initialize the index vector
  for(int k=0; k<N-1; k++){ 			// for each pivot
    for(int i=k+1; i<N; i++){ 			// for each equation
      float m=A[i*N+k] / A[k*N+k];		// calculate scale
      for(int j=k; j<N; j++){ 			// for each coefficient
        A[i*N+j] -= m * A[k*N+j];		// subtract
      }
      b[i] -= m * b[k];					// forward elimination
    }
  }
}

Partial Pivoting Algorithm

David Mayerich

STIM Laboratory, University of Houston

void ppivot(double* A, double* b, int N){
  int* L;
  init(L);								// initialize the index vector
  for(int k=0; k<N-1; k++){ 			// for each pivot
  	int p = select_pivot(A, L, k, N);	// select the best pivot
    for(int i=k+1; i<N; i++){ 			// for each equation
      float m=A[i*N+k] / A[k*N+k];		// calculate scale
      for(int j=k; j<N; j++){ 			// for each coefficient
        A[i*N+j] -= m * A[k*N+j];		// subtract
      }
      b[i] -= m * b[k];					// forward elimination
    }
  }
}

Partial Pivoting Algorithm

David Mayerich

STIM Laboratory, University of Houston

void ppivot(double* A, double* b, int N){
  int* L;
  init(L);								// initialize the index vector
  for(int k=0; k<N-1; k++){ 			// for each pivot
  	int p = select_pivot(A, L, k, N);	// select the best pivot
    swap(L, p, k);						// swap indices
    for(int i=k+1; i<N; i++){ 			// for each equation
      float m=A[i*N+k] / A[k*N+k];		// calculate scale
      for(int j=k; j<N; j++){ 			// for each coefficient
        A[i*N+j] -= m * A[k*N+j];		// subtract
      }
      b[i] -= m * b[k];					// forward elimination
    }
  }
}

Partial Pivoting Algorithm

David Mayerich

STIM Laboratory, University of Houston

void ppivot(double* A, double* b, int N){
  int* L;
  init(L);								// initialize the index vector
  for(int k=0; k<N-1; k++){ 			// for each pivot
  	int p = select_pivot(A, L, k, N);	// select the best pivot
    swap(L, p, k);						// swap indices
    for(int i=k+1; i<N; i++){ 			// for each equation
      float m=A[L[i]*N+k] / A[L[k]*N+k];// calculate scale
      for(int j=k; j<N; j++){ 			// for each coefficient
        A[L[i]*N+j] -= m * A[L[k]*N+j];	// subtract
      }
      b[L[i]] -= m * b[L[k]];			// forward elimination
    }
  }
}

Pivoting Strategies

Scaled Partial Pivoting

Full Pivoting

David Mayerich

STIM Laboratory, University of Houston

Scaling Challenges

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{ccc|c} 8.00 & 1.00 & 1.00 & 10.0\\ 9.00 & 2000 & 16700 & 18700\\ 7.00 & 1.00 & 3.00 & 11.0 \end{array} \right]

  • Using partial pivoting, \(r_2\) will be selected as the first pivot, resulting in:

\left[ \begin{array}{ccc|c} 9.00 & 2000 & 16700 & 18700\\ 0.00 & -1780 & -14800 & -16600\\ 0.00 & 0.00 & 100 & 0.00 \end{array} \right]
100x_3 = 0.00
-1780 x_2 - 14800 (0.00) = -16600
9.00x_1+2000(9.33)+16700(0.00)=18700
\text{where}\quad x_1 = x_2 = x_3 = 1.00
x_3 = 0.00
\epsilon_r = 100\%
x_2 = 9.33
\epsilon_r = 833\%
x_1=0.00
\epsilon_r = 100\%

Scaled Partial Pivoting

  • Divide each equation by its largest coefficient:

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{ccc|c} 8.00 & 1.00 & 1.00 & 10.0\\ 9.00 & 2000 & 16700 & 18700\\ 7.00 & 1.00 & 3.00 & 11.0 \end{array} \right]

  • Applying partial pivoting to this scaled matrix yields:

\left[ \begin{array}{ccc|c} 1.00 & 0.125 & 0.125 & 1.25\\ 0.00 & 0.118 & 0.999 & 1.10\\ 0.00 & 0.018 & 0.304 & 0.32 \end{array} \right]
\left[ \begin{array}{ccc|c} 1.00 & 0.125 & 0.125 & 1.25\\ 0.00 & 0.118 & 0.999 & 1.10\\ 0.00 & 0.00 & 0.151 & 0.152 \end{array} \right]
x_1=1.03
x_2=0.763
x_3=1.03
\left[ \begin{array}{ccc|c} 1.00 & 0.125 & 0.125 & 1.25\\ 0.00536 & 0.119 & 1.00 & 1.11\\ 1.00 & 0.143 & 0.429 & 1.57 \end{array} \right]
\frac{r_1}{8.00}
\frac{r_2}{2000}
\frac{r_3}{7.00}
\epsilon_r=3\%
\epsilon_r=23\%
\epsilon_r=3\%

Is Scaling Necessary?

  • Use the scaled coefficient to select the pivot, but don't actually apply scaling:

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{ccc|c} 8.00 & 1.00 & 1.00 & 10.0\\ 9.00 & 2000 & 16700 & 18700\\ 7.00 & 1.00 & 3.00 & 11.0 \end{array} \right]
\begin{bmatrix} 1.00\\ 0.00536\\ 1.00 \end{bmatrix}
\left[ \begin{array}{ccc|c} 8.00 & 1.00 & 1.00 & 10.0\\ 0.00 & 2000 & 16700 & 18700\\ 0.00 & 0.125 & 2.13 & 2.25 \end{array} \right]
\begin{bmatrix} 0.120\\ 0.0587 \end{bmatrix}
\left[ \begin{array}{ccc|c} 8.00 & 1.00 & 1.00 & 10.0\\ 0.00 & 2000 & 16700 & 18700\\ 0.00 & 0.00 & 1.09 & 1.08 \end{array} \right]
x_1=0.989
x_2=1.10
x_3=0.991
\epsilon_r=1.1\%
\epsilon_r=10\%
\epsilon_r=0.9\%
x_1=1.03
x_2=0.763
x_3=1.03
\epsilon_r=3\%
\epsilon_r=23\%
\epsilon_r=3\%

w/ scaling

Full Pivoting

  • Swap rows and columns to optimize for the largest coefficient:

David Mayerich

STIM Laboratory, University of Houston

\left[ \begin{array}{ccc|c} 8.00 & 1.00 & 1.00 & 10.0\\ 9.00 & 2000 & 16700 & 18700\\ 7.00 & 1.00 & 3.00 & 11.0 \end{array} \right]
\left[ \begin{array}{ccc|c} 9.00 & 2000 & 16700 & 18700\\ 8.00 & 1.00 & 1.00 & 10.0\\ 7.00 & 1.00 & 3.00 & 11.0 \end{array} \right]
\left[ \begin{array}{ccc|c} 16700 & 2000 & 9.00 & 18700\\ 1.00 & 1.00 & 8.00 & 10.0\\ 3.00 & 1.00 & 7.00 & 11.0 \end{array} \right]
\left[ \begin{array}{ccc|c} 16700 & 2000 & 9.00 & 18700\\ 0.00 & 0.88 & 8.00 & 8.88\\ 0.00 & 0.64 & 7.00 & 7.63 \end{array} \right]
\left[ \begin{array}{ccc|c} 16700 & 2000 & 9.00 & 18700\\ 0.00 & 8.00 & 0.88 & 8.88\\ 0.00 & 0.00 & -0.13 & -0.14 \end{array} \right]
\begin{bmatrix} x_3\\ x_2\\ x_1 \end{bmatrix}
\begin{bmatrix} x_3\\ x_2\\ x_1 \end{bmatrix}
\begin{bmatrix} x_3\\ x_1\\ x_2 \end{bmatrix}
x_1=0.991
x_2=1.08
x_3=1.00
\epsilon_r=0.9\%
\epsilon_r=8\%
\epsilon_r=0\%

Pivoting Strategies

  • Partial pivoting:

David Mayerich

STIM Laboratory, University of Houston

  • Scaling with partial pivoting:

  • Scaled partial pivoting (use the scale to select the pivot):

  • Full pivoting:

x_3 = 0.00
\epsilon_r = 100\%
x_2 = 9.33
\epsilon_r = 833\%
x_1=0.00
\epsilon_r = 100\%
x_3 = 1.03
\epsilon_r = 3\%
x_2 = 0.763
\epsilon_r = 23\%
x_1=1.03
\epsilon_r = 3\%
x_3 = 0.989
\epsilon_r = 1.1\%
x_2 = 1.10
\epsilon_r = 10\%
x_1=0.991
\epsilon_r = 0.9\%
x_3 = 0.991
\epsilon_r = 0.9\%
x_2 = 1.08
\epsilon_r = 8\%
x_1=1.00
\epsilon_r = 0\%

LU Decomposition

Solving Multiple Right-Hand Sides

Matrix Factorization

Benefits for LU Decomposition

David Mayerich

STIM Laboratory, University of Houston

Row Reduction as Matrix Multiplication

  • Consider row reduction of the following matrix:

David Mayerich

STIM Laboratory, University of Houston

\mathbf{A} = \begin{bmatrix} 3 & 6 & 9\\ 2 & 5 & -2\\ 1 & 3 & -1 \end{bmatrix}
\begin{bmatrix} 3 & 6 & 9\\ 0 & 1 & -8\\ 0 & 1 & -4 \end{bmatrix}
m=2/3
m=1/3
\begin{bmatrix} 3 & 6 & 9\\ 0 & 1 & -8\\ 0 & 0 & 4 \end{bmatrix}
m=1

  • This can be encoded in the matrix multiplication:

\begin{bmatrix} 1 & 0 & 0\\ -\frac{2}{3} & 1 & 0\\ -\frac{1}{3} & -1 & 1 \end{bmatrix} \begin{bmatrix} 3 & 6 & 9\\ 2 & 5 & -2\\ 1 & 3 & -1 \end{bmatrix}
= \begin{bmatrix} 3 & 6 & 9\\ 0 & 1 & -8\\ 0 & 0 & 4 \end{bmatrix}

unitriangular

matrix

Unitriangular Matrix

  • A matrix \(\mathbf{A}\) is unitriangular if it is:

    • all diagonal elements are one

    • all lower-left or upper-right are zero

David Mayerich

STIM Laboratory, University of Houston

\begin{bmatrix} 1 & 0 & \cdots & 0\\ & 1 & \ddots & \vdots\\ & & \ddots & 0\\ & & & 1 \end{bmatrix}
= \begin{bmatrix} \mathbf{T}_1 & \mathbf{0}\\ \mathbf{M} & \mathbf{T}_2 \end{bmatrix}
  • \(\mathbf{0}\) is a matrix of all zeros
  • \(\mathbf{M}\) is a general matrix
  • \(\mathbf{T}\) is a smaller unitriangular matrix

  • The inverse of a unitriangular matrix is:

\begin{bmatrix} \mathbf{T}_1 & \mathbf{0}\\ \mathbf{M} & \mathbf{T}_2 \end{bmatrix}^{-1} = \begin{bmatrix} \mathbf{T}_1^{-1} & \mathbf{0}\\ -\mathbf{M} & \mathbf{T}_2^{-1} \end{bmatrix}

(negate all non-zero off-diagonal elements)

\begin{bmatrix} 1 & 0 & 0\\ -\frac{2}{3} & 1 & 0\\ -\frac{1}{3} & -1 & 1 \end{bmatrix}^{-1}
= \begin{bmatrix} 1 & 0 & 0\\ \frac{2}{3} & 1 & 0\\ \frac{1}{3} & 1 & 1 \end{bmatrix}

Matrix Factorization

  • We can perform row reduction using a unidiagonal matrix of scale factors:

David Mayerich

STIM Laboratory, University of Houston

\begin{bmatrix} 3 & 6 & 9\\ 2 & 5 & -2\\ 1 & 3 & -1 \end{bmatrix}
\begin{bmatrix} 1 & 0 & 0\\ -\frac{2}{3} & 1 & 0\\ -\frac{1}{3} & -1 & 1 \end{bmatrix}
= \begin{bmatrix} 3 & 6 & 9\\ 0 & 1 & -8\\ 0 & 0 & 4 \end{bmatrix}

  • And invert that process using the inverse matrix:

\mathbf{T}
\mathbf{A}
\mathbf{U}
=
\begin{bmatrix} 1 & 0 & 0\\ \frac{2}{3} & 1 & 0\\ \frac{1}{3} & 1 & 1 \end{bmatrix}
\begin{bmatrix} 3 & 6 & 9\\ 0 & 1 & -8\\ 0 & 0 & 4 \end{bmatrix}
= \begin{bmatrix} 3 & 6 & 9\\ 2 & 5 & -2\\ 1 & 3 & -1 \end{bmatrix}
\mathbf{L}
\mathbf{U}
\mathbf{A}
=
\mathbf{L}= \mathbf{T}^{-1}

LU Decomposition

LU Decomposition is the process of factorizing a matrix \(\mathbf{A}=\mathbf{L}\mathbf{U}\) into a lower unitriangular matrix \(\mathbf{L}\) and an upper triangular matrix \(\mathbf{U}\)

David Mayerich

STIM Laboratory, University of Houston

  • Given a matrix \(\mathbf{A}\), calculate \(\mathbf{L}\) and \(\mathbf{U}\):

    1. Initialize \(\mathbf{L}_0 = \mathbf{I}\) with the identity matrix

    2. Perform Gaussian elimination of \(\mathbf{A}\), which produces \(\mathbf{U}\)

    3. For each iteration \(n\) of Gaussian elimination, store the scale factors in \(\mathbf{L}_n\)

LU Decomposition Algorithm

David Mayerich

STIM Laboratory, University of Houston

lu_decomposition(double* A, double* L, int N){
  identity(L, N);					// initialize L
  for(int k=0; k<N-1; k++){ 		// for each pivot
    for(int i=k+1; i<N; i++){ 		// for each equation
      float m=A[i*N+k] / A[k*N+k]; 	// calculate scale
      L[i*N+k] = m;					// store the scale factor
      for(int j=k; j<N; j++){ 		// for each coefficient
        A[i*N+j] -= m * A[k*N+j]; 	// subtract
      }
	}
  }
}

Solving Multiple Right-Hand Sides

  • Solve the following two linear systems:

David Mayerich

STIM Laboratory, University of Houston

\begin{bmatrix} 8.00 & 1.00 & 1.00\\ 9.00 & 2000 & 16700\\ 7.00 & 1.00 & 3.00 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix} 10.0\\ 18700\\ 11.0 \end{bmatrix}
\begin{bmatrix} 8.00 & 1.00 & 1.00\\ 9.00 & 2000 & 16700\\ 7.00 & 1.00 & 3.00 \end{bmatrix} \begin{bmatrix} x_1\\ x_2\\ x_3 \end{bmatrix} = \begin{bmatrix} 13.0\\ 54100\\ 18.0 \end{bmatrix}

  • Using a direct method would require \(O(n)\) time for each system

  • What is the complexity for solving a linear system \(\mathbf{A}\mathbf{x}=\mathbf{b}\) if \(\mathbf{A}\) is triangular?

  • Factorize \(\mathbf{A}\) and solve:

O(n^2)
\mathbf{A} \mathbf{x}_1=\mathbf{b}_1
\mathbf{A} \mathbf{x}_2=\mathbf{b}_2
\mathbf{LU} \mathbf{x}_1=\mathbf{b}_1
\mathbf{LU} \mathbf{x}_2=\mathbf{b}_2

LU decomposition:

\mathbf{L} \mathbf{z}_1=\mathbf{b}_1
\mathbf{L} \mathbf{z}_2=\mathbf{b}_2

Solve for \(\mathbf{z}=\mathbf{Ux}\):

Solve for \(\mathbf{x}\):

\mathbf{U} \mathbf{x}_1=\mathbf{z}_1
\mathbf{U} \mathbf{x}_2=\mathbf{z}_2

Solving Multiple Right-Hand Sides

David Mayerich

STIM Laboratory, University of Houston

\mathbf{A} = \begin{bmatrix} 8.00 & 1.00 & 1.00\\ 9.00 & 2000 & 16700\\ 7.00 & 1.00 & 3.00 \end{bmatrix}
\mathbf{LU} = \begin{bmatrix} 8.00 & 1.00 & 1.00\\ 9.00 & 2000 & 16700\\ 7.00 & 1.00 & 3.00 \end{bmatrix}