chap02_8up - Existence Uniqueness and Conditioning Solving...

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Outline Scientific Computing: An Introductory Survey 1 Existence, Uniqueness, and Conditioning 2 Solving Linear Systems Department of Computer Science University of Illinois at Urbana-Champaign 3 Special Types of Linear Systems Copyright c 2002. Reproduction permitted for noncommercial, educational use only. 4 Software for Linear Systems Chapter 2 – Systems of Linear Equations Prof. Michael T. Heath Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Michael T. Heath 1 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Systems of Linear Equations Scientific Computing 2 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Singularity and Nonsingularity Given m × n matrix A and m-vector b, find unknown n-vector x satisfying Ax = b n × n matrix A is nonsingular if it has any of following equivalent properties System of equations asks “Can b be expressed as linear combination of columns of A?” 1 2 rank(A) = n 4 Solution may or may not exist, and may or may not be unique det(A) = 0 3 If so, coefficients of linear combination are given by components of solution vector x Inverse of A, denoted by A−1 , exists For any vector z = 0, Az = 0 For now, we consider only square case, m = n Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems 3 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Existence and Uniqueness Scientific Computing 4 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Geometric Interpretation In two dimensions, each equation determines straight line in plane Existence and uniqueness of solution to Ax = b depend on whether A is singular or nonsingular Solution is intersection point of two lines Can also depend on b, but only in singular case If two straight lines are not parallel (nonsingular), then intersection point is unique If b ∈ span(A), system is consistent A nonsingular singular singular b arbitrary b ∈ span(A) b ∈ span(A) / Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems If two straight lines are parallel (singular), then lines either do not intersect (no solution) or else coincide (any point along line is solution) # solutions one (unique) infinitely many none In higher dimensions, each equation determines hyperplane; if matrix is nonsingular, intersection of hyperplanes is unique solution 5 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Example: Nonsingularity Scientific Computing 6 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Example: Singularity 2 × 2 system 2 × 2 system 2x1 + 3x2 = b1 5x1 + 4x2 = b2 Ax = or in matrix-vector notation Ax = x1 b = 1 =b x2 b2 With b = 4 7 T , then x = 1 2 Scientific Computing T , there is no solution T is nonsingular regardless of value of b Michael T. Heath x1 b = 1 =b x2 b2 is singular regardless of value of b 23 54 For example, if b = 8 13 solution 23 46 T T With b = 4 8 , x = γ (4 − 2γ )/3 is solution for any real number γ , so there are infinitely many solutions is unique 7 / 88 Michael T. Heath Scientific Computing 8 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Singularity and Nonsingularity Norms Condition Number Error Bounds Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Vector Norms Singularity and Nonsingularity Norms Condition Number Error Bounds Example: Vector Norms Drawing shows unit sphere in two dimensions for each norm Magnitude, modulus, or absolute value for scalars generalizes to norm for vectors We will use only p-norms, defined by 1/p n x |xi |p = p i=1 for integer p > 0 and n-vector x Important special cases n i=1 |xi | n 2 1/2 i=1 |xi | 1-norm: x 1 = 2-norm: x 2 = ∞-norm: x ∞ Norms have following values for vector shown = maxi |xi | Michael T. Heath x 9 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds 2 Michael T. Heath = 2.0 x ∞ = 1.6 Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Equivalence of Norms 10 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Properties of Vector Norms In general, for any vector x in Rn , x However, we also have √ x 1 ≤ n x 2, x x < interactive example > Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems = 2.8 1 2 ≤ √ nx 1 ≥x ∞, x For any vector norm 2 ≥x 1 ≤n x ∞ x > 0 if x = 0 γ x = |γ | · x for any scalar γ x+y ≤ x + y In more general treatment, these properties taken as definition of vector norm Thus, for given n, norms differ by at most a constant, and hence are equivalent: if one is small, they must all be proportionally small. Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems (triangle inequality) ∞ Scientific Computing Useful variation on triangle inequality | x − y |≤ x−y 11 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Matrix Norms 12 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Matrix Norms Matrix norm corresponding to vector 1-norm is maximum absolute column sum Matrix norm corresponding to given vector norm is defined by Ax A = maxx=0 x n A 1 |aij | = max j i=1 Matrix norm corresponding to vector ∞-norm is maximum absolute row sum Norm of matrix measures maximum stretching matrix does to any vector in given vector norm n A ∞ |aij | = max i j =1 Handy way to remember these is that matrix norms agree with corresponding vector norms for n × 1 matrix Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Scientific Computing 13 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Properties of Matrix Norms Scientific Computing 14 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Condition Number Condition number of square nonsingular matrix A is defined by cond(A) = A · A−1 Any matrix norm satisfies A > 0 if A = 0 By convention, cond(A) = ∞ if A is singular γ A = |γ | · A for any scalar γ Since A+B ≤ A + B Matrix norms we have defined also satisfy A · A−1 = AB ≤ A · B Ax ≤ A · x for any vector x max x=0 Ax x · min x=0 Ax x −1 condition number measures ratio of maximum stretching to maximum shrinking matrix does to any nonzero vectors Large cond(A) means A is nearly singular Michael T. Heath Scientific Computing 15 / 88 Michael T. Heath Scientific Computing 16 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Singularity and Nonsingularity Norms Condition Number Error Bounds Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Properties of Condition Number Singularity and Nonsingularity Norms Condition Number Error Bounds Computing Condition Number Definition of condition number involves matrix inverse, so it is nontrivial to compute For any matrix A, cond(A) ≥ 1 Computing condition number from definition would require much more work than computing solution whose accuracy is to be assessed For identity matrix, cond(I ) = 1 For any matrix A and scalar γ , cond(γ A) = cond(A) For any diagonal matrix D = diag(di ), cond(D ) = max |di | min |di | In practice, condition number is estimated inexpensively as byproduct of solution process Matrix norm A is easily computed as maximum absolute column sum (or row sum, depending on norm used) < interactive example > Estimating A−1 at low cost is more challenging Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Scientific Computing 17 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Computing Condition Number, continued Scientific Computing 18 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Error Bounds Condition number yields error bound for computed solution to linear system From properties of norms, if Az = y , then ˆ Let x be solution to Ax = b, and let x be solution to ˆ Ax = b + ∆b z ≤ A−1 y ˆ If ∆x = x − x, then and bound is achieved for optimally chosen y ˆ b + ∆b = A(x) = A(x + ∆x) = Ax + A∆x Efficient condition estimators heuristically pick y with large ratio z / y , yielding good estimate for A−1 which leads to bound ∆x ∆b ≤ cond(A) x b Good software packages for linear systems provide efficient and reliable condition estimator for possible relative change in solution x due to relative change in right-hand side b < interactive example > Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Scientific Computing 19 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Error Bounds, continued Scientific Computing 20 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Error Bounds – Illustration Similar result holds for relative change in matrix: if ˆ (A + E )x = b, then In two dimensions, uncertainty in intersection point of two lines depends on whether lines are nearly parallel ∆x E ≤ cond(A) ˆ x A If input data are accurate to machine precision, then bound for relative error in solution x becomes ˆ x−x ≤ cond(A) x mach Computed solution loses about log10 (cond(A)) decimal digits of accuracy relative to accuracy of input Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Scientific Computing < interactive example > 21 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Error Bounds – Caveats Scientific Computing 22 / 88 Singularity and Nonsingularity Norms Condition Number Error Bounds Residual ˆ Residual vector of approximate solution x to linear system Ax = b is defined by Normwise analysis bounds relative error in largest components of solution; relative error in smaller components can be much larger ˆ r = b − Ax Componentwise error bounds can be obtained, but somewhat more complicated ˆ In theory, if A is nonsingular, then x − x = 0 if, and only if, r = 0, but they are not necessarily small simultaneously Conditioning of system is affected by relative scaling of rows or columns Since Ill-conditioning can result from poor scaling as well as near singularity Rescaling can help the former, but not the latter Michael T. Heath Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Scientific Computing ∆x r ≤ cond(A) ˆ ˆ x A·x small relative residual implies small relative error in approximate solution only if A is well-conditioned 23 / 88 Michael T. Heath Scientific Computing 24 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Singularity and Nonsingularity Norms Condition Number Error Bounds Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Residual, continued Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Solving Linear Systems ˆ If computed solution x exactly satisfies To solve linear system, transform it into one whose solution is same but easier to compute ˆ (A + E )x = b then r ˆ x A ≤ What type of transformation of linear system leaves solution unchanged? E A so large relative residual implies large backward error in matrix, and algorithm used to compute solution is unstable We can premultiply (from left) both sides of linear system Ax = b by any nonsingular matrix M without affecting solution Stable algorithm yields small relative residual regardless of conditioning of nonsingular system Solution to M Ax = M b is given by x = (M A)−1 M b = A−1 M −1 M b = A−1 b Small residual is easy to obtain, but does not necessarily imply computed solution is accurate Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 25 / 88 Example: Permutations Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 26 / 88 Example: Diagonal Scaling Permutation matrix P has one 1 in each row and column and zeros elsewhere, i.e., identity matrix with rows or columns permuted Row scaling: premultiplying both sides of system by nonsingular diagonal matrix D , DAx = Db, multiplies each row of matrix and right-hand side by corresponding diagonal entry of D , but solution x is unchanged Note that P −1 = P T Premultiplying both sides of system by permutation matrix, P Ax = P b, reorders rows, but solution x is unchanged Column scaling: postmultiplying A by D , ADx = b, multiplies each column of matrix by corresponding diagonal entry of D , which rescales original solution Postmultiplying A by permutation matrix, AP x = b, reorders columns, which permutes components of original solution x = (AD )−1 b = D −1 A−1 b x = (AP )−1 b = P −1 A−1 b = P T (A−1 b) Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 27 / 88 Triangular Linear Systems Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 28 / 88 Triangular Matrices What type of linear system is easy to solve? Two specific triangular forms are of particular interest If one equation in system involves only one component of solution (i.e., only one entry in that row of matrix is nonzero), then that component can be computed by division lower triangular : all entries above main diagonal are zero, aij = 0 for i < j upper triangular : all entries below main diagonal are zero, aij = 0 for i > j If another equation in system involves only one additional solution component, then by substituting one known component into it, we can solve for other component Successive substitution process described earlier is especially easy to formulate for lower or upper triangular systems If this pattern continues, with only one new solution component per equation, then all components of solution can be computed in succession. Any triangular matrix can be permuted into upper or lower triangular form by suitable row permutation System with this property is called triangular Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 29 / 88 Forward-Substitution Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy xi = bi − ij xj / ii , n i = 2, . . . , n xn = bn /unn , xi = bi − j =1 for j = 1 to n if j j = 0 then stop xj = bj / jj for i = j + 1 to n bi = bi − ij xj end end Michael T. Heath 30 / 88 Back-substitution for upper triangular system U x = b i− 1 11 , Scientific Computing Back-Substitution Forward-substitution for lower triangular system Lx = b x1 = b1 / Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems uij xj / uii , i = n − 1, . . . , 1 j =i+1 { loop over columns } { stop if matrix is singular } { compute solution component } { update right-hand side } Scientific Computing 31 / 88 for j = n to 1 if ujj = 0 then stop xj = bj /ujj for i = 1 to j − 1 bi = bi − uij xj end end Michael T. Heath { loop backwards over columns } { stop if matrix is singular } { compute solution component } { update right-hand side } Scientific Computing 32 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Example: Triangular Linear System Elimination To transform general linear system into triangular form, we need to replace selected nonzero entries of matrix by zeros 2 4 −2 x1 2 0 1 1 x2 = 4 00 4 x3 8 This can be accomplished by taking linear combinations of rows Using back-substitution for this upper triangular system, last equation, 4x3 = 8, is solved directly to obtain x3 = 2 Consider 2-vector a = Next, x3 is substituted into second equation to obtain x2 = 2 Michael T. Heath Scientific Computing a1 a2 If a1 = 0, then Finally, both x3 and x2 are substituted into first equation to obtain x1 = −1 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 1 0 −a2 /a1 1 Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 33 / 88 a1 a =1 a2 0 Elementary Elimination Matrices Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 34 / 88 Elementary Elimination Matrices, continued More generally, we can annihilate all entries below k th position in n-vector a by transformation 1 ··· 0 0 ··· 0 a1 a1 . .. . . . . .. . . . . . . . . . . . . . . . 0 · · · 1 0 · · · 0 ak ak Mk a = 0 · · · −mk+1 1 · · · 0 ak+1 = 0 . .. . . .. . . . . . . . . . . . . . . . . . 0 ··· −mn 0 · · · 1 an 0 Matrix Mk , called elementary elimination matrix, adds multiple of row k to each subsequent row, with multipliers mi chosen so that result is zero Mk is unit lower triangular and nonsingular Mk = I − mk eT , where mk = [0, . . . , 0, mk+1 , . . . , mn ]T k and ek is k th column of identity matrix − − Mk 1 = I + mk eT , which means Mk 1 = Lk is same as k Mk except signs of multipliers are reversed where mi = ai /ak , i = k + 1, . . . , n Divisor ak , called pivot, must be nonzero Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 35 / 88 Elementary Elimination Matrices, continued Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 36 / 88 Example: Elementary Elimination Matrices 2 For a = 4, −2 If Mj , j > k , is another elementary elimination matrix, with vector of multipliers mj , then Mk Mj 100 2 2 M1 a = −2 1 0 4 = 0 1 0 1 −2 0 = I − mk eT − mj eT + mk eT mj eT j j k k = I − mk eT − mj eT j k which means product is essentially “union,” and similarly for product of inverses, Lk Lj Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems and Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 37 / 88 L1 = Michael T. Heath Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 38 / 88 To reduce general linear system Ax = b to upper triangular form, first choose M1 , with a11 as pivot, to annihilate first column of A below first row L2 = − M2 1 1 0 0 0 1 0 = 0 −1/2 1 System becomes M1 Ax = M1 b, but solution is unchanged Next choose M2 , using a22 as pivot, to annihilate second column of M1 A below second row and 100 M1 M2 = −2 1 0 , 1 1 /2 1 Scientific Computing Gaussian Elimination Note that 100 2 1 0 , = −1 0 1 Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Example, continued − M1 1 100 2 2 M2 a = 0 1 0 4 = 4 0 1 /2 1 −2 0 System becomes M2 M1 Ax = M2 M1 b, but solution is still unchanged 1 0 0 1 0 L1 L 2 = 2 −1 −1/2 1 Scientific Computing Process continues for each successive column until all subdiagonal entries have been zeroed 39 / 88 Michael T. Heath Scientific Computing 40 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Gaussian Elimination, continued Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy LU Factorization Product Lk Lj is unit lower triangular if k < j , so − −1 L = M −1 = M1 1 · · · Mn−1 = L1 · · · Ln−1 Resulting upper triangular linear system Mn−1 · · · M1 Ax = Mn−1 · · · M1 b is unit lower triangular M Ax = M b By design, U = M A is upper triangular can be solved by back-substitution to obtain solution to original linear system Ax = b So we have Process just described is called Gaussian elimination with L unit lower triangular and U upper triangular A = LU Thus, Gaussian elimination produces LU factorization of matrix into triangular factors Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 41 / 88 LU Factorization, continued Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 42 / 88 Example: Gaussian Elimination Use Gaussian elimination to solve linear system 2 4 −2 x1 2 9 −3 x2 = 8 = b Ax = 4 −2 −3 7 x3 10 Having obtained LU factorization, Ax = b becomes LU x = b, and can be solved by forward-substitution in lower triangular system Ly = b, followed by back-substitution in upper triangular system U x = y To annihilate subdiagonal entries of first column of A, 100 2 4 −2 2 4 −2 9 −3 = 0 1 1 , M1 A = −2 1 0 4 101 −2 −3 7 01 5 100 2 2 M1 b = −2 1 0 8 = 4 101 10 12 Note that y = M b is same as transformed right-hand side in Gaussian elimination Gaussian elimination and LU factorization are two ways of expressing same solution process Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 43 / 88 Example, continued Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 44 / 88 Example, continued We have reduced original system to equivalent upper triangular system 2 4 −2 x1 2 0 1 1 x2 = 4 = M b Ux = 00 4 x3 8 To annihilate subdiagonal entry of second column of M1 A, 1 00 2 4 −2 2 4 −2 1 0 0 1 1 = 0 1 1 = U , M2 M1 A = 0 0 −1 1 01 5 00 4 which can now be solved by back-substitution to obtain −1 x = 2 2 1 00 2 2 1 0 4 = 4 = M b M2 M1 b = 0 0 −1 1 12 8 Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 45 / 88 Example, continued Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 46 / 88 Row Interchanges Gaussian elimination breaks down if leading diagonal entry of remaining unreduced matrix is zero at any stage To write out LU factorization explicitly, 100 100 100 2 1 0 0 1 0 = 2 1 0 = L L1 L2 = −1 0 1 0 1 1 −1 1 1 Easy fix: if diagonal entry in column k is zero, then interchange row k with some subsequent row having nonzero entry in column k and then proceed as usual If there is no nonzero on or below diagonal in column k , then there is nothing to do at this stage, so skip to next column so that 2 4 −2 100 2 4 −2 4 = 2 1 0 0 1 9 −3 1 = LU A= −2 −3 7 −1 1 1 00 4 Zero on diagonal causes resulting upper triangular matrix U to be singular, but LU factorization can still be completed Subsequent back-substitution will fail, however, as it should for singular matrix Michael T. Heath Scientific Computing 47 / 88 Michael T. Heath Scientific Computing 48 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Partial Pivoting Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy LU Factorization with Partial Pivoting With partial pivoting, each Mk is preceded by permutation Pk to interchange rows to bring entry of largest magnitude into diagonal pivot position Still obtain M A = U , with U upper triangular, but now In principle, any nonzero value will do as pivot, but in practice pivot should be chosen to minimize error propagation To avoid amplifying previous rounding errors when multiplying remaining portion of matrix by elementary elimination matrix, multipliers should not exceed 1 in magnitude M = Mn−1 Pn−1 · · · M1 P1 M −1 L= is still triangular in general sense, but not necessarily lower triangular Alternatively, we can write This can be accomplished by choosing entry of largest magnitude on or below diagonal as pivot at each stage PA = LU Such partial pivoting is essential in practice for numerically stable implementation of Gaussian elimination for general linear systems < interactive example > Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems where P = Pn−1 · · · P1 permutes rows of A into order determined by partial pivoting, and now L is lower triangular Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 49 / 88 Complete Pivoting Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 50 / 88 Example: Pivoting Complete pivoting is more exhaustive strategy in which largest entry in entire remaining unreduced submatrix is permuted into diagonal pivot position Requires interchanging columns as well as rows, leading to factorization P AQ = L U Need for pivoting has nothing to do with whether matrix is singular or nearly singular For example, with L unit lower triangular, U upper triangular, and P and Q permutations Numerical stability of complete pivoting is theoretically superior, but pivot search is more expensive than for partial pivoting Numerical stability of partial pivoting is more than adequate in practice, so it is almost always used in solving linear systems by Gaussian elimination Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems is nonsingular yet has no LU factorization unless rows are interchanged, whereas is singular yet has LU factorization 51 / 88 Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Example: Small Pivots Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 52 / 88 Example, continued To illustrate effect of small pivots, consider Using small pivot, and correspondingly large multiplier, has caused loss of information in transformed matrix If rows interchanged, then pivot is 1 and multiplier is − , so 1 11 A= where is positive number smaller than mach If rows are not interchanged, then pivot is and multiplier is −1/ , so 1 0 10 M= , L= , −1/ 1 1/ 1 1 0 1 − 1/ = 0 1 1 0 −1/ U= 1 0 −1/ Michael T. Heath Scientific Computing 0 , 1 1 1 0 1− in floating-point arithmetic Thus, 10 LU = 1 1 = =A 10 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems 1 − M= in floating-point arithmetic, but then 1 LU = 1/ 11 11 A= Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy U= 01 10 A= L= = 10 , 1 11 01 11 11 = 01 1 which is correct after permutation Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 53 / 88 Pivoting, continued Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 54 / 88 Residual ˆ ˆ Residual r = b − Ax for solution x computed using Gaussian elimination satisfies Although pivoting is generally required for stability of Gaussian elimination, pivoting is not required for some important classes of matrices r A Diagonally dominant ≤ E ≤ ρ n2 A mach where E is backward error in matrix A and growth factor ρ is ratio of largest entry of U to largest entry of A n |aij | < |ajj |, ˆ x j = 1, . . . , n i=1, i=j Without pivoting, ρ can be arbitrarily large, so Gaussian elimination without pivoting is unstable Symmetric positive definite A = AT and xT Ax > 0 for all x = 0 Michael T. Heath Scientific Computing With partial pivoting, ρ can still be as large as 2n−1 , but such behavior is extremely rare 55 / 88 Michael T. Heath Scientific Computing 56 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Residual, continued Example: Small Residual Use 3-digit decimal arithmetic to solve There is little or no growth in practice, so r A Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy ˆ x ≤ E ≈n A 0.641 0.242 0.321 0.121 mach Gaussian elimination with partial pivoting yields triangular system which means Gaussian elimination with partial pivoting yields small relative residual regardless of conditioning of system 0.641 0.242 0 0.000242 Thus, small relative residual does not necessarily imply computed solution is close to “true” solution unless system is well-conditioned Michael T. Heath Scientific Computing x1 0.883 = x2 −0.000383 Back-substitution then gives solution ˆ x = 0.782 1.58 T Exact residual for this solution is −0.000622 ˆ r = b − Ax = −0.000202 Complete pivoting yields even smaller growth factor, but additional margin of stability usually is not worth extra cost Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems x1 0.883 = x2 0.442 Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 57 / 88 Example, continued Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 58 / 88 Implementation of Gaussian Elimination Gaussian elimination has general form of triple-nested loop Residual is as small as we can expect using 3-digit arithmetic, but exact solution is x = 1.00 1.00 for for T for aij = aij − (aik /akk )akj end end end so error is almost as large as solution Cause of this phenomenon is that matrix is nearly singular (cond(A) > 4000) Division that determines x2 is between two quantities that are both on order of rounding error, and hence result is essentially arbitrary When arbitrary value for x2 is substituted into first equation, value for x1 is computed so that first equation is satisfied, yielding small residual, but poor solution Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Indices i, j , and k of for loops can be taken in any order, for total of 3! = 6 different arrangements These variations have different memory access patterns, which may cause their performance to vary widely on different computers Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 59 / 88 Uniqueness of LU Factorization Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 60 / 88 Elementary elimination matrices Mk , their inverses Lk , and permutation matrices Pk used in formal description of LU factorization process are not formed explicitly in actual implementation Provided row pivot sequence is same, if we have two LU ˆˆ ˆ ˆ factorizations P A = LU = LU , then L−1 L = U U −1 = D is both lower and upper triangular, hence diagonal U overwrites upper triangle of A, multipliers in L overwrite strict lower triangle of A, and unit diagonal of L need not be stored ˆ If both L and L are unit lower triangular, then D must be ˆ ˆ identity matrix, so L = L and U = U Uniqueness is made explicit in LDU factorization P A = LDU , with L unit lower triangular, U unit upper triangular, and D diagonal Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Row interchanges usually are not done explicitly; auxiliary integer vector keeps track of row order in original locations Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 61 / 88 Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 62 / 88 Even with many right-hand sides b, inversion never overcomes higher initial cost, since each matrix-vector multiplication A−1 b requires n2 operations, similar to cost of forward- and back-substitution Forward- and back-substitution for single right-hand-side vector together require about n2 multiplications and similar number of additions Inversion gives less accurate answer; for example, solving 3x = 18 by division gives x = 18/3 = 6, but inversion gives x = 3−1 × 18 = 0.333 × 18 = 5.99 using 3-digit arithmetic Can also solve linear system by matrix inversion: x = A−1 b Matrix inverses often occur as convenient notation in formulas, but explicit inverse is rarely required to implement such formulas Computing A−1 is tantamount to solving n linear systems, requiring LU factorization of A followed by n forward- and back-substitutions, one for each column of identity matrix Scientific Computing Scientific Computing Inversion vs. Factorization LU factorization requires about n3 /3 floating-point multiplications and similar number of additions n3 , Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Complexity of Solving Linear Systems Michael T. Heath Scientific Computing Storage Management Despite variations in computing it, LU factorization is unique up to diagonal scaling of factors Operation count for inversion is about expensive as LU factorization Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems For example, product A−1 B should be computed by LU factorization of A, followed by forward- and back-substitutions using each column of B three times as 63 / 88 Michael T. Heath Scientific Computing 64 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Gauss-Jordan Elimination Gauss-Jordan Elimination, continued Gauss-Jordan elimination requires about n3 /2 multiplications and similar number of additions, 50% more expensive than LU factorization In Gauss-Jordan elimination, matrix is reduced to diagonal rather than triangular form Row combinations are used to annihilate entries above as well as below diagonal Elimination matrix used for given column vector a is of form 1 . . . 0 0 0 . . . 0 ··· .. . ··· ··· ··· .. . ··· 0 . . . −m 1 . . . 0 . . . 1 0 0 . . . 0 −mk−1 1 −mk+1 . . . −mn 0 0 1 . . . 0 ··· .. . ··· ··· ··· .. . ··· During elimination phase, same row operations are also applied to right-hand-side vector (or vectors) of system of linear equations 0 a1 0 . . . . . . . . . 0 ak−1 0 0 ak = ak ak+1 0 0 . . . . . . . . . 0 an 1 Once matrix is in diagonal form, components of solution are computed by dividing each entry of transformed right-hand side by corresponding diagonal entry of matrix Latter requires only n divisions, but this is not enough cheaper to offset more costly elimination phase < interactive example > where mi = ai /ak , i = 1, . . . , n Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 65 / 88 Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Solving Modified Problems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 66 / 88 Sherman-Morrison Formula Sometimes refactorization can be avoided even when matrix does change If right-hand side of linear system changes but matrix does not, then LU factorization need not be repeated to solve new system Sherman-Morrison formula gives inverse of matrix resulting from rank-one change to matrix whose inverse is already known Only forward- and back-substitution need be repeated for new right-hand side (A − uv T )−1 = A−1 + A−1 u(1 − v T A−1 u)−1 v T A−1 This is substantial savings in work, since additional triangular solutions cost only O(n2 ) work, in contrast to O(n3 ) cost of factorization Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems where u and v are n-vectors Evaluation of formula requires O(n2 ) work (for matrix-vector multiplications) rather than O(n3 ) work required for inversion Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 67 / 88 = A−1 b + A−1 u(1 − v T A−1 u)−1 v T A−1 b 68 / 88 (with 3, 2 entry changed) of system whose LU factorization was computed in earlier example One way to choose update vectors is 0 0 0 and v = 1 u= −2 0 which can be implemented by following steps −1 Solve Az = u for z , so z = A u Solve Ay = b for y , so y = A−1 b Compute x = y + ((v T y )/(1 − v T z ))z If A is already factored, procedure requires only triangular solutions and inner products, so only O(n2 ) work and no explicit inverses Scientific Computing Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Consider rank-one modification 2 4 −2 x1 2 4 9 −3 x2 = 8 −2 −1 7 x3 10 x = (A − uv T )−1 b Michael T. Heath Scientific Computing Example: Rank-One Updating of Solution To solve linear system (A − uv T )x = b with new matrix, use Sherman-Morrison formula to obtain Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Rank-One Updating of Solution so matrix of modified system is A − uv T Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 69 / 88 Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Example, continued Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 70 / 88 Scaling Linear Systems In principle, solution to linear system is unaffected by diagonal scaling of matrix and right-hand-side vector Using LU factorization of A to solve Az = u and Ay = b, −3/2 −1 1/2 and y = 2 z= −1/2 2 In practice, scaling affects both conditioning of matrix and selection of pivots in Gaussian elimination, which in turn affect numerical accuracy in finite-precision arithmetic Final step computes updated solution −1 −3/2 −7 2 vT y 2 + 1/2 = 4 x=y+ z= 1 − vT z 1 − 1/2 2 −1/2 0 It is usually best if all entries (or uncertainties in entries) of matrix have about same size Sometimes it may be obvious how to accomplish this by choice of measurement units for variables, but there is no foolproof method for doing so in general We have thus computed solution to modified system without factoring modified matrix Michael T. Heath Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Scientific Computing Scaling can introduce rounding errors if not done carefully 71 / 88 Michael T. Heath Scientific Computing 72 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Example: Scaling Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy Iterative Refinement Given approximate solution x0 to linear system Ax = b, compute residual r0 = b − Ax0 Linear system 10 0 x1 1 = x2 Now solve linear system Az0 = r0 and take has condition number 1/ , so is ill-conditioned if is small x1 = x0 + z0 If second row is multiplied by 1/ , then system becomes perfectly well-conditioned as new and “better” approximate solution, since Ax1 = A(x0 + z0 ) = Ax0 + Az0 Apparent ill-conditioning was due purely to poor scaling = ( b − r0 ) + r0 = b In general, it is usually much less obvious how to correct poor scaling Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Process can be repeated to refine solution successively until convergence, potentially producing solution accurate to full machine precision Triangular Systems Gaussian Elimination Updating Solutions Improving Accuracy 73 / 88 Iterative Refinement, continued Michael T. Heath Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems 74 / 88 Symmetric Systems Banded Systems Iterative Methods Special Types of Linear Systems Iterative refinement requires double storage, since both original matrix and its LU factorization are required Work and storage can often be saved in solving linear system if matrix has special properties Due to cancellation, residual usually must be computed with higher precision for iterative refinement to produce meaningful improvement Examples include Symmetric : A = AT , aij = aji for all i, j For these reasons, iterative improvement is often impractical to use routinely, but it can still be useful in some circumstances Positive definite : xT Ax > 0 for all x = 0 Band : aij = 0 for all |i − j | > β , where β is bandwidth of A Sparse : most entries of A are zero For example, iterative refinement can sometimes stabilize otherwise unstable algorithm Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Scientific Computing 75 / 88 Symmetric Systems Banded Systems Iterative Methods Symmetric Positive Definite Matrices l11 l21 0 l22 implies a11 , l21 = a21 /l11 , Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems l22 = 2 a22 − l21 Scientific Computing 77 / 88 Symmetric Systems Banded Systems Iterative Methods Scientific Computing 78 / 88 Symmetric Systems Banded Systems Iterative Methods Symmetric Indefinite Systems Features of Cholesky algorithm for symmetric positive definite matrices For symmetric indefinite A, Cholesky factorization is not applicable, and some form of pivoting is generally required for numerical stability All n square roots are of positive numbers, so algorithm is well defined No pivoting is required to maintain numerical stability Only lower triangle of A is accessed, and hence upper triangular portion need not be stored Only n3 /6 multiplications and similar number of additions are required Factorization of form P AP T = LDLT with L unit lower triangular and D either tridiagonal or block diagonal with 1 × 1 and 2 × 2 diagonal blocks, can be computed stably using symmetric pivoting strategy Thus, Cholesky factorization requires only about half work and half storage compared with LU factorization of general matrix by Gaussian elimination, and also avoids need for pivoting In either case, cost is comparable to that of Cholesky factorization < interactive example > Scientific Computing Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Cholesky Factorization, continued Michael T. Heath Symmetric Systems Banded Systems Iterative Methods for j = 1 to n for k = 1 to j − 1 for i = j to n aij = aij − aik · ajk end end √ ajj = ajj for k = j + 1 to n akj = akj /ajj end end where L is lower triangular with positive diagonal entries Algorithm for computing it can be derived by equating corresponding entries of A and LLT In 2 × 2 case, for example, √ 76 / 88 One way to write resulting general algorithm, in which Cholesky factor L overwrites original matrix A, is A = L LT l11 = Scientific Computing Cholesky Factorization If A is symmetric and positive definite, then LU factorization can be arranged so that U = LT , which gives Cholesky factorization a11 a21 l 0 = 11 a21 a22 l21 l22 Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems 79 / 88 Michael T. Heath Scientific Computing 80 / 88 Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Symmetric Systems Banded Systems Iterative Methods Band Matrices Symmetric Systems Banded Systems Iterative Methods Tridiagonal Matrices Consider tridiagonal matrix Gaussian elimination for band matrices differs little from general case — only ranges of loops change b1 a2 A=0 . . . Typically matrix is stored in array by diagonals to avoid storing zero entries If pivoting is required for numerical stability, bandwidth can grow (but no more than double) 0 For fixed small bandwidth, band solver can be extremely simple, especially if pivoting is not required for stability Michael T. Heath Scientific Computing 0 b2 .. . c2 .. . ··· .. . .. . .. . ··· an−1 0 bn − 1 an 0 . . . 0 cn−1 bn Gaussian elimination without pivoting reduces to d1 = b1 for i = 2 to n mi = ai /di−1 di = bi − mi ci−1 end General purpose solver for arbitrary bandwidth is similar to code for Gaussian elimination for general matrices Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems c1 81 / 88 Michael T. Heath Symmetric Systems Banded Systems Iterative Methods Tridiagonal Matrices, continued Scientific Computing Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems 82 / 88 Symmetric Systems Banded Systems Iterative Methods General Band Matrices LU factorization of A is then given by 1 m2 L= 0 . . . 0 0 1 .. . .. . ··· ··· .. . .. . mn−1 0 ··· .. . 1 mn 0 . . . . , . . 0 1 Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems d1 0 . U =. . . . . 0 c1 0 d2 .. . c2 .. . .. ··· . ··· ··· .. . .. . dn−1 0 In general, band system of bandwidth β requires O(βn) storage, and its factorization requires O(β 2 n) work 0 . . . 0 cn−1 dn Scientific Computing Compared with full system, savings is substantial if β 83 / 88 Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems Symmetric Systems Banded Systems Iterative Methods Iterative Methods for Linear Systems LAPACK is more recent replacement for LINPACK featuring higher performance on modern computer architectures, including some parallel computers For some types of problems, iterative methods have significant advantages over direct methods Both LINPACK and LAPACK are available from Netlib We will study specific iterative methods later when we consider solution of partial differential equations 85 / 88 Michael T. Heath Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems LINPACK and LAPACK BLAS Basic Linear Algebra Subprograms Level 1 Work O(n) BLAS encapsulate basic operations on vectors and matrices so they can be optimized for given computer architecture while high-level routines that call them remain portable 2 O(n2 ) 3 O(n3 ) Higher-level BLAS encapsulate matrix-vector and matrix-matrix operations for better utilization of memory hierarchies such as cache and virtual memory with paging Generic Fortran versions of BLAS are available from Netlib, and many computer vendors provide custom versions optimized for their particular systems Scientific Computing Scientific Computing 86 / 88 LINPACK and LAPACK BLAS Examples of BLAS High-level routines in LINPACK and LAPACK are based on lower-level Basic Linear Algebra Subprograms (BLAS) Michael T. Heath LINPACK and LAPACK BLAS Solving linear systems of such fundamental importance in scientific computing that LINPACK has become standard benchmark for comparing performance of computers In theory, it might take infinite number of iterations to converge to exact solution, but in practice iterations are terminated when residual is as small as desired Scientific Computing 84 / 88 LINPACK is software package for solving wide variety of systems of linear equations, both general dense systems and special systems, such as symmetric or banded Iterative methods begin with initial guess for solution and successively improve it until desired accuracy attained Michael T. Heath Scientific Computing LINPACK and LAPACK Gaussian elimination is direct method for solving linear system, producing exact solution in finite number of steps (in exact arithmetic) Existence, Uniqueness, and Conditioning Solving Linear Systems Special Types of Linear Systems Software for Linear Systems n Examples saxpy sdot snrm2 sgemv strsv sger sgemm strsm ssyrk Function Scalar × vector + vector Inner product Euclidean vector norm Matrix-vector product Triangular solution Rank-one update Matrix-matrix product Multiple triang. solutions Rank-k update Level-3 BLAS have more opportunity for data reuse, and hence higher performance, because they perform more operations per data item than lower-level BLAS 87 / 88 Michael T. Heath Scientific Computing 88 / 88 ...
View Full Document

{[ snackBarMessage ]}

Ask a homework question - tutors are online