• Contents
• Previous
• Next

A. FMS Theory

This appendix contains a mathematical description of the matrix algebra functions performed by FMS.

Introduction

Simultaneous Equations

The calculation required is to determine the values of N unknowns X(N) from the array of known coefficients A(N,N) and vector B(N). Because the time required to solve this system increases as the cube of the number of equations N, this calculation often dominates analysis programs, resulting in limitations being placed on the analysis resolution, N.

In discussing techniques used to solve simultaneous equations, it is convenient to use the following matrix notation:

[A]{X} = {B}

where

The symbol [A] denotes the N-by-N matrix of known coefficients. A single term in the matrix [A] is referred to as A_ij, where the first subscript i is the row and the second subscript j is the column. The set of N known values in {B} is called the right-hand side vector while the set of N unknown values represented by {X} is called the solution vector.

Systems of simultaneous equations are often used to solve boundary value problems in the fields of structural mechanics, thermal analysis, fluid mechanics, diffusion, and electromagnetics. For these applications, the vector {X} usually represents unknown values of the field variable at sample points distributed throughout the domain of the problem. The vector {B} contains the known boundary condition values at the sample points. The matrix [A] serves as a mathematical model of the domain, relating the known boundary conditions {B} to the unknown field variable values {X}. The number of sample points, N, generally determines the resolution of the analysis.

As an example, suppose a structural analysis is performed. The matrix [A], known as the stiffness matrix, serves as a mathematical model of the object being analyzed. The vector {B} contains loads that are applied to the object [A], resulting in the displacements {X}.

Several techniques are currently used for generating simultaneous equation models of boundary value problems. Among the most popular are finite elements, finite difference, boundary integrals and the moment method. Simultaneous equations can also be generated directly from a model of discrete elements, as occurs in power network and electrical circuit analysis.

Related Systems of Equations

FMS includes special features for parallel processing of multiple right-hand side and solution vectors that provide extremely fast processing of related systems of equations. Therefore, the symbols {X} and {B} are used throughout this manual to mean one or more vectors.

Linear and Nonlinear Systems of Equations

Matrix Symmetry and Data Type

Many boundary value problems produce symmetric matrices because the influence of one point i on another point j is the same as the influence of point j on point i. The property of symmetry is used in the equation solving algorithm to eliminate almost half of the matrix data storage and computation. The more general case of nonsymmetrical matrices, which provides for arbitrary coefficients A_ij, is also included.

The data type can be real or complex. Real data requires only a single floating-point value to represent each term, while complex data requires both the real and imaginary parts. In addition to requiring twice as much storage, complex data requires four times as much processing because cross-products between the real and imaginary parts must be evaluated. For example, it requires one multiply(*) and one add (+) to compute the following for real data:

However, for complex data, the calculation becomes the following:

Uppercase letters denote the real part of the data, lowercase letters denote the imaginary part of the data, and i is the square root of (-1). The real and imaginary parts of the data are computed as follows:

A = B + C*D - c*d
a = b + C*d + c*D

To optimize the use of memory pipelines, complex data is stored as the real part followed by the imaginary part. This storage is the default for variables typed COMPLEX in FORTRAN programs.

Profile Matrix Storage

[A] is stored by rows, starting at the first nonzero term and proceeding across the row up to, but not including, the diagonal. Similarly, [AU] is stored by columns, starting at the first nonzero term and proceeding down the column to, but not including, the diagonal. The diagonal is stored as a separate vector.

A profile vector {LOWEQ} is used, which completely describes the profile sparsity of the matrix. This vector contains the lowest coupled equation for each equation. For the lower triangle matrix [AL], this is the column number of the first nonzero term in each row. For the upper triangle matrix [AU], this is the row number of the first nonzero term in each column. The number of off-diagonal terms in row I of [AL] or column I of [AU] is I - LOWEQ(I). For a full matrix, all terms in the vector {LOWEQ} would be 1.

Nonsymmetric Matrices

FMS profile storage includes full matrices, which always start at row and column 1, and banded matrices, which start at a constant distance from the diagonal. As discussed during the factoring process, the matrix factors may always be overlaid on the profile storage described above.

Symmetric Matrices

Advantages of Profile Storage

Mathematical Advantages

Computational Advantages

Vector Inner-Products

In FORTRAN, the vector inner-product is implemented by a function as follows:

       DOUBLE PRECISION FUNCTION DOT (A, B, N)
       REAL*8 A(N), B(N)
       DOT = 0.0D0
       DO I=1,N
          DOT = DOT + A(I) * B(I)
       END DO
       RETURN
       END

This calculation, which uses unit address increments for the vectors {A} and {B}, achieves maximum memory performance. Several modern processors contain a multiply-accumulate instruction which implements the above loop in hardware.

Factoring

[A] = [B][C]

A typical term of the matrix [A], A_ij is computed from the inner-product of the i-th row of matrix [B] and the j-th column of matrix [C],

Nonsymmetric Factoring

[A] = [L][U]

To obtain an expression for computing a typical term U_ij in [U], consider its corresponding term A_ij in [A]. Based on the matrix multiplication process described above, the term A_ij is computed from the inner-product of the i-th row of [L] and the j-th column of [U],

The value of the first part of the inner-product is zero, because either or both of L_ik and U_kj will be outside the matrix profile. For this reason profile sparsity is used on large matrices that are to be factored. The terms L_ik and U_kj for k=1,(IJLOW-1) do not contribute to the matrix factors and remain zero during the factoring process.

The last part of the inner-product involving the sum from (i+1) to N will also be zero because L_ik is zero for terms in the upper triangle.

In computing the factors [L] and [U], there are N(N+1) terms to be determined from the N² terms in [A]. Therefore, N terms of [L] or [U] may be chosen arbitrarily. The choice that results in the most efficient algorithm is to select the diagonal terms of [L] to be 1. Using L_ii = 1, the above equation may be solved for U_ij as follows:

Terms in the Lower Triangular Matrix

[A] = [L][U]

Expressions for computing terms in the lower triangular factor [L] can be derived in an identical way. Consider a term A_ij that is below the diagonal,

The values of the first and fourth parts of the inner-product are zero as before because they involve terms outside the matrix profile or terms in the lower triangle of [U]. Solving the equation for the term L_ij gives the following:

Therefore, terms in the lower triangular factor [L] are also computed by subtracting an inner-product from the corresponding term in the original matrix [A]. However, the result is then multiplied by the diagonal factor reciprocal 1/U_jj.

During the factoring and vector solution processes, the diagonals of [U] are only used as denominators. Because multiplication is faster than division, the terms U_ii are inverted immediately after being computed and the reciprocal values 1/U_ii are overlaid on the original matrix diagonals A_ii in the diagonal vector.

Term-by-term Computation

Symmetric Factoring

[A] = [L][D][L]^T

(LD)_ik = L_ik D_kk

This inner-product can be divided into the following four parts:

The value of the first part of the inner-product is zero because (LD)_ik or L_jk is outside the matrix profile. The last part of the inner-product is zero because terms in the upper triangle of [L] are zero. Noting that L_jj = 1, the above equation can be solved to give the following:

The final values of L_ij for row i are obtained by multiplying the terms (LD)_ik by the diagonal factor reciprocals 1/D_kk.

L_ik = (LD)_ik (1/D_kk)

Because the diagonal factors [D] are only used as denominators, they are inverted and stored in the diagonal vector as reciprocals.

Before completing the calculation for L_ij the value in (LD)_ij must be used to compute the diagonal factor term D_ii. Letting j = i in the equation for (LD) is as follows:

(LD)_ii = L_ii D_ii

and that L_ii = 1, the equation for D_ii becomes the following:

       DOT = 0.0D0
       DO K = ILOW,(I-1)
          LDSAVE = LD(I,K)
          L(I,K) = DINV(K) * LD(I,K)
          DOT    = DOT + LDSAVE * L(I,K)
       END DO
       D(I,I) = A(I,I) - DOT
       DINV(I) = 1.0D0/D(I,I)

When computing [L] and [D], the order of computation must begin at the upper left corner of the matrix A₁₁ and end at the lower right A_nn. The only restriction is that terms required to compute each L_ij and D_ii be known from previous calculations. Unless additional storage is to be used, the calculation of the diagonal factor D_ii and the multiplication of the i-th row by the diagonal factor reciprocals must occur together, as discussed above. The actual FMS order is based on maximizing the reuse of high-speed registers once they have been loaded with data.

Vector Solution

[L][D][L]^T{X} = {B} (symmetric)

or

[L][U]{X} = {B} (nonsymmetric)

The first step is to solve the intermediate problem for the vector {Y} as follows:

[L]{Y} = {B}

For the data structure and algorithms in FMS, this is exactly the same for symmetric and nonsymmetric problems. In fact, if a symmetric matrix is processed with the nonsymmetric module, the matrix factor [L] is identical to the matrix factor [L] that would have been obtained by using the symmetric module.

An expression for computing a typical term in vector {Y}, Y_i, may be obtained from the vector calculation for B_i

Forward Reduction

FMS stores the matrix [L] by rows and the terms Y_k as a vector. Therefore the inner-product used for forward reduction produces incremental addressing that can be pipelined.

Diagonal Scaling

[D]{Z} = {Y}

Recalling that the diagonal factor [D] was stored in reciprocal form, this calculation becomes the following:

Z_i = (1/D_ii) Y_i

As each multiplication is performed, the values of the scaled vector Z_i are overlaid on the reduced vector Y_i.

Back Substitution in Nonsymmetric Problems

[U]{X} = {Y}

To obtain an expression for computing {X}, consider the vector calculation to compute Y_i, as follows:

(XU)_i = X_i U_ii

the above equation can be solved for the following:

X_i = (XU)_i (1/U_ii)

As the calculations proceed in a backwards direction of decreasing i, the terms in X_k are known when they are required to compute X_i. For this reason, this calculation is usually called back substitution.

The inner-product

X_i = (XU)_i (1/U_ii)

Back Substitution in Symmetric Problems

Matrix and Vector Multiplication

{C} = {D} - [A]{B}
{C} =       [A]{B}
{C} = {D} + [A]{B}

{E} = {B} - [A]{X}

The resulting error can then be used in interative refinement to compute a correction to the solution {X}.

A similar set of multiply routines are provided which use submatrix files instead of the assembled matrix:

{C} = {D} - [Si]{B}
{C} =       [Si]{B}
{C} = {D} + [Si]{B}

FMS also contains the following multiplication routines:

{C} = {D} - {A}{B}
{C} =       {A}{B}
{C} = {D} + {A}{B}

{C} = {D} - {A}T{B}
{C} =       {A}T{B}
{C} = {D} + {A}T{B}

[C] = [D] - {A}T{B}
[C] = [D]   {A}T{B}
[C] = [D] + {A}T{B}

[F] = {X}^T{Y}

where [F] is a full matrix and {X} and {Y} are FMS vector files containing one or more vectors. For symmetric problems, only the lower triangle of [F] is computed.

These multiply operations have several practical applications in computing Eigenvalues of the matrix [A]. For example, if {X} is a known Eigenvector, the corresponding Eigenvalue can be computed from the following:

e = ({X}^T[A]{X})/({X}^T{X})

This calculation uses the matrix-vectors multiply operation once and the vectors-vectors multiply operation twice.

These multiply operations are also used to project the matrix [A] into a subspace [A^*] where its Eigenvalues can be computed more efficiently. This projection involves the following computation:

[A^*] = {X}^T[A]{X}

where {X} are the vectors that define the subspace.

In the special case where the matrix is diagonal [D], FMS directly computes the following quadratic form:

[F] = {X}^T[D]{X} (symmetric)

[F] = {X}^T[D]{Y} (nonsymmetric)

FMS also includes the vectors matrix multiply calculation

{Y} = {X}[F]

This calculation may be used to project subspace eigenvectors back into the problem space or to compute new vectors {Y} from a linear combination of old vectors {X}.