Performance Improvements in Maple 2023

Rational Linear Solver Speedup


•

The underlying engine that solve uses to solve systems of linear equations with rational coefficients has been improved for many types of systems. The recently added SolveTools:LinearSolvers:RationalDense is now used for most systems of equations with more than 5 percent of their coefficients nonzero that are also overdetermined or highly underdetermined with user specified free variables.

•

Here is an example of a highly overdetermined system with many solutions where the new method performs much better.

>

N := 200: n := 100: r := rand(1000..1000):

>

M := Matrix(N+n, N, (i,j)>ifelse(j<(in), 0, r())):

>

V1 := M . Vector([seq(r(),i=1..N)]):V2 := M . Vector([seq(r(),i=1..N)]):

>

vars := [seq(cat(`x__`,i),i=1..N+1)]:

>

sys := (lhsrhs)~(LinearAlgebra:GenerateEquations( <V1  V2  M>, vars)):

>

CodeTools:Usage(SolveTools:LinearSolvers:Rational(sys, indets(sys),dense=false) ):

memory used=64.06MiB, alloc change=66.62MiB, cpu time=733.00ms, real time=719.00ms, gc time=75.30ms
 
•

The new default is significantly faster on this system

>

CodeTools:Usage(SolveTools:LinearSolvers:Rational(sys, indets(sys)) ):

memory used=27.11MiB, alloc change=2.94MiB, cpu time=185.00ms, real time=168.00ms, gc time=36.64ms
 


Hardware FloatingPoint Evaluation of Numeric Matrices


•

When initializing hardwarefloat matrix, vector, and array data structures the kernel will opt for computing with numeric implementations of the initializer when possible. This results in dramatic speedups in some cases. The following example used to take several seconds and now returns instantly:

>

signal:=Vector(10^5,i>sin(i),datatype=float[8]):

•

A simple initialization of a float[8] rtable with integer values is now about 3.5 times faster.

>

v:=Vector(10^6, i>i, datatype=float[8]):

•

Evaluation of sums, products and functions now use less memory and a more direct path to hardwarefloat evaluation. The following example is more than two times faster in Maple 2023.

>

X := RandomVector(N, generator = 0 .. 1.0, datatype = float[8]):

>

Y := RandomVector(N, generator = 0 .. 1.0, datatype = float[8]):

>

total := Vector(N, i > ifelse(X[i]^2 + Y[i]^2 < 1, 1, 0), datatype=float[8]):

•

Type checking of datatype=float[8] rtables has also improved. The following example is hundreds of times faster in Maple 2023.

>

A := [seq(Array(1..1,[i],datatype=float[8]),i=1..10^4)]:

>

map(type,A,'rtable'(datatype=float[8])):

•

When a hardwarebased rtable's elements are extracted they remain as a doubleprecision HFloat datastructure. Support for direct operation on these structures has been extended to include integeroutput commands such as floor, round, and ceil. This gives a small improvement that adds up over millions of computations. Timings in the following example are now almost half of what they were in the previous version.

>

v := LinearAlgebra:RandomVector(n,datatype=float[8]):
tt := time():
for i from 1 to numelems(v) do
v[i] := floor(v[i]/10);
end do:



The Units Package


•

The Units package, and in particular its Units:Simple subpackage have received several upgrades that make them faster, sometimes several orders of magnitude faster.

•

The Units:Simple package works by calling the command Units:TestDimensions every time an arithmetic operation is performed. This command has been reimplemented to be much faster. In usersubmitted worksheets that heavily use units, we often see speedups between 3 times and 100 times for the whole worksheet.

•

The following example shows a simple forward Euler discretization of constant acceleration for 2 seconds in 2000 steps.

>

make_sys := proc(x :: symbol, v :: symbol, dt :: symbol, a :: symbol, n :: nonnegint, dtval :: algebraic, aval :: algebraic, $)
local i;
return [x[0] = 0, v[0] = 0, a = aval, dt = dtval,
seq(v[i+1] = v[i] + dt * a, i = 0 .. n1),
seq(x[i+1] = x[i] + dt * v[i], i = 0 .. n1)];
end proc;

${{\mathrm{make\_sys}}{:=}{\mathbf{proc}}\left({x}{::}{\mathrm{symbol}}{\,}{v}{::}{\mathrm{symbol}}{\,}{\mathrm{dt}}{::}{\mathrm{symbol}}{\,}{a}{::}{\mathrm{symbol}}{\,}{n}{::}{\mathrm{nonnegint}}{\,}{\mathrm{dtval}}{::}{\mathrm{algebraic}}{\,}{\mathrm{aval}}{::}{\mathrm{algebraic}}{\,}{\$}\right)\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{\mathbf{local}}\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{i}{\;}\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{\mathbf{return}}\phantom{\rule[0.0ex]{0.5em}{0.0ex}}\left[{x}{\[}{0}{\]}{\=}{0}{\,}{v}{\[}{0}{\]}{\=}{0}{\,}{a}{\=}{\mathrm{aval}}{\,}{\mathrm{dt}}{\=}{\mathrm{dtval}}{\,}{\mathrm{seq}}{}\left({v}{\[}{i}{\+}{1}{\]}{\=}{v}{\[}{i}{\]}{\+}{\mathrm{dt}}{\*}{a}{\,}{i}{\=}{0}{..}{n}{}{1}\right){\,}{\mathrm{seq}}{}\left({x}{\[}{i}{\+}{1}{\]}{\=}{x}{\[}{i}{\]}{\+}{\mathrm{dt}}{\*}{v}{\[}{i}{\]}{\,}{i}{\=}{0}{..}{n}{}{1}\right)\right]\phantom{\rule[0.0ex]{0.5em}{0.0ex}}{\mathbf{end\; proc}}}}$
 (6) 
>

sys := make_sys(x, v, dt, a, 2000, 0.001*Unit(s), 9.8*Unit(m/s^2)):

•

The TestDimensions command without further options verifies that this is a dimensionally consistent set of equations.

>

Units:TestDimensions(sys);

•

With the option below, it computes the dimensions of ${x}_{0}$ and ${v}_{2000}$.

>

Units:TestDimensions(sys, output=units(x[0], v[2000]));

$\left\{{{v}}_{{2000}}{::}\left(\u27e6\frac{{m}}{{s}}\u27e7\right){\,}{{x}}_{{0}}{::}\u27e6{m}\u27e7\right\}$
 (9) 
•

This example is about 200 times faster in Maple 2023 than before.

•

The Units:Simple package has been optimized to improve performance of base operations. In order to provide automatic units simplifications and verify consistency, the Units:Simple package overloads core operations such as add, multiply and exponentiations. In previous versions this could add significant overhead even when calculating with expressions that do not contain units. In Maple 2023 this overhead has been greatly reduced.

•

For example, the following example computes a numerical approximation to pi using a Monte Carlo method.

>

tt := time[real]():
with(Units[Simple]):
N := 10^6:
with(LinearAlgebra):
X := RandomVector(N, generator = 0 .. 1.0, datatype = float[8]):
Y := RandomVector(N, generator = 0 .. 1.0, datatype = float[8]):
total := Vector(N, i > ifelse(X[i]^2 + Y[i]^2 < 1, 1, 0), datatype=float[8]):
add(total[i], i = 1..N) * 4 / N;
time[real]()tt:

${3.14037200000000}$
 (10) 
•

On a representative computer, this example now completes in 7 seconds, compared to more than 60 seconds in prior versions.



Performance Improvements for Plots


•

In Maple 2022, we introduced a new adaptive plotting engine. This engine has been sped up significantly in Maple 2023. For example, the following simple plot has been sped up by about a factor of 4, and it uses about 4 times less memory:

•

This example runs about 15 times faster in Maple 2023, and it uses about 10 times less memory.

>

plot(sin(ln(exp(x) + 1)), x=5..5);



Faster Water and Steam Properties


•

The ThermophysicalData:CoolProp subpackage has been updated to CoolProp 6.4.1. Most notably, the update introduces IAPWSIF97 water and steam propertiesthese are faster to compute than those given by IAPWS95 (which is used by the default Helmholtz equation of state). This is at the expense of slightly less accuracy and applicability over a smaller envelope of temperatures and pressures.

•

IAPWSIF97 is recommended when water properties need to be computed quickly and continuously, such as in the simulation and optimization of steam power cycles.

•

Specific heat capacity of water using IAPWIF97

>

ThermophysicalData:CoolProp:Property(C, "IF97::Water", temperature=500*Unit('K'), pressure=1*Unit('atm'));

${1981.542297}{}\u27e6\frac{{J}}{{\mathrm{kg}}{}{K}}\u27e7$
 (11) 
•

Specific heat capacity of water using the Helmholtz equation of state (which can be specified by using "HEOS::Water" or "Water" as the fluid).

>

ThermophysicalData:CoolProp:Property(C, "HEOS::Water", temperature=500*Unit('K'), pressure=Unit('atm'));

${1981.609716}{}\u27e6\frac{{J}}{{\mathrm{kg}}{}{K}}\u27e7$
 (12) 
•

Evaluate the specific heat capacity of water N times using IAPWIF97 and the Helmholtz equation of state

>

tt:=time[real]():
for i to N do
ThermophysicalData:CoolProp:Property(C, "IF97::Water", pressure=1*Units:Unit(atm), enthalpy=Unit(('kJ'/'kg'))):
end do:
time[real]()tt;

>

tt:=time[real]():
for i to N do
ThermophysicalData:CoolProp:Property(C, "HEOS::Water", pressure=1*Units:Unit(atm), enthalpy=Unit(('kJ'/'kg'))):
end do:
time[real]()tt;



evalhf Hardware FloatingPoint Subsystem


•

A call to evalhf evaluates an expression to a numerical value using the floatingpoint hardware of the underlying system. The evaluation is done in double precision. The argument passed to evalhf must be an expression that evaluates to float[8] or complex[8] values as singletons or in an array, matrix, or vector. In Maple 2023, the evalhf subsystem has been extended to handle all builtin (kernel) procedures that process arbitrary inputs and return hardware values.

•

Consider the following example which uses the builtin function, rhs. Evaluation of rhs and its arguments will be performed outside of evalhf so the nonhardware range structure returned by rtable_dims can be computed, and subsequently given to rhs, which will return a hardware compatible structure back to the evalhf subsystem.

>

A := Array(1..2,datatype=float[8]):

>

evalhf( rhs(rtable_dims(A,1) ) );

•

In this example, the given loopsum procedure uses the builtin, numelem, which was formerly not supported by evalhf.

>

loopsum := proc( A )
local s := 0;
local n := numelems(A);
for local i from 1 to n do
s := s + A[i];
end do;
s;
end proc:

>

v := LinearAlgebra:RandomVector(100,datatype=float[8]):

•

These additions extend the reach of evalhf allowing a wider variety of code to be used inside evalhf.



Importing Data from a File


•

The ImportVector command can be used to load data from a variety of different file formats. Performance of importing commaseparated value (CSV) files has been improved on Windows. The following example completes in about half the time it used to.

>

data:=ImportVector("C:\\data.csv", source=csv[standard],datatype=float[8]):



BitShift Operations


•

Two new commands integerdivq2exp and integermul2exp provide what are considered hardware bitshift operations. They are fast methods for dividing or multiplying an integer by a power of 2. integermul2exp computes i*2^pow, and integerdivq2exp calculates trunc(i / 2^pow).

>

integermul2exp( 10, 4 )  10*2^4;

>

integerdivq2exp(255,4)  trunc(255/2^4);


