程序代写代做代考 Fortran algorithm handoutA.dvi
handoutA.dvi
ECS130 Scientific Computation Handout A January 11, 2017
1. Review of vector and matrix algebra
A±B, A ·B, AT (= A′), A · x, xT · y, x · yT
where x and y are vectors, and A and B are matrices.
2. BLAS = Basic Linear Algebra Subprograms, a de facto application programming interface stan-
dard to perform basic linear algebra operations such as vector and matrix multiplication.
• Level 1: vector operations, such as y ← αx+ y, …
• Level 2: matrix-vector operations, such y ← αAx+ βy, ….
• Level 3: matrix-matrix operations, such C ← αAB + βC, ….
Highly optimized implementations of the BLAS have been developed by hardware vendors such
as by MKL of Intel, ESSL of IBM, cuBLAS of NVIDIA.
3. Linear system of equations
Ax = b
where A is a given square matrix of order n, b is a given column vector of n components, and x
is an unknown column vector of n components.
The following statements are equivalent:
• for any vector b, the linear system has a solution x.
• If a solution exists, it is unique.
• If Ax = 0, then x = 0.
• The columns (rows) of A are linearly independent.
• There is a matrix A−1 such that A−1 ·A = A ·A−1 = I.
(A is called nonsingular or invertible, A−1 is called the inverse of A)
• det(A) 6= 0.
4. Lower triangular linear system of equations:
Lx = b
where L = (ℓij) is an n×n lower triangular, i.e., ℓij = 0 if i < j, and nonsingular (invertible), i.e., ℓii 6= 0 for i = 1, . . . , n. 5. Forward substitution algorithm – row-oriented • Algorithm: for i = 1, 2, . . . , n, xi = (bi − li,1x1 − li,2x2 − · · · − li,i−1xi−1) /li,i • Flops: n2 • M-scripts in componentwise form: 1 for i = 1:n x(i) = b(i); for j = 1:i-1 x(i) = x(i) - L(i,j)*x(j); end x(i) = x(i)/L(i,i); end • M-scripts in vectorized form: x(1) = b(1)/L(1,1); for i = 2:n x(i) = (b(i) - L(i,1:i-1)*x(1:i-1))/L(i,i); end 6. Forward substitution algorithm – column-oriented • Algorithm: By the partition [ 1 n− 1 1 ℓ11 n− 1 L21 L22 ] [ x1 x(2:n) ] = [ b1 b(2:n) ] we have ℓ11x1 = b1 L21x1 + L22x(2:n) = b(2:n) Therefore, x1 = b1/ℓ11, and then after updating b̂(2:n) = b(2:n) − L21x1, we solve solve the same-kind of lower triangular system with dimension n− 1: L22x(2:n) = b̂(2:n). The procedure is repeated until finding all entries of x. • M-scripts in vectorized form: x = zeros(n,1); for j = 1:n-1 x(j) = b(j)/L(j,j); b(j+1:n) = b(j+1:n) - L(j+1:n,j)*x(j); end x(n) = b(n)/L(n,n); 7. Row-oriented vs. column-oriented forward substitution Language considerations: • C stores double subscripted arrays by rows, • Fortran stores by columns. Memory considerations • virtual memory • page-hit/miss 2 • locality of reference. 8. Triangular systems with multiple right-hand sides: LX = B, where B is n×m. Algorithm 1: solve Lxk = bk for k = 1 : m Algorithm 2: vectorized on m X = zeros(n,m); for j = 1:n-1 X(j,:) = B(j,:)/L(j,j); B(j+1:n,:) = B(j+1:n,:) - L(j+1:n,:)*X(j,:); end X(n,:) = B(n,:)/L(n,n); 9. Exercise: Upper triangular linear system of equations Ux = b can be treated similarly using back substitution, where U is an n×n upper triangular (i.e., uij = 0 for i > j) and nonsingular (i.e., uii = 0 for i = 1, 2, . . . , n)
3
