Skip to content

Commit

Permalink
240312.030622.HKT support half-precision in MATLAB
Browse files Browse the repository at this point in the history
  • Loading branch information
zaikunzhang committed Mar 11, 2024
1 parent 272f637 commit 3941491
Show file tree
Hide file tree
Showing 7 changed files with 76 additions and 31 deletions.
13 changes: 10 additions & 3 deletions matlab/interfaces/private/getmax.m
Original file line number Diff line number Diff line change
Expand Up @@ -25,12 +25,19 @@
% maxfloat, minfloat, and maxpow10 are intended to the return of the Fortran intrinsics HUGE, TINY,
% and RANGE corresponding to the given precision. They are the largest finite value, the smallest
% positive normalized value, and the decimal exponent range of the floating point numbers, respectively.
if strcmpi(precision, 'single')
switch lower(precision)
case 'half' % IEEE 754 half precision (float16); followed by nagfor 7.1+ when providing REAL16.
maxfloat = 2^16 - 2^5; %65504;
minfloat = 2^(-14); %6.1035e-05;
case 'single'
maxfloat = realmax('single');
minfloat = realmin('single');
else % Even if precision = 'quadruple'
case {'double', 'quadruple'}
maxfloat = realmax('double');
minfloat = realmin('double');
otherwise
% Private/unexpected error
error(sprintf('%s:InvalidInput', funname), '%s: UNEXPECTED ERROR: invalid precision received.', funname);
end
% According to Sec. 16.9.170 of Fortran 2023 Interpretation Document J3/24-007, the Fortran intrinsic
% RANGE returns the following value.
Expand All @@ -44,7 +51,7 @@
case {'bound'}
maxnum = 0.25 * maxfloat;
case {'fun', 'func', 'function', 'con', 'constr', 'constraint'}
maxnum = 10^min(30, floor(maxpow10 / 2));
maxnum = 10^max(4, min(30, floor(maxpow10 / 2)));
otherwise
% Private/unexpected error
error(sprintf('%s:InvalidInput', funname), '%s: UNEXPECTED ERROR: invalid data_type received.', funname);
Expand Down
4 changes: 2 additions & 2 deletions matlab/interfaces/private/preprima.m
Original file line number Diff line number Diff line change
Expand Up @@ -117,8 +117,8 @@
probinfo.raw_data = struct('objective', fun, 'x0', x0, 'Aineq', Aineq, 'bineq', bineq, ...
'Aeq', Aeq, 'beq', beq, 'lb', lb, 'ub', ub, 'nonlcon', nonlcon, 'options', options);

% Decide the precision ('single', 'double', or 'quadruple') of the real calculation within the
% Fortran solvers. This is needed ONLY by `boundmax`, `funcmax`, and `constrmax` defined below by
% Decide the precision ('half', 'single', 'double', or 'quadruple') of the real calculation within
% the Fortran solvers. This is needed ONLY by `boundmax`, `funcmax`, and `constrmax` defined below by
% calling `getmax`. These three numbers will be used in `pre_x0`, `pre_fun`, and `pre_nonlcon`
% respectively in the sequel. Note the following.
% 1. `precision` takes effect only if Fortran solvers are called (i.e., when options.fortran = true).
Expand Down
2 changes: 1 addition & 1 deletion matlab/setup_tools/all_precisions_possible.m
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
% the Fortran solvers in this package; `default_precision`, if requested, is the name of the default
% precision.

precision_list = {'double', 'single', 'quadruple'};
precision_list = {'half', 'single', 'double', 'quadruple'};

default_precision = 'double';

Expand Down
12 changes: 12 additions & 0 deletions matlab/setup_tools/compile.m
Original file line number Diff line number Diff line change
Expand Up @@ -241,17 +241,29 @@ function prepare_header(header_file, precision, debug_flag)
%PREPARE_HEADER prepares `header_file` for the compilation according to `precision` and `debug_flag`.

switch precision
case {'h', 'half'}
rep_str(header_file, '#define PRIMA_REAL_PRECISION 32', '#define PRIMA_REAL_PRECISION 16');
rep_str(header_file, '#define PRIMA_REAL_PRECISION 64', '#define PRIMA_REAL_PRECISION 16');
rep_str(header_file, '#define PRIMA_REAL_PRECISION 128', '#define PRIMA_REAL_PRECISION 16');
rep_str(header_file, '#define PRIMA_HP_AVAILABLE 0', '#define PRIMA_HP_AVAILABLE 1');
rep_str(header_file, '#define PRIMA_QP_AVAILABLE 1', '#define PRIMA_QP_AVAILABLE 0');
case {'s', 'single'}
rep_str(header_file, '#define PRIMA_REAL_PRECISION 16', '#define PRIMA_REAL_PRECISION 32');
rep_str(header_file, '#define PRIMA_REAL_PRECISION 64', '#define PRIMA_REAL_PRECISION 32');
rep_str(header_file, '#define PRIMA_REAL_PRECISION 128', '#define PRIMA_REAL_PRECISION 32');
rep_str(header_file, '#define PRIMA_HP_AVAILABLE 1', '#define PRIMA_HP_AVAILABLE 0');
rep_str(header_file, '#define PRIMA_QP_AVAILABLE 1', '#define PRIMA_QP_AVAILABLE 0');
case {'q', 'quadruple'}
rep_str(header_file, '#define PRIMA_REAL_PRECISION 16', '#define PRIMA_REAL_PRECISION 128');
rep_str(header_file, '#define PRIMA_REAL_PRECISION 32', '#define PRIMA_REAL_PRECISION 128');
rep_str(header_file, '#define PRIMA_REAL_PRECISION 64', '#define PRIMA_REAL_PRECISION 128');
rep_str(header_file, '#define PRIMA_HP_AVAILABLE 1', '#define PRIMA_HP_AVAILABLE 0');
rep_str(header_file, '#define PRIMA_QP_AVAILABLE 0', '#define PRIMA_QP_AVAILABLE 1');
otherwise
rep_str(header_file, '#define PRIMA_REAL_PRECISION 16', '#define PRIMA_REAL_PRECISION 64');
rep_str(header_file, '#define PRIMA_REAL_PRECISION 32', '#define PRIMA_REAL_PRECISION 64');
rep_str(header_file, '#define PRIMA_REAL_PRECISION 128', '#define PRIMA_REAL_PRECISION 64');
rep_str(header_file, '#define PRIMA_HP_AVAILABLE 1', '#define PRIMA_HP_AVAILABLE 0');
rep_str(header_file, '#define PRIMA_QP_AVAILABLE 1', '#define PRIMA_QP_AVAILABLE 0');
end

Expand Down
12 changes: 10 additions & 2 deletions matlab/setup_tools/create_all_precisions.m
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
function create_all_precisions(options_or_directory)
%CREATE_ALL_PRECISIONS creates `all_precisions.m` under the directory containing this script
% according to `options_or_directory`. `all_precisions.m` should return a cell array containing
% the names of all the precisions ('double', 'single', 'quadruple') available for the Fortran
% the names of all the precisions ('half', 'single', 'double', 'quadruple') available for the Fortran
% solvers in this package. It is created in the following way.
%
% 0. We assume that 'double' is always available.
Expand All @@ -14,8 +14,9 @@ function create_all_precisions(options_or_directory)
% 2. If `options_or_directory` is a structure (or empty), then it will be interpreted as compilation
% options. The return of `all_precisions.m` will reflect the precisions available after the compilation.

% Default values for the availability of 'single' and 'quadruple'. They are used only if
% Default values for the availability of 'half', 'single' and 'quadruple'. They are used only if
% `options_or_directory` is a structure (i.e., it is indeed the compilation options).
half_precision = false;
single_precision = true;
quadruple_precision = false;

Expand Down Expand Up @@ -45,12 +46,16 @@ function create_all_precisions(options_or_directory)
elseif ischarstr(options_or_directory)

directory = options_or_directory;
half_precision = isavailable(directory, 'half');
single_precision = isavailable(directory, 'single');
quadruple_precision = isavailable(directory, 'quadruple');

elseif isa(options_or_directory, 'struct')

options = options_or_directory;
if isfield(options, 'half') && islogicalscalar(options.half)
half_precision = options.half;
end
if isfield(options, 'single') && islogicalscalar(options.single)
single_precision = options.single;
end
Expand All @@ -64,6 +69,9 @@ function create_all_precisions(options_or_directory)

% Decide the precision list string.
precision_list_string = '''double''';
if half_precision
precision_list_string = [precision_list_string, ', ''half'''];
end
if single_precision
precision_list_string = [precision_list_string, ', ''single'''];
end
Expand Down
59 changes: 38 additions & 21 deletions matlab/tests/pdv.m
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,12 @@ function pdv()
options.compile = true;
test_dir = prepare_test_dir(fake_solver_name, funname, options);

is_mac_silicon = false;
if ismac
[~, result] = system('uname -v');
is_mac_silicon = contains(result, 'arm64', 'IgnoreCase', true);
end

exception = [];

try
Expand All @@ -35,26 +41,31 @@ function pdv()
clear('setup');
opt=struct();
opt.verbose = true;
opt.single=true;
opt.quadruple=true;
opt.debug=true;
opt.classical=true;
opt.half = is_mac_silicon;
opt.single = true;
opt.double = true;
opt.quadruple = true;
opt.debug = true;
%opt.classical = true;
opt.classical = false;

tic

setup(opt);
testprima
toc

solvers = {'cobyla', 'uobyqa', 'newuoa', 'bobyqa', 'lincoa'};
precisions = {'double', 'single', 'quadruple'};
precisions = {'half', 'single', 'double', 'quadruple'};
precisions = precisions([opt.half, opt.single, opt.double, opt.quadruple]);
debug_flags = {true, false};
variants = {'modern', 'classical'};
%variants = {'modern', 'classical'};
variants = {'modern'};

% Show current path information.
showpath(solvers);

% Test the solvers.
fun = @sin;
x0 = 1;
fun = @chrosen;
x0 = [-1; -1];
for isol = 1 : length(solvers)
solver = str2func(solvers{isol});
solver
Expand All @@ -65,27 +76,33 @@ function pdv()
options.debug = debug_flags{idbg};
for ivar = 1 : length(variants)
options.classical = strcmp(variants{ivar}, 'classical');
if ismac && strcmp(func2str(solver), 'cobyla') && strcmp(options.precision, 'half') && options.classical
% Skip the classical cobyla in half precision on macOS, as it will encounter an infinite cycling.
continue;
end
if ismac && strcmp(func2str(solver), 'bobyqa') && strcmp(options.precision, 'quadruple') && options.classical
% Skip the classical bobyqa in quadruple precision on macOS, as it will encounter a segmentation fault.
continue;
end
for iprint = -4 : 4
options.output_xhist = true;
options.iprint = iprint;
options
format long
[x, f, exitflag, output] = solver(fun, x0, options)
if (ismember(solvers{isol}, matlab_implemented))
options_mat = options;
options_mat.fortran = false;
[x, f, exitflag, output] = solver(fun, x0, options_mat)
end
options.output_xhist = true;
options.maxfun = 100*length(x0);
options.rhoend = 1.0e-3;
options.iprint = randi([-4, 4]);
options
format long
[x, f, exitflag, output] = solver(fun, x0, options)
if (ismember(solvers{isol}, matlab_implemented))
options_mat = options;
options_mat.fortran = false;
[x, f, exitflag, output] = solver(fun, x0, options_mat)
end
end
end
end
end

toc

% Show current path information again at the end of test.
showpath(solvers);

Expand Down
5 changes: 3 additions & 2 deletions setup.m
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@ function setup(varargin)
% setup(solver_name, options) % Compile a solver with `options`
%
% Possible options for compilation:
% - half: whether to compile the half precision of the Fortran solvers (default: false)
% - single: whether to compile the single precision of the Fortran solvers (default: true)
% - quadruple: whether to compile the quadruple precision of the Fortran solvers (default: false)
% - classical: whether to compile the classical variant of the Fortran solvers (default: true)
Expand Down Expand Up @@ -123,7 +124,7 @@ function setup(varargin)
if strcmp(action, 'path')
add_save_path(interfaces, package_name);
% Create `all_precisions.m` and `all_variants.m` under `tools` according to the content of
% `mexdir`. They reflect the precisions ('double', 'single', 'quadruple') and variants
% `mexdir`. They reflect the precisions ('half', 'single', 'double', 'quadruple') and variants
% ('modern', 'classical') of the solvers available under `mexdir`.
create_all_precisions(mexdir);
create_all_variants(mexdir);
Expand All @@ -145,7 +146,7 @@ function setup(varargin)


% Create `all_precisions.m` and `all_variants.m` under `tools` according to `options`.
% They reflect the precisions ('double', 'single', 'quadruple') and variants ('modern','classical')
% They reflect the precisions ('half', 'single', 'double', 'quadruple') and variants ('modern','classical')
% of the solvers available after the compilation.
create_all_precisions(options);
create_all_variants(options);
Expand Down

0 comments on commit 3941491

Please sign in to comment.