Dense Matrix Computations
 |
For a detailed description of MPI Capability see Dense Matrix Parallelism Using MPI. |
This section is concerned with methods for computing with dense matrices. Consider a Fortran 90 code fragment that solves a linear system of algebraic equations, Ay = b, then computes the residual r = b − Ay. A standard mathematical notation is often used to write the solution,
A user thinks: “matrix and right-hand side yields solution.” The code shows the computation of this mathematical solution using a defined Fortran operator “.ix.”, and random data obtained with the function, rand. This operator is read “inverse matrix times.” The residuals are computed with another defined Fortran operator “.x.”, read “matrix times vector.” Once a user understands the equivalence of a mathematical formula with the corresponding Fortran operator, it is possible to write this program with little effort. The last line of the example before end is discussed below.
USE linear_operators
integer,parameter :: n=3; real A(n,n), y(n), b(n), r(n)
A=rand(A); b=rand(b); y = A .ix. b
r = b - (A .x. y ) ! Parentheses are needed
end
The IMSL Fortran Numerical Library provides additional lower-level software that implements the operation “.ix.”, the function rand, matrix multiply “.x.”, and others not used in this example. Standard matrix products and inverse operations of matrix algebra are shown in the following table:
Defined Array Operation | Matrix Operation | Alternative in Fortran 90 |
|---|---|---|
A .x. B | matmul(A, B) | |
.i. A | lin_sol_gen lin_sol_lsq | |
.t. A, .h. A | transpose(A) conjg(transpose(A)) | |
A .ix. B | lin_sol_gen lin_sol_lsq | |
B .xi. A | lin_sol_gen lin_sol_lsq | |
A .tx. B, or (.t. A) .x. B A .hx. B, or (.h. A) .x. B | matmul(transpose (A), B) matmul(conjg(transpose(A)), B) | |
B .xt. A, or B .x. (.t. A) B .xh. A, or B .x. (.h. A) | matmul(B, transpose(A)) matmul(B, conjg(transpose(A))) |
The IMSL operators apply generically to all standard precisions and floating-point data types – real and complex – and to objects that are broader in scope than arrays with a fixed number of dimensions. For example, the matrix product “.x.” applies to matrix times vector and matrix times matrix represented as Fortran 90 arrays. It also applies to “independent matrix products.” For this, use the notion: a box of problems to refer to independent linear algebra computations, of the same kind and dimension, but different data. The racks of the box are the distinct problems. In terms of Fortran 90 arrays, a rank-3, assumed-shape array is the data structure used for a box. The first two dimensions are the data for a matrix; the third dimension is the rack number. Each problem is independent of other problems in consecutive racks of the box. We use parallelism of an underlying network of processors, and MPI, when computing these disjoint problems.
In addition to the operators .ix., .xi., .i., and .x., additional operators .t., .h., .tx., .hx., .xt., and .xh. are provided for complex matrices. Since the transpose matrix is defined for complex matrices, this meaning is kept for the defined operations. In order to write one defined operation for both real and complex matrices, use the conjugate-transpose in all cases. This will result in only real operations when the data arrays are real.
For sums and differences of vectors and matrices, the intrinsic array operations “+” and “‑” are available. It is not necessary to have separate defined operations. A parsing rule in Fortran 90 states that the result of a defined operation involving two quantities has a lower precedence than any intrinsic operation. This explains the parentheses around the next-to-last line containing the sub-expression “A .x. y” found in the example. Users are advised to always include parentheses around array expressions that are mixed with defined operations, or whenever there is possible confusion without them. The next-to-last line of the example results in computing the residual associated with the solution, namely r = b ‑ Ay. Ideally, this residual is zero when the system has a unique solution. It will be computed as a non-zero vector due to rounding errors and conditioning of the problem.
