Efficient Numeric Linear Algebra Computation
Computations on large Matrices and Vectors that contain floating-point data -- both hardware float data and arbitrary precision software float data -- can be performed very efficiently in Maple by taking advantage of a built-in 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
Specifying the Shape, Storage, and Order Options
The Floating-Point Computation Environment
NAG Routine Availability
More on Shape and Storage
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 floating-point 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 floating-point entries, with each entry converted to the nearest floating-point value (as determined by Digits in the software floating-point 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.
See the section Execution-Time Feedback for information on how to determine if any copy operations were required.
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.)
See the section Execution-Time Feedback for information on how to determine if copy operations were required.
The UseHardwareFloats environment variable is a Boolean flag that controls the environment in which floating-point computation is carried out. While it will eventually control this environment selection for all floating-point computation, it is only used to choose between the hardware-precision NAG library (UseHardwareFloats = true) and the arbitrary-precision 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 floating-point 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, meaning that the entries are 8-byte 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 hardware-precision or arbitrary-precision 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 hardware-precision NAG routine, returning a result that is accurate only to, at most, hardware precision (approximately 15 base-10 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.
See the section Execution-Time Feedback for information on how to determine if such copy operations were required.
The incorporation of the NAG library into Maple is an ongoing project. At least one general-purpose 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 Execution-Time Feedback for information on how to determine if a NAG routine is used to complete a LinearAlgebra computation.
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.
In general, LinearAlgebra functions use NAG routines to perform computations on input objects with datatype=sfloat or datatype=float 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 user-defined 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.
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 floating-point 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.
with( LinearAlgebra ):
A := <<2.2,4.0,-6.1>|<0,1.8,3>|<-3.2,-5,7.4>>;
b := <2.1,5.9,3>;
LinearSolve( A, b );
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' );
x := LinearSolve( [p, lu], b );
x2 := LinearSolve( [p, lu], <3.1,1.3,-.27> );
By preserving the packed form of the factorization, unnecessary copy operations and repeated factorizations are not performed in this sequence of computations.
Details of the LinearAlgebra Package
Download Help Document