Efficient Numeric Linear Algebra Computation

Computations on large Matrices and Vectors that contain floatingpoint data  both hardware float data and arbitrary precision software float data  can be performed very efficiently in Maple by taking advantage of a builtin library of numeric linear algebra routines. A number of these routines are provided through the alliance between Maplesoft and the Numerical Algorithms Group (NAG). In addition, parts of the CLAPACK and optimized ATLAS libraries have been integrated.


The selection of a subroutine to carry out a requested operation in the LinearAlgebra package is based on several factors:

–

The type of data in the input objects

–

The shape, storage, and order options of the input object(s)

–

The setting of the UseHardwareFloats environment variable

–

The availability of a NAG routine to complete the requested operation


These selection criteria are detailed in the following sections.


The Type of Data



Maple only considers looking for a NAG routine to carry out a particular LinearAlgebra routine if the input objects (Matrices, Vectors, or scalars) consist of only numeric data (i.e., the scalars (if any) and the entries in any Matrix or Vector inputs must be of type complex(extended_numeric)). Also, at least one of the numeric entries or scalars must be a floatingpoint number. To determine if these properties are satisfied for Matrix or Vector inputs, Maple first checks the datatype option on Matrix and Vector inputs. If this option indicates that the data is some form of float (including complex float), no further investigation of the corresponding entries is done. If the datatype option does not provide this information, the entries are inspected.


If it is determined that all the input objects are comprised of only numeric data, and at least one value is a float, copies are made of all inputs which do not have a datatype specifying floatingpoint entries, with each entry converted to the nearest floatingpoint value (as determined by Digits in the software floatingpoint case) as the copy is made. While this operation is fast, it is clearly overhead that can be avoided by ensuring that the input objects have the appropriate (float) datatype, as specified by the datatype option to the Matrix or Vector constructor when the objects are created.



Specifying the Shape, Storage, and Order Options



NAG routines expect input Matrices to be stored in Fortran (column major) order. This is the default for Matrices in Maple, but it can be overridden by using the order option when creating a Matrix. If necessary, an input Matrix will be copied and put into Fortran order before a NAG routine is called.


In some cases, several different NAG routines are available to carry out a particular LinearAlgebra operation, where the distinctions are based on the data structure used to represent the input Matrix or Vector objects. The intention here is that if a large subset of the entries in a Matrix or Vector are known to always be 0, then these entries do not have to be stored, and a more compact data structure can be used in place of the standard rectangular structure.


The Matrix and Vector constructors take options shape and storage, which allow you to tailor your Matrices and Vectors to the optimal NAG routine for the task at hand. In general, Maple chooses a storage that is compatible with whatever shape you specify, so specifying a shape is usually sufficient. See Matrix and Vector for details. (See also MatrixOptions and VectorOptions.)


If a NAG routine is not available to carry out the requested operation on the input Matrix and Vector objects with their given shapes and storages, but a NAG routine is available for the operation when applied to more generally shaped or stored Matrix or Vector objects, then copies will be made as necessary and that more general NAG routine will be invoked to carry out the operation. (Since the importing of NAG routines into Maple is an ongoing project, future releases will include the ability to handle a larger selection of input shapes and storage modes.)



The FloatingPoint Computation Environment



The UseHardwareFloats environment variable is a Boolean flag that controls the environment in which floatingpoint computation is carried out. While it will eventually control this environment selection for all floatingpoint computation, it is only used to choose between the hardwareprecision NAG library (UseHardwareFloats = true) and the arbitraryprecision NAG library (UseHardwareFloats = false), and to set the specific value of the datatype option on Matrices and Vectors when they are created with datatype = float. The arbitrary precision NAG library is a specially built version that is capable of working directly with Maple software's floatingpoint numbers.


If a Matrix or Vector is created (by a call to the appropriate constructor) with the option datatype = float, the UseHardwareFloats flag is checked to see if the Matrix or Vector should have hardware float entries (UseHardwareFloats = true) or software float entries (UseHardwareFloats = false). In the former case, the Matrix or Vector is created with datatype = float[8], meaning that the entries are 8byte hardware floats, while in the latter case it is created with datatype = sfloat ("software float").


If the inputs to a LinearAlgebra operation are all numeric, with at least one value being a float (so that the datatype check passes), and a NAG routine is available to carry out the task, then the UseHardwareFloats flag determines which of the hardwareprecision or arbitraryprecision versions of the NAG routine is used. Note:The datatypes of the input objects and the datatype of the result object are not considered here; the computation environment is controlled solely by the UseHardwareFloats flag. This means, for example, that if A and B are Matrices with datatype = sfloat, Digits = 100 and UseHardwareFloats = true, then the computation A+B is done by a hardwareprecision NAG routine, returning a result that is accurate only to, at most, hardware precision (approximately 15 base10 digits).


As mentioned above, setting a datatype value in the outputoptions parameter of the routine being invoked does not override the UseHardwareFloats flag. However, if this value does not agree with the UseHardwareFloats flag, then a final copy operation is required to produce the object that is returned.



NAG Routine Availability



The incorporation of the NAG library into Maple is an ongoing project. At least one generalpurpose NAG routine is available for every LinearAlgebra routine that warrants it, and many LinearAlgebra routines can dispatch to one of several NAG routines, depending on the particulars of the input objects for the operation.


See the section ExecutionTime Feedback for information on how to determine if a NAG routine is used to complete a LinearAlgebra computation.



ExecutionTime Feedback



Use the infolevel table to request user feedback messages from the LinearAlgebra routines.


If infolevel[LinearAlgebra] is 1 or higher, a userinfo message will be produced whenever an attempt is made to use a NAG routine to carry out a LinearAlgebra computation. The name of the NAG routine which is invoked will be preceded by "hw_" or "sw_". The "hw_" means that the hardware precision NAG routine is called. The "sw_" means that the arbitrary precision ("software") NAG routine is called. If for some reason the NAG routine is unable to complete this computation, a userinfo message to this effect is displayed.


If infolevel[LinearAlgebra] is 2 or higher, a userinfo message is produced whenever a copy operation is required to change any of the shape, storage, order or datatype of an input or output Matrix or Vector.



More on Shape and Storage



In general, LinearAlgebra functions use NAG routines to perform computations on input objects with datatype=sfloat or datatype=float[8] and which have only the default shape (none) and storage (rectangular). Some LinearAlgebra routines can also use NAG routines for input objects with other shapes (for example, triangular[upper]) and compatible storage (for example, if shape=triangular[upper] and storage=triangular[upper] or storage=rectangular, it may be that a NAG routine can be used without a copy operation being required), or for input objects with sparse storage.


If an input object has an unrecognized shape, or more than one shape, then the LinearAlgebra routines will copy the object to one with a format for which there is a NAG routine available (if there is one). This is because the NAG routine receives only the underlying physically stored data; it does not receive the shape information of the object. (The shape information is stored in a Maple routine called an indexing function, and this indexing function is not passed to the NAG routine.)


When more than one shape, or an unrecognized shape (for example, a userdefined shape) is attached to the input object, Maple assumes that the underlying physically stored data cannot be properly interpreted without the assistance of that shape information.


Extending the LinearAlgebra capabilities by adding NAG support for a wider range of input objects is an ongoing project.



Packed Formats



Some of the routines in the LinearAlgebra package can return objects that are stored in special "packed" formats. For example, LinearAlgebra:LUDecomposition(..) can return the factored form of its input in the form of a Vector (the pivots) and a single Matrix (containing both the "L" and "U" factors together). This is the form in which the underlying NAG routine returns its results. Therefore, in the case of floatingpoint factorizations (where the NAG routine is called as a subroutine), and where the factorization is used in subsequent operations by NAG routines, it is most efficient to request this packed form.


For example, suppose that A is a square Matrix and b is a Vector of matching dimension. Solving the system A . x = b can be completed by a single call to LinearSolve.

>

A := <<2.2,4.0,6.1><0,1.8,3><3.2,5,7.4>>;

${A}{:=}\left[\begin{array}{ccc}{2.2}& {0}& {}{3.2}\\ {4.0}& {1.8}& {}{5}\\ {}{6.1}& {3}& {7.4}\end{array}\right]$
 (1) 
${b}{:=}\left[\begin{array}{c}{2.1}\\ {5.9}\\ {3}\end{array}\right]$
 (2) 
$\left[\begin{array}{c}{}{1.79059829059829}\\ {2.01442307692308}\\ {}{1.88728632478632}\end{array}\right]$
 (3) 

However, if preserving the factored form for later use is required (for example, to solve a set of systems having the same coefficient Matrix but different right hand sides), then solving the system can be completed by the following sequence of steps.

>

p, lu := LUDecomposition( A, output='NAG' );

${p}{\,}{\mathrm{lu}}{:=}\left[\begin{array}{r}{3}\\ {2}\\ {3}\end{array}\right]{\,}\left[\begin{array}{ccc}{}{6.10000000000000}& {3.}& {7.40000000000000}\\ {}{0.655737704918033}& {3.76721311475410}& {}{0.147540983606556}\\ {}{0.360655737704918}& {0.287206266318538}& {}{0.488772845953002}\end{array}\right]$
 (4) 
>

x := LinearSolve( [p, lu], b );

${x}{:=}\left[\begin{array}{c}{}{1.79059829059829}\\ {2.01442307692308}\\ {}{1.88728632478632}\end{array}\right]$
 (5) 
>

x2 := LinearSolve( [p, lu], <3.1,1.3,.27> );

${\mathrm{x2}}{:=}\left[\begin{array}{c}{}{6.56666666666667}\\ {0.0833333333333346}\\ {}{5.48333333333334}\end{array}\right]$
 (6) 

By preserving the packed form of the factorization, unnecessary copy operations and repeated factorizations are not performed in this sequence of computations.






