Version 1.1 of ai05s/ai05-0047-1.txt
!standard G.3.1 07-04-04 AI05-0047-1/01
!class binding interpretation 07-04-04
!status work item 07-04-04
!status received 07-04-04
!subject Annoyances in the Generic_Real_Arrays package
I am using this AI to record a number of annoyances or bugs that I encountered during
the implementation of the Generic_Real_Arrays package.
1 - The multiplication of matrices is defined to "involve inner products" (G.3.1(56/2)).
In strict mode, this requires that each component of the result comply with the
requirements of G.3.1(83/2-84/2). Technically this is known as a "componentwise"
error bound. An algorithm that gives a componentwise error bound is necessarily
cubic (i.e., in O(N**3)). There exist however algorithms for multiplying matrices
that are sub-cubic (e.g., Strassen's method) but these algorithms give normwise
error bounds, meaning that the error on one component may spill over other components.
Some variants of BLAS and other widely-used linear algebra libraries use fast matrix
multiplication algorithms, so the accuracy requirements makes it impossible to
implement Generic_Real_Arrays by interfacing to these libraries.
2 - Generic_Real_Arrays doesn't provide an operation to transpose a matrix. Transposition
is an operation that shows up very often, particularly in conjunction with orthogonal
matrices (and therefore with eigensystems). True, implementing your own transposition
is not too hard, but the same could be said of many of the operations provided by this
3 - The definition of Eigensystem does not impose any constraint on the length of the
out parameters Values and Vectors. Instead it has the mysterious sentence "The index
ranges of the parameter Vectors are those of A". It is written as if Eigensystem had
a way to change the constraints of Vectors, which is evidently false. It would seem
that it should require that Values'Length, Vectors'Length(1) and Vectors'Length(2) be
all equal and equal to A'Length(1).
4 - For some matrices, the QR iteration used to compute the eigenvalues will just not
converge (no deflation will happen). While it is always possible to let it run forever,
it is typical in this case to give up after some number of iterations. The RM makes
no provision for raising an exception in this case. The values computed for (some of)
the eigenvalues may be really bogus, so it would be better to raise an exception than
to return garbage to the user.
5 - The index subtype for types Real_Vector and Real_Matrix is Integer. Presumably this was
intended to provide maximum flexibility in selecting the index range. It has however
unpleasant consequences. Consider:
Identity: constant Real_Matrix := ((1.0, 0.0), (0.0, 1.0));
Anyone trying to evaluate Identity'First - 1 won't like the result. Why would anyone
do that? Maybe to initialize a variable that will be used to iterate through the rows
or column of the matrix. This is sure to bite many users, and using a slightly narrower
index subtype would have been much wiser.
Questions? Ask the ACAA Technical Agent