Lecture 9: Op Count for Gaussian Elimination;Pivoting for StabilityAMath 352Fri., Apr. 161 / 17Operation Count for GE% 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); % 1 div.A(i,j:n) = A(i,j:n) - mult*A(j,j:n); % (n-j+1) mults. and (n-j+1) adds.b(i) = b(i) - mult*b(j); % 1 mult. and 1 add.end;end;No. of operations =Pn−1j=1Pni=j+1(no. of ops done inside thedouble loop) .Ops done inside the double loop: 2(n − j) + 5Total ops:n−1Xj=1nXi=j+1[2(n − j) + 5] =n−1Xj=1[2(n − j)2+ 5(n − j)].2 / 17Operation Count for GE% 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); % 1 div.A(i,j:n) = A(i,j:n) - mult*A(j,j:n); % (n-j+1) mults. and (n-j+1) adds.b(i) = b(i) - mult*b(j); % 1 mult. and 1 add.end;end;No. of operations =Pn−1j=1Pni=j+1(no. of ops done inside thedouble loop) .Ops done inside the double loop: 2(n − j) + 5Total ops:n−1Xj=1nXi=j+1[2(n − j) + 5] =n−1Xj=1[2(n − j)2+ 5(n − j)].2 / 17Operation Count for GE% 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); % 1 div.A(i,j:n) = A(i,j:n) - mult*A(j,j:n); % (n-j+1) mults. and (n-j+1) adds.b(i) = b(i) - mult*b(j); % 1 mult. and 1 add.end;end;No. of operations =Pn−1j=1Pni=j+1(no. of ops done inside thedouble loop) .Ops done inside the double loop: 2(n − j) + 5Total ops:n−1Xj=1nXi=j+1[2(n − j) + 5] =n−1Xj=1[2(n − j)2+ 5(n − j)].2 / 17Operation Count for GE% 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); % 1 div.A(i,j:n) = A(i,j:n) - mult*A(j,j:n); % (n-j+1) mults. and (n-j+1) adds.b(i) = b(i) - mult*b(j); % 1 mult. and 1 add.end;end;No. of operations =Pn−1j=1Pni=j+1(no. of ops done inside thedouble loop) .Ops done inside the double loop: 2(n − j) + 5Total ops:n−1Xj=1nXi=j+1[2(n − j) + 5] =n−1Xj=1[2(n − j)2+ 5(n − j)].2 / 17Operation Count for GE, Cont.n−1Xj=1[2(n − j)2+ 5(n − j)].mXi=1i = 1+2+. . .+m =m(m + 1)2, (no. of terms times avg term).mXi=1i2=m(m + 1)(2m + 1)6Thereforen−1Xj=1[2(n−j)2+5(n−j)] = 2(n − 1)n(2n − 1)6+5(n − 1)n2=23n3+O(n2).3 / 17Operation Count for GE, Cont.n−1Xj=1[2(n − j)2+ 5(n − j)].mXi=1i = 1+2+. . .+m =m(m + 1)2, (no. of terms times avg term).mXi=1i2=m(m + 1)(2m + 1)6Thereforen−1Xj=1[2(n−j)2+5(n−j)] = 2(n − 1)n(2n − 1)6+5(n − 1)n2=23n3+O(n2).3 / 17If you forget the exact formula:mXi=1ip≈Zm0xpdx =xp+1p + 1m0=mp+1p + 1.More precisely,mXi=1ip=mp+1p + 1+ O(mp).4 / 17Op Count for Solving Upper Triangular Systemu11u12. . . u1n0 u22. . . u2n.........0 0 . . . unnx1x2...xn=y1y2...yn.xn= yn/unn,xn−1= (yn−1− un−1,nxn)/un−1,n−1.In genl.,xi= (yi−nXj=i+1ui,jxj)/uii, i = n, n − 1, . . . , 1.Op count:Pni=1(work to compute xi):nXi=1(1 + 2(n − i)) = n + 2(n − 1)n2= n2+ O(n).5 / 17Solving Linear Systems with Multiple Right-Hand SidesIWork for GE =23n3+ O(n2).If it takes 1 second to solve a 1000 by 1000 linear system,it takes 23= 8 seconds to solve a 2000 by 2000 linearsystem.IWork for solving a triangular system (upper or lower) =n2+ O(n).If you want to solve A~x1=~b1, A~x2=~b2, . . ., factorA = LU once, then solve: L~yj=~bj, U~xj=~yj, j = 1, 2, . . ..Work is about23n3+ 2n2(no. of right-hand sides).6 / 17PivotingThe GE codes that we have looked at so far may fail! What if thedenominator in mult = A(i,j)/A(j,j) is 0?One possible fix: If A(j,j) = 0, search through that column to finda nonzero, say, A(k,j), k > j. Interchange row k and row j and dothe elimination. If no nonzero is found in the column, return withan error message: ”Matrix is singular.”Why is this matrix singular?a11a12a13a140 a22a23a240 0 0 a340 0 0 a44.7 / 17PivotingThe GE codes that we have looked at so far may fail! What if thedenominator in mult = A(i,j)/A(j,j) is 0?One possible fix: If A(j,j) = 0, search through that column to finda nonzero, say, A(k,j), k > j. Interchange row k and row j and dothe elimination. If no nonzero is found in the column, return withan error message: ”Matrix is singular.”Why is this matrix singular?a11a12a13a140 a22a23a240 0 0 a340 0 0 a44.7 / 17PivotingThe GE codes that we have looked at so far may fail! What if thedenominator in mult = A(i,j)/A(j,j) is 0?One possible fix: If A(j,j) = 0, search through that column to finda nonzero, say, A(k,j), k > j. Interchange row k and row j and dothe elimination. If no nonzero is found in the column, return withan error message: ”Matrix is singular.”Why is this matrix singular?a11a12a13a140 a22a23a240 0 0 a340 0 0 a44.7 / 17Pivoting, Cont.In exact arithmetic, this is fine. But not a good idea in finiteprecision arithmetic. Suppose a pivot would be 0 in exact arith.,but we computed . We would use this tiny (and incorrect) entryto eliminate by subtracting huge multiples of this row from others.Example.10−2011 1x1x2=12.Soln. is very nearly x1= x2= 1; precisely,x1= 1 + 1/(1020− 1), x2= 1 −1/(1020− 1).8 / 17Pivoting, Cont.In exact arithmetic, this is fine. But not a good idea in finiteprecision arithmetic. Suppose a pivot would be 0 in exact arith.,but we computed . We would use this tiny (and incorrect) entryto eliminate by subtracting huge multiples of this row from others.Example.10−2011 1x1x2=12.Soln. is very nearly x1= x2= 1; precisely,x1= 1 + 1/(1020− 1), x2= 1 −1/(1020− 1).8 / 17Pivoting, Cont.Use GE, rounding to 16 decimal places:10−201 | 11 1 | 2−→10−201 | 10 round(1 −1020) | round(2 −1020)=10−201 | 10 −1020| −1020.x2= 1, 10−20x1+ 1 = 1 ⇒ x1= 0. WRONG!9 / 17Pivoting, Cont.To avoid this, interchange rows:1 1 | 210−201 | 1−→1 1 | 20 round(1 −10−20) | round(1 −2 ·10−20)=1 1 | 20 1 | 1=x2= 1, x1+ 1 = 2 ⇒ x1= 1. RIGHT.10 / 17Partial PivotingUse the largest entry in abs. value on or below the diagonal in thecolumn as the pivot element.1 2 34 5 67 8 0P1−→7 8 04 5 61 2 3L1−→7 8 00 3/7 60 6/7 3P2−→7 8 00 6/7 30 3/7 6L2−→7 8 00 6/7 30 0 9/2.P1=0 0 10 1 01 0