• Choose a category
• All
• Classic-Fiction

Figure 14-1. LU factorization

A tool often used to teach Gaussian elimination is MATLAB. It features a scripting language (also called MATLAB) that makes developing matrix algorithms very simple. The language might seem very unusual to people familiar with other scripting languages because it is oriented to process multidimensional arrays. The unique features of the language that we use in the example code are:

Transposition operator for vectors and matrices: ' (single quote)

Matrix indexing specified as:

— Simple integer values: A(m, k)

— Ranges: A(k:n, k)

— Other matrices: A([k m], : )

Built-in matrix functions such as size (returns matrix dimensions), tril (returns the lower triangular portion of the matrix), triu (returns the upper triangular portion of the matrix), and eye (returns an identity matrix, which contains only zero entries, except for the diagonal, which is all ones)

Example 14-1 shows the simple implementation.

Example 14-1. Simple variant (MATLAB coding)

Code View: Scroll / Show All

function [L,U,p] = lutx(A)

%LUTX Triangular factorization, textbook version

% [L,U,p] = lutx(A) produces a unit lower triangular matrix L,

% an upper triangular matrix U, and a permutation vector p,

% so that L*U = A(p,:)

[n,n] = size(A);

p = (1:n)';

for k = 1:n-1

% Find index 'm' of largest element 'r' below diagonal in k-th column

[r,m] = max(abs(A(k:n,k)));

m = m+k-1; % adjust 'm' so it becomes a global index

% Skip elimination if column is zero

if (A(m,k) ~= 0)

% Swap pivot row

if (m ~= k)

A([k m],:) = A([m k],:); % swap rows 'k' and 'm' of 'A'

p([k m]) = p([m k]); % swap entrix 'k' and 'm' of 'p'

end

% Compute multipliers

i = k+1:n;

A(i,k) = A(i,k)/A(k,k);

% Update the remainder of the matrix

j = k+1:n;

A(i,j) = A(i,j) - A(i,k)*A(k,j);

end

end

% Separate result

L = tril(A,-1) + eye(n,n);

U = triu(A);

The algorithm presented in Example 14-1 is row-oriented, in the sense that we are taking a scalar multiple of the "pivot" row and adding it to the rows below to introduce zeros below the diagonal. The beauty of the algorithm lies in its similarity to the mathematical notation. Hence, this is the preferred way of teaching the algorithm for the first time so that students can quickly turn formulas into running code.

This beauty, however, has its price. In the 1970s, Fortran was the language for scientific computations. Fortran stores two-dimensional arrays by column. Accessing the array in a row-wise fashion within the matrix could involve successive memory reference to locations separated from each other by a large increment, depending on the size of the declared array. The situation was further complicated by the operating system's use of memory pages to effectively control memory usage. With a large matrix and a row-oriented algorithm in a Fortran environment, an excessive number of page swaps might be generated in the process of running the software. Cleve Moler pointed this out in the 1970s (see "Further Reading" at the end of this chapter).

To avoid this situation, one needed simply to interchange the order of the innermost nested loops on i and j. This simple change resulted in more than 30 percent savings in wall-clock time to solve problems of size 200 on an IBM 360/67. Beauty was thus traded for efficiency by using a less obvious ordering of loops and a much more obscure (by today's standard) language.

How Elegant Code Evolves with Hardware The Case of Gaussian Elimination > LINPACK's DGEFA Subroutine

14.4. LINPACK's DGEFA Subroutine

The performance issues with the MATLAB version of the code continued as, in the mid-1970s, vector architectures became available for scientific computations. Vector architectures exploit pipeline processing by running mathematical operations on arrays of data in a simultaneous or pipelined fashion. Most algorithms in linear algebra can be easily vectorized. Therefore, in the late 1970s there was an effort to