In symbolic matrix algebra, we use expressions like C = A + B, C =
A*B and
. If an algorithm could be presented to a compiler
in this way, it would make the program easier to read, and easier to
maintain. The numerical package MatLab [] is already doing
this, and so, to some extent, is High Performance Fortran. With C++ it
is possible to implement this kind of syntax for almost any structure,
although one should be careful when doing so.
The mathematical operations of interest to matrix algebra are the
assignment (=), addition (+), subtraction (-) and multiplication
(*) operators. There are also the update operators
,
and
which can be used instead of C = C + A, C
= C - A and C = C * A respectively.
The division operator might have been used for solving equations, but
there is the problem of stability for this operation, which depends on
the values of the matrix. Such solving is further made difficult
because the solving may be done from the left side (
), or the
right side (
). Also, solving systems usually involves the
selection of a method suited to the problem, and inversions of a
matrix is both time consuming and may be numerically unstable. For
these reasons, it is better to specify the method explicitly in the
application. How SymPLA handles this is discussed in Section
. It might also be feasible to let division be
performed as an implicit QR solver step, or to let the algorithm
select a method as the default solver method.
The only division operator supported by SymPLA is used to divide
all elements of the matrix with a scalar value. The update division
is also supported for this method of scaling.
The arithmetic operations on matrices do, of course, check that the operations are made between operands of compatible dimensions, and, if that is not the case, report the error through the C++ exception handling system.
The transpose of a matrix,
, is often used in expressions, and a
number of libraries have used a syntax like A.t(),
A.trans() or transp(A) to get the transpose of the matrix
A. While this kind of function-call may or may not involve
duplication of data, this syntax gives the impression that to use the
transpose of a matrix you have to use a function with duplication of
data first. In LAPACK++, the use of implicit transpose is supported by
functions like Blas_Mat_Trans_Mat_Mult for multiplication
involving one transposed matrix (
), while in BLAS this is
specified with flags.
In SymPLA, the transpose of a matrix is considered to be an attribute
of the matrix, that may be accessed by A.T, which might
(theoretically) have been implemented with duplication of data by a
conversion operator, but instead the same class is handling transposed
and nontransposed matrices. A flag in the class is used to decide whether
the operation shall use the transposed or nontransposed matrix in its
execution, and optimization of the operation is left to the internal
routines in the library. Using this convention
can be written
A * B.T.
enum MatrixTranspose {NONTRANSPOSED, TRANSPOSED};
class Matrix;
class MatrixHandler{
private:
MatrixTranspose transposed;
public:
Matrix operator +(MatrixHandler &);
Matrix operator /(double);
MatrixHandler &operator =(MatrixHandler &);
MatrixHandler &operator +=(MatrixHandler &);
MatrixHandler &operator /=(double);
};
class Matrix : public MatrixHandler{
public:
MatrixHandler T;
Matrix &operator =(Matrix &);
Matrix &operator =(MatrixHandler &);
};
The skeleton class definition above shows how operations and handling of transposed matrices are implemented in this library.