Engineering
Computation
Lecture 6
E. T. S. I. Caminos, Canales y Puertos
1
Gauss elimination (operation counting)
Forward Elimination:
DO k = 1 to n–1
DO i = k+1 to n
r = A(i,k)/A(k,k)
DO j = k+1 to n
A(i,j)=A(i,j) – r*A(k,j)
ENDDO
B(i) = B(i) – r*B(k)
ENDDO
ENDDO
E. T. S. I. Caminos, Canales y Puertos
2
Gauss elimination (operation counting)
Operation Counting for Gaussian Elimination
Back Substitution:
X(n) = B(n)/A(n,n)
DO i = n–1 to 1 by –1
SUM = 0
DO j = i+1 to n
SUM = SUM + A(i,j)*X(j)
ENDDO
X(i) = [B(i) – SUM]/A(i,i)
ENDDO
E. T. S. I. Caminos, Canales y Puertos
3
Gauss elimination (operation counting)
Operation Counting for Gaussian Elimination
Forward Elimination
n
Inner loop:

1 = n - (k + 1 ) + 1 = n - k
j= k +1
n
Second loop:

 2 + (n - k )    ( 2  n )  k  ( n  k )
i = k +1
= (n2 + 2n) – 2(n + 1)k + k2
E. T. S. I. Caminos, Canales y Puertos
4
Gauss elimination (operation counting)
Operation Counting for Gaussian Elimination
Forward Elimination (cont'd)
n 1
Outer loop
 [( n  2 n )  2( n  1) k  k ]
2
=
2
k 1
n 1
 (n
2
 2n )

n 1
1  2 ( n  1)
k 1
 (n
2

k 1
 2 n )( n  1)  2 ( n  1)
n 1
k

k
2
k 1
( n  1))( n )
2

(n  1)(n )(2 n  1)
6
E. T. S. I. Caminos, Canales y Puertos
n3
=
+ O (n 2 )
3
5
Gauss elimination (operation counting)
Operation Counting for Gaussian Elimination
Back Substitution
n
Inner Loop:
 1  n  (i  1)  1  n - i
j i  1
Outer Loop:
n 1
n 1
n 1
i 1
i 1
i 1
 1  ( n  i)   (1  n )  1   i
 (1  n )( n  1) 
( n  1) n
2
=
n
2
+ O (n )
2
E. T. S. I. Caminos, Canales y Puertos
6
Gauss elimination (operation counting)
Total flops = Forward Elimination + Back Substitution
=
n3/3 + O (n2)
+ n2/2 + O (n)

n3/3 + O (n2)
To convert (A,b) to (U,b') requires n3/3, plus terms of order n2 and smaller,
flops.
To back solve requires:
1 + 2 + 3 + 4 + . . . + n = n (n+1) / 2 flops;
Grand Total: the entire effort requires n3/3 + O(n2) flops altogether.
E. T. S. I. Caminos, Canales y Puertos
7
Gauss-Jordan Elimination
Diagonalization by both forward and backward elimination in
each column.  a 11 a 12 a 13 ... a 1n   x 1   b 1 

a
 21
 a 31


a
 n1
a 22
a 23
...
a 32
a 33
...
...
a n2
a n3
...


x2




a 3n   x 3  =




a nn  
 x nn 

a 2n
 
b
 2
 
 b3 
 
 

bn 

Perform elimination both backwards and forwards until:
1

0

0


 0
0
0
1
0
0
1
0
0
0

0

0


1 
 x1 
 x1 




x2
x2








=
x
x
 3
 3










xn 

xn 

3
Operation count for Gauss-Jordan is: n  O ( n 2 )
(slower than Gauss elimination)
2
E. T. S. I. Caminos, Canales y Puertos
8
Gauss-Jordan Elimination
Example (two-digit arithmetic):
50
 1

 2
1
40
6
1
0

0
2
4
30
1
2

3
0
1
0
0 .0 3 8
0 .1
29
x1 = 0.015
x2 = 0.041
x3 = 0.093
E. T. S. I. Caminos, Canales y Puertos
1
0

0
0 .0 2
40
6
0 .0 1 9 
0 .0 5 

2 .7 
0 .0 2 
2

3
0 .0 4
4
30
1
0

0
0
1
0
0
0
1
0 .0 1 5 
0 .0 4 1 

0 .0 9 3 
(vs. 0.016, et = 6.3%)
(vs. 0.041, et = 0%)
(vs. 0.091, et = 2.2%)
9
Gauss-Jordan Matrix Inversion
The solution of: [A]{x} = {b} is:
{x} = [A]-1{b}
where [A]-1 is the inverse matrix of [A]
Consider:
[A] [A]-1 = [ I ]
1) Create the augmented matrix: [ A | I ]
2) Apply Gauss-Jordan elimination:
==> [ I | A-1 ]
E. T. S. I. Caminos, Canales y Puertos
10
Gauss-Jordan Matrix Inversion
Gauss-Jordan Matrix Inversion (with 2 digit arithmetic):
 A I  =
1
0
0   0


1
0
50
 1

 2
1
40
6
1

0

 0
0
0 .0 3 8
0 .0 2
 0 .0 0 5
1
0 .1
 0 .0 0 0 5
0 .0 2 5
0
28
 0 .0 3 7
 0 .1 5
1

0

 0
2
4
30
1
0
0
0
1
0
0 .0 2
40
6
0
0
0 .0 2
 0 .0 0 0 2 9
1
0
 0 .0 0 0 3 7
0 .0 2 6
0
1
 0 .0 0 1 3
 0 .0 0 5 4
0 .0 4
4
30
0 .0 2
-0 .0 2
-0 .0 4
0
1
0
0
0

1
0

0

1 
 0 .0 0 1 4 

 0 .0 0 3 6

0 .0 3 6 
MATRIX INVERSE [A-1]
E. T. S. I. Caminos, Canales y Puertos
11
Gauss-Jordan Matrix Inversion
CHECK:
[ A ] [ A ]-1 = [ I ]
50
2
 2
1
40
6
2 
4 
30 

 0 .0 2 0
 -0 .0 0 0 3 7
 -0 .0 0 1 3
-0 .0 0 0 2 9
0 .0 2 6
-0 .0 0 5 4
-0 .0 0 1 4 
-0 .0 0 3 6  
0 .0 3 6 

 0 .9 9 7
 0 .0 0 0
  0 .0 0 1
 0 .1 3
1 .0 1 6
 0 .0 1 2
 0 .0 0 2 
 0 .0 0 1 
1 .0 5 6 

[ A ]-1 { b } = { x }
 0 .0 2 0
 -0 .0 0 0 3 7
 -0 .0 0 1 3
-0 .0 0 2 9
0 .0 2 6
-0 .0 0 5 4
-0 .0 0 1 4 
-0 .0 0 3 6 
0 .0 3 6 

 1   0 .0 1 5 
 2    0 .0 3 3 
 3   0 .0 9 9 
x tru e
 0 .0 1 6 
  0 .0 4 1 
 0 .0 9 1 
 0 .0 1 6 
x   0 .0 4 0 
 0 .0 9 3 
Gaussian
Elimination
E. T. S. I. Caminos, Canales y Puertos
12
LU decomposition
•
LU decomposition - The LU decomposition is a method that uses the
elimination techniques to transform the matrix A in a product of triangular
matrices. This is specially useful to solve systems with different vectors b,
because the same decomposition of matrix A can be used to evaluate in an
efficient form, by forward and backward sustitution, all cases.
A =
E. T. S. I. Caminos, Canales y Puertos
L U
13
LU decomposition
A =
L U
Decomposition
U X =
D
A X =
B
Initial system
L U X =
B
Transformed system 1
Substitution
L D =
B
Transformed system 2
X
Backward sustitution
D
Forward sustitution
E. T. S. I. Caminos, Canales y Puertos
14
LU decomposition
•
LU decomposition is very much related to Gauss method, because the upper
triangular matrix is also looked for in the LU decomposition. Thus, only the
lower triangular matrix is needed.
•
Surprisingly, during the Gauss elimination procedure, this matrix L is obtained,
but one is not aware of this fact. The factors we use to get zeroes below the
main diagonal are the elements of this matrix L.
Substract
E. T. S. I. Caminos, Canales y Puertos
15
LU decomposition
E. T. S. I. Caminos, Canales y Puertos
16
LU decomposition (Complexity)
Basic Approach
Consider [A]{x} = {b}
a) Gauss-type "decomposition" of [A] into [L][U]
n3/3 flops
[A]{x} = {b} becomes [L] [U]{x} = {b}; let [U]{x}  {d}
b) First solve [L] {d} = {b} for {d} by forward subst.
n2/2 flops
c) Then solve [U]{x} = {d} for {x} by back substitution
n2/2 flops
E. T. S. I. Caminos, Canales y Puertos
17
LU Decompostion: notation
 a11
L =  a
 ij
0 
a n n 
 a11
U =  0

 1
L
=
 1  a
 ij
0
1 
1
U
=
 1  0

a ij 
1 
 0
L 0  = a
 ij
0
0 
0
 U 0  = 0

a ij 
0 
1
D
=
  0

0
1 
E. T. S. I. Caminos, Canales y Puertos
a ij 
a n n 
[A] = [L] + [U0]
[A] = [L0] + [U]
[A] = [L0] + [U0] + [D]
[A] = [L1] [U]
[A] = [L] [U1]
18
LU decomposition
LU Decomposition Variations
Doolittle [L1][U]
General [A]
Crout
[L][U1]
General [A]
Cholesky [L][L] T
Pos. Def. Symmetric [A]
Cholesky works only for Positive Definite Symmetric matrices
Doolittle versus Crout:
• Doolittle just stores Gaussian elimination factors where Crout
uses a different series of calculations (see C&C 10.1.4).
• Both decompose [A] into [L] and [U] in n3/3 FLOPS
• Different location of diagonal of 1's
• Crout uses each element of [A] only once so the same array
can be used for [A] and [L\U] saving computer memory!
E. T. S. I. Caminos, Canales y Puertos
19
LU decomposition
Matrix Inversion
Definition of a matrix inverse:
[A] [A]-1 = [ I ]
==> [A] {x} = {b}
[A]-1 {b} = {x}
First Rule:
Don’t do it.
(numerically unstable calculation)
E. T. S. I. Caminos, Canales y Puertos
20
LU decomposition
Matrix Inversion
If you really must -1) Gaussian elimination: [A | I ] –> [U | B'] ==> A-1
2) Gauss-Jordan: [A | I ] ==> [I | A-1 ]
Inversion will take
n3 + O(n2)
flops if one is careful about where zeros are (taking advantage of
the sparseness of the matrix)
Naive applications (without optimization) take 4n3/3 + O(n2) flops.
For example, LU decomposition requires n3/3 + O(n2) flops. Back
solving twice with n unit vectors ei:
2 n (n2/2) = n3 flops.
Altogether: n3/3 + n3 = 4n3/3 + O(n2) flops
E. T. S. I. Caminos, Canales y Puertos
21
FLOP Counts for Linear Algebraic Equations
Summary
FLOP Counts for Linear Algebraic Equations, [A]{x} =
{b}
Gaussian Elimination (1 r.h.s)
n3/3 + O (n2)
Gauss-Jordan (1 r.h.s)
n3/2 + O (n2)
LU decomposition
n3/3 + O (n2)
Each extra LU right-hand-side
n2
Cholesky decomposition (symmetric A) n3/6 + O (n2)
Inversion (naive Gauss-Jordan)
4n3/3 +O (n2)
Inversion (optimal Gauss-Jordan)
Solution by Cramer's Rule
E. T. S. I. Caminos, Canales y Puertos
n3 + O (n2)
n!
22
Descargar

Lecture 6 (10/20)