Using LU Decomposition to Optimize the modconcen.m Routine

Using LU Decomposition to Optimize the modconcen.m Routine

Using LU Decomposition to Optimize the modconcen.m Routine Matt Tornowske April 1, 2002 Review of Gaussian Elimination

Find the augmented matrix of the system of linear equations. Find an echelon form of the augmented matrix using elementary row operations. Find the system of equations corresponding to the echelon form. Use back substitution to arrive at the solution.

Example using Gaussian Elimination Solve the following system of linear equations using the method of Gaussian Elimination: Starting with the augmented matrix create zeros below the pivot in the first column: Gaussian Elimination : Echelon Form At this stage we create a zero only below the pivot:

The corresponding system of equations is: Gaussian Elimination : Back Substitution We now solve the system using back substitution. Since we have x4 = 2 we substitute this into the row above: Substituting x4 = 2 and x3 = -5 into the first equation: And were left with:

Computational Cost of Gaussian Elimination for an n x n Matrix For multiplication: For addition: Total operations:

For 1532 x 1532 matrix: Method of LU Decomposition Find the LU decomposition of the matrix

Solve LY=B using forward substitution Solve UX=Y by back substitution Hopefully, we save computational time Finding the LU Factors

The LU Factorization is done by finding the Upper and Lower triangle matrices of an original matrix. The upper matrix, U, is the reduced echelon form of the original matrix, A. The lower matrix, L, is the inverse of the elementary matrices obtained from the operations to obtain U. LU decomposition: Getting

U Solve the following system of equations using LU decomposition: First reduce the associated matrix into upper triangular form: LU Decomposition: Elementary Matrices The elementary matrices that correspond to these row operations are: The inverses of these matrices are:

LU Decomposition: Getting L Thus, LU Decomposition: Solving LUX = B We now solve the given system LUX = B by solving the two subsystems LY = B and UX = Y: This lower triangular system has solution y1 = -1, y2 = 7, y3 = 2.

The solution to the upper triangular system is x1 = 5, x2 = -8, x3 = -1, which happens to be the solution to A. Computational Cost of LU Decomposition for an n x n Matrix

Total number of arithmetic operations to solve a system of equations using LU decomposition is exactly the same as Gaussian Elimination. The method is the most efficient method available for solving many systems all having the same matrix of coefficients, when the right hand sides are not all available in advance. Since the factorization only has to be done once, it uses (4n3 3n2 n)/6 operations. The solutions of the two triangular systems can be carried

out in 2n2 n operations compared to (4n3 + 9n2 7n)/6 operations for full Gaussian elimination each time. Ideally, for 1528 x 1528 matrix, this is a savings of nearly 99.8% per system! Permutation Matrix

Our LU decomposition wont work properly if there are any pivots needed We can use the lu() function in Matlab to create a permutation matrix P to keep track of these pivots Well simply change the algorithm slightly to solve LY=PB instead of LY=B Implementation in Matlab We use the Matlab lu() function to give us our L, U, and P matrices We will adjust the code in Jennifers

modified concentration routine to use LU factorization. Modified Concentration Routine with LU Implemented In the last outer for block (modifications in red): % do the LU factorization here. [L,U,P]=lu(w); clear w; for m=1:ntimes m

rhs = zeros(1528,1); for k = 1:1528 N = nb(k,find((nbor(k,:) == 1))); W = nb(k,find((nbor(k,:) == 3))); S = nb(k,find((nbor(k,:) == 5))); E = nb(k,find((nbor(k,:) == 7))); if(~isempty(N) & ~isempty(W) & ~isempty(S) & ~isempty(E)) Modified Routine, continued ux = scaledvelocity(k,1); uy = scaledvelocity(k,2);

rhs(k,1) = -concentration(k,m)/deltat; else rhs(k,1) = 0; end end %xjen=w\rhs; replace this line with %backward/forward substitution rhs=P*rhs; xjen2=L\rhs; xjen2=U\xjen2; Modified Routine, continued

concentration(:,m+1)=xjen2; end Note: We use the left matrix divide instead of inverse operation. While some improvement can Be gained by performing some inversions outside of the for loop, the increase in speed can be far offset in the larger memory requirements. Performance Gains on Pentium 4 1400 Mhz with 384 MB of RAM With deltat=0.5, old routine took 73 seconds while the new routine took

11 seconds Adjusting deltat for 0.25, old = 145 seconds, new = 17 seconds Deltat for .1, old = 350 seconds, new = 38 seconds Conclusions Depending on the timestep, weve created a speedup of 7-9 times the original routine Use LU decomposition, all the cool

people do!

Recently Viewed Presentations

  • Annex CHAPTER 1

    Annex CHAPTER 1

    ordre. public), the protection of public health or morals or the protection of the rights and freedoms of others. This article shall not prevent the imposition of lawful restrictions on members of the armed forces and of the police in...
  • What I wished I had learned in my TESOL program

    What I wished I had learned in my TESOL program

    Carolyn Kristj√°nsson,TWU & Nathan Kielstra, TWU. What I wished I had learned. in my TESOL program. TESL Canada Conference2012,Kamloops, BCOct. 11-13. ... - Assessment as its own module or elective course? Response. cont... Karen Rauser, UBCO.
  • Food Chains - Trinity University

    Food Chains - Trinity University

    Food Chain Song. Nature's food chain. Nature's food chain. Living things get energy. Through nature's food chain. Sung to the tune: Farmer in the Dell song 1st 4 slides. Plants Need the sunPlants Need the sunLiving things get energyThrough nature's...
  • Air Flow and Circulation

    Air Flow and Circulation

    Fronts After drawing isobars (and isotherms), we can now identify where low pressure, and high pressure, centers are located and where cold air is moving in, and areas where warm air is moving into Please remember: Wind is named by...
  • Word Problems Word Problems I think of a

    Word Problems Word Problems I think of a

    How many 1/2 litres in 5 1/2 litres? 5.5 √∑ 0.5 = 11 jugs Word Problems Dad bought a 2kg bag of carrots. He used 400grams to make a carrot cake. How many grams of carrots were left. Change the...
  • Title of Program: Medicine Grand Rounds Title of

    Title of Program: Medicine Grand Rounds Title of

    Ting "Christina" He. Weill Cornell Medicine. Jonathan Hingre. State University of New York Upstate Medical University. Arjun . Janardhan. Robert Larner, M.D., College of Medicine at the University of Vermont. Aurasch Moaven. University of Maryland School of Medicine. Alison Kim....
  • PROPERTIES OF MATTER - Valley Central High School

    PROPERTIES OF MATTER - Valley Central High School

    Extensive Properties: depends on how much or amount of substance. Ex. Mass, weight, volume, length. Intensive Properties: does not depend on the amount. Ex. Color, odor, density, Physical Change: size, shape, phase. Chemical Property/Change. The composition of the substance DOES...
  • Fixed Prosthodontics

    Fixed Prosthodontics

    Nonimpregnated Impregnated Final Impression Elastomeric impression materials are used to create these extremely accurate impressions. Mixing and application of light-bodied material around the prepared tooth. Mixing and loading of tray with heavy-bodied material to be seated onto the quadrant or...