Lecture 8: Gaussian EliminationAMath 352Wed., Apr. 141 / 13Gaussian Elimination as LU FactorizationSolve A~x =~b. A is n by n nonsingular. Example:1 2 34 5 67 8 0x1x2x3=102.1 2 3 | 14 5 6 | 07 8 0 | 2→1 2 3 | 10 −3 −6 | −40 −6 −21 | −5;that is, ifL1=1 0 0−4 1 0−7 0 1,thenL1A =1 0 0−4 1 0−7 0 11 2 34 5 67 8 0=1 2 30 −3 −60 −6 −21.2 / 13Gaussian Elimination as LU FactorizationSolve A~x =~b. A is n by n nonsingular. Example:1 2 34 5 67 8 0x1x2x3=102.1 2 3 | 14 5 6 | 07 8 0 | 2→1 2 3 | 10 −3 −6 | −40 −6 −21 | −5;that is, ifL1=1 0 0−4 1 0−7 0 1,thenL1A =1 0 0−4 1 0−7 0 11 2 34 5 67 8 0=1 2 30 −3 −60 −6 −21.2 / 13Gaussian Elimination as LU Factorization, Cont.L1=1 0 0−4 1 0−7 0 1, L1A =1 2 30 −3 −60 −6 −21.1 2 3 | 10 −3 −6 | −40 −6 −21 | −5→1 2 3 | 10 −3 −6 | −40 0 −9 | 3;that is, ifL2=1 0 00 1 00 −2 1,thenL2L1A =1 0 00 1 00 −2 1,1 2 30 −3 −60 −6 −21=1 2 30 −3 −60 0 −9.3 / 13Gaussian Elimination as LU Factorization, Cont.L1=1 0 0−4 1 0−7 0 1, L1A =1 2 30 −3 −60 −6 −21.1 2 3 | 10 −3 −6 | −40 −6 −21 | −5→1 2 3 | 10 −3 −6 | −40 0 −9 | 3;that is, ifL2=1 0 00 1 00 −2 1,thenL2L1A =1 0 00 1 00 −2 1,1 2 30 −3 −60 −6 −21=1 2 30 −3 −60 0 −9.3 / 13Gaussian Elimination as LU Factorization, Cont.L2L1A = U ⇒ L1A = L−12U ⇒ A = L−11L−12U.Claim: We have factored A in the form LU, where L is lowertriangular with 1’s on its diagonal and U is upper triangular.To get L−11and L−12from L1and L2, just change the signs of theoff-diagonal entries:L1=1 0 0−4 1 0−7 0 1⇒ L−11=1 0 04 1 07 0 1,L2=1 0 00 1 00 −2 1⇒ L−12=1 0 00 1 00 2 1.To compute the product L−11L−12, just put 1’s on the diagonal andtake the off-diagonal entries of both L−11and L−12:L := L−11L−12=1 0 04 1 07 2 1.4 / 13Gaussian Elimination as LU Factorization, Cont.L2L1A = U ⇒ L1A = L−12U ⇒ A = L−11L−12U.Claim: We have factored A in the form LU, where L is lowertriangular with 1’s on its diagonal and U is upper triangular.To get L−11and L−12from L1and L2, just change the signs of theoff-diagonal entries:L1=1 0 0−4 1 0−7 0 1⇒ L−11=1 0 04 1 07 0 1,L2=1 0 00 1 00 −2 1⇒ L−12=1 0 00 1 00 2 1.To compute the product L−11L−12, just put 1’s on the diagonal andtake the off-diagonal entries of both L−11and L−12:L := L−11L−12=1 0 04 1 07 2 1.4 / 13Gaussian Elimination as LU Factorization, Cont.Check:LU =1 0 04 1 07 2 11 2 30 −3 −60 0 −9=1 2 34 5 67 8 0.Use the LU factors to solve linear systems with different right-handsides:L~y =~b, U~x =~y ⇒ L(U~x) = (LU)~x = A~x = L~y =~b.5 / 13Gaussian Elimination as LU Factorization, Cont.Check:LU =1 0 04 1 07 2 11 2 30 −3 −60 0 −9=1 2 34 5 67 8 0.Use the LU factors to solve linear systems with different right-handsides:L~y =~b, U~x =~y ⇒ L(U~x) = (LU)~x = A~x = L~y =~b.5 / 13Generalize to n by n Matricesa11a12. . . a1n| b1a21a22. . . a2n| b2......... |...an1an2. . . ann| bn−→a11a12. . . a1n| b10 ˜a22. . . ˜a2n|˜b2......... |...0 ˜an2. . . ˜ann|˜bn.% First step of Gaussian Elimination.L = eye(n); % Initialize L to the identity.for i=2:n,mult = A(i,1)/A(1,1);L(i,1) = mult;A(i,:) = A(i,:) - mult*A(1,:);b(i) = b(i) - mult*b(1);end;6 / 13Generalize to n by n Matricesa11a12. . . a1n| b1a21a22. . . a2n| b2......... |...an1an2. . . ann| bn−→a11a12. . . a1n| b10 ˜a22. . . ˜a2n|˜b2......... |...0 ˜an2. . . ˜ann|˜bn.% First step of Gaussian Elimination.L = eye(n); % Initialize L to the identity.for i=2:n,mult = A(i,1)/A(1,1);L(i,1) = mult;A(i,:) = A(i,:) - mult*A(1,:);b(i) = b(i) - mult*b(1);end;6 / 13Gaussian Elimination Without Pivotinga11a12. . . a1n| b10 ˜a22. . . ˜a2n|˜b2......... |...0 ˜an2. . . ˜ann|˜bn−→a11a12. . . a1n| b10 ˜a22. . . ˜a2n|˜b2...... |...0 0 . . .˜˜ann|˜˜bn.% Gaussian elimination without pivoting.for j=1:n-1, % Loop over columns.for i=j+1:n, % Loop over rows below row j.mult = A(i,j)/A(j,j);L(i,j) = mult;A(i,:) = A(i,:) - mult*A(j,:);b(i) = b(i) - mult*b(j);end;end;7 / 13Gaussian Elimination Without Pivotinga11a12. . . a1n| b10 ˜a22. . . ˜a2n|˜b2......... |...0 ˜an2. . . ˜ann|˜bn−→a11a12. . . a1n| b10 ˜a22. . . ˜a2n|˜b2...... |...0 0 . . .˜˜ann|˜˜bn.% Gaussian elimination without pivoting.for j=1:n-1, % Loop over columns.for i=j+1:n, % Loop over rows below row j.mult = A(i,j)/A(j,j);L(i,j) = mult;A(i,:) = A(i,:) - mult*A(j,:);b(i) = b(i) - mult*b(j);end;end;7 / 13Gaussian Elimination Without Pivoting, Cont.Some improvements:% Gaussian elimination without pivoting.for j=1:n-1, % Loop over columns.for i=j+1:n, % Loop over rows below row j.mult = A(i,j)/A(j,j);% L(i,j) = mult;% Save storage by storing the below-diagonal entries of L in A.A(i,j) = mult;% A(i,:) = A(i,:) - mult*A(j,:);% Columns to the left of column j just have 0’s. 0 - mult*0 is still 0!A(i,j:n) = A(i,j:n) - mult*A(j,j:n);b(i) = b(i) - mult*b(j);end;end;8 / 13Speed of a CodeMany things affect the speed of a code:IHow many floating point operations (⊕, , ⊗, ) areperformed.IHow many memory references.To do an operation, operands are fetched from memory(slow), operated on in registers, and stored back to mem-ory (slow). When you bring numbers into the registers, doas much as you can with them before storing the resultback to memory.IHow well the code vectorizes.IHow well the code parallelizes.We will talk mostly about the first, but say a little about theothers, as they may actually be more important.9 / 13Speed of a CodeMany things affect the speed of a code:IHow many floating point operations