Skip to content

Commit

Permalink
Add new mex to run inner loop of spinwave.m in parallel (#163)
Browse files Browse the repository at this point in the history 
* Add new mex to run inner loop of spinwave.m in parallel

The new mex file uses Eigen for linear algebra and standard
C++11 threads. Eigen operations are numerically slightly
different from the BLAS/LAPACK/MKL routines used by Matlab
and other mex files so there will be discrepancies.

Also, Eigen eig/chol operations are serial unlike MKL
so will be much slower for large system. So, the new code
has a check and will use the old mex files for these systems

* Fix sw_mex.m issue downloading Eigen

* Compute incommensurate phase in mex

* swloop: improve error handling; add unit tests

* Add field and biquadratic mex systemtests

* Add twin rotation to swloop mex

* PR changes

* Add comments to swloop input fields. Updated docstring
* Removed erroneous "Very baaaad!" error message
* Change nThreads parameter in C++ to default to number of cores
  if input Matlab parameter (pref.nthread) is set to -1
* Add pref.nspinlarge parameter to determine whether to use
  new swloop (good for small systems) or old mex files.

* Add unit test for swpref.nthread and nspinlarge
  • Loading branch information
@mducle
mducle authored Apr 9, 2024
1 parent ce8dd43 commit f11cf75
Show file tree
Hide file tree
Showing 12 changed files with 994 additions and 31 deletions.
8 changes: 7 additions & 1 deletion +sw_tests/+system_tests/systemtest_spinwave.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,9 @@ function get_reference_data(testCase)
testCase.reference_data = load(fname);
end
end
function disable_mex_setup(testCase)
swpref.setpref('usemex', false);
end
end

methods (TestClassTeardown)
Expand All @@ -39,6 +42,9 @@ function save_reference_data(testCase)
save(fname, '-struct', 'ref_dat');
end
end
function disable_mex_teardown(testCase)
swpref.setpref('usemex', false);
end
end

methods (Static)
Expand Down Expand Up @@ -110,7 +116,7 @@ function save_reference_data(testCase)
function fieldname = get_fieldname(testCase, pars)
if isempty(pars)
fieldname = 'data';
elseif ischar(pars);
elseif ischar(pars)
fieldname = pars;
else
fieldname = ['d' reshape(dec2hex(testCase.get_hash(pars)),1,[])];
Expand Down
9 changes: 7 additions & 2 deletions +sw_tests/+system_tests/systemtest_spinwave_biquadratic.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,10 @@
reference_data_file = 'systemstest_spinwave_biquadratic.mat';
end

properties (TestParameter)
mex = {0, 1};
end

methods (TestMethodSetup)
function prepareForRun(testCase)
% From Tutorial 28, to test the biquadratic interactions functionality, theory from PHYSICAL REVIEW B 85, 054409 (2012)
Expand All @@ -28,12 +32,13 @@ function prepareForRun(testCase)
end

methods (Test)
function test_biquadratic(testCase)
function test_biquadratic(testCase, mex)
swpref.setpref('usemex', mex);
fcc = testCase.swobj;
spec = fcc.spinwave({[1 0 0] [0 0 0] [1/2 1/2 0] [1/2 1/2 1/2] [0 0 0] 50});
spec = sw_egrid(spec);
spec = sw_omegasum(spec,'zeroint',1e-5,'emptyval',0,'tol',1e-4);
testCase.generate_or_verify(spec, {}, struct(), 'approxSab', 0.01);
testCase.generate_or_verify(spec, {}, struct(), 'approxSab', 0.02);
end
end

Expand Down
5 changes: 5 additions & 0 deletions +sw_tests/+system_tests/systemtest_spinwave_yb2ti2o7.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,11 @@ function test_yto_twin(testCase)
testCase.swobj.addtwin('axis', [1 -1 0], 'phid', 90);
testCase.test_yto(4, {[-0.5 -0.5 -0.5] [2 2 2]});
end
function test_yto_mex(testCase, B)
% Tests mex with just first Q setting
swpref.setpref('usemex', 1)
testCase.test_yto(B, {[-0.5 -0.5 -0.5] [2 2 2]});
end
end

end
53 changes: 51 additions & 2 deletions +sw_tests/+unit_tests/unittest_spinw_spinwave.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
% Test directions and literal qpts work
qpts_h5 = {{[0 0 0], [1 0 0], 5}, ...
[0:0.25:1; zeros(2,5)]};
mex = {0, 1};
mex = {0, 'old', 1};
end

methods (TestClassSetup)
Expand Down Expand Up @@ -156,13 +156,15 @@ function test_sw_qh5_tid(testCase)
end
methods (Test)
function test_sw_qh5(testCase, qpts_h5, mex)
testCase.assumeNotEqual(mex, 1); % swloop outputs c.c. Sab so fails here
swpref.setpref('usemex', mex);
testCase.disable_warnings('spinw:spinwave:NonPosDefHamiltonian');
sw_out = testCase.swobj.spinwave(qpts_h5);
expected_sw = testCase.get_expected_sw_qh5();
testCase.verify_spinwave(sw_out, expected_sw);
end
function test_sw_qh5_sortmode(testCase, mex)
testCase.assumeNotEqual(mex, 1); % swloop outputs c.c. Sab so fails here
swpref.setpref('usemex', mex);
testCase.disable_warnings('spinw:spinwave:NonPosDefHamiltonian');
sw_out = testCase.swobj.spinwave(testCase.qh5, 'sortMode', false);
Expand All @@ -188,6 +190,7 @@ function test_sw_qh5_nformula(testCase, mex)
testCase.verify_spinwave(sw_out_nformula, expected_sw);
end
function test_sw_qh5_periodic(testCase, mex)
testCase.assumeNotEqual(mex, 1); % swloop outputs c.c. Sab so fails here
swpref.setpref('usemex', mex);
% Test qpts in different BZ give same omega, Sab
qpts = testCase.qh5 + 1;
Expand All @@ -199,6 +202,7 @@ function test_sw_qh5_periodic(testCase, mex)
testCase.verify_spinwave(sw_out, expected_sw);
end
function test_sw_qh5_perpendicular(testCase, mex)
testCase.assumeNotEqual(mex, 1); % swloop outputs c.c. Sab so fails here
swpref.setpref('usemex', mex);
% Test qpts in perpendicular direction give flat modes
qpts = [zeros(1, 5); 0:0.25:1; 0:0.25:1];
Expand Down Expand Up @@ -228,6 +232,7 @@ function test_sw_qh5_saveH_saveV(testCase, mex)
testCase.verify_spinwave(sw_out, expected_sw);
end
function test_sw_qh5_title(testCase, mex)
testCase.assumeNotEqual(mex, 1); % swloop outputs c.c. Sab so fails here
swpref.setpref('usemex', mex);
title = 'Example title';
testCase.disable_warnings('spinw:spinwave:NonPosDefHamiltonian');
Expand Down Expand Up @@ -392,6 +397,7 @@ function test_twin(testCase, mex)
testCase.verify_spinwave(sw_out, expected_sw);
end
function test_notwin(testCase, mex)
testCase.assumeNotEqual(mex, 1); % swloop outputs c.c. Sab so fails here
swpref.setpref('usemex', mex);
% Create copy to avoid changing obj for other tests
swobj_twin = copy(testCase.swobj);
Expand All @@ -417,8 +423,10 @@ function test_cmplxBase_equivalent_with_tri(testCase, mex)
end
function test_cmplxBase_fails_with_chain(testCase, mex)
swpref.setpref('usemex', mex);
if mex
if ischar(mex)
err ='chol_omp:notposdef';
elseif mex == 1
err ='swloop:notconverge';
else
err= 'spinw:spinwave:NonPosDefHamiltonian';
end
Expand Down Expand Up @@ -515,6 +523,7 @@ function test_hermit(testCase, mex)
testCase.assertGreaterThan(sum(abs(imag(spec.omega(:)))), 0);
end
function test_sw_qh5_tol(testCase, mex)
testCase.assumeNotEqual(mex, 1); % swloop outputs c.c. Sab so fails here
swpref.setpref('usemex', mex);
tol = 5e-4;
qpts = testCase.qh5;
Expand All @@ -537,6 +546,7 @@ function test_sw_qh5_tol(testCase, mex)
testCase.verify_spinwave(sw_out, expected_sw);
end
function test_sw_qh5_omega_tol(testCase, mex)
testCase.assumeNotEqual(mex, 1); % swloop outputs c.c. Sab so fails here
swpref.setpref('usemex', mex);
if mex
err ='chol_omp:notposdef';
Expand Down Expand Up @@ -575,6 +585,45 @@ function test_no_magstr_causes_error(testCase)
@() swobj.spinwave(testCase.qh5), ...
'spinw:spinwave:NoMagneticStr');
end
function test_mex_prefs_nthreads(testCase)
% Tests that when nthreads set to -1 mex will parallelise over all cores
hkl = {[0 0 0] [1 0 0] [0 1 0] 2000};
nCores = maxNumCompThreads();
if nCores > 2
% If only 2 cores, times may not be different enough to notice
mexpref = swpref.getpref('usemex');
nthrprf = swpref.getpref('nthread');
swpref.setpref('usemex', 1);
swpref.setpref('nthread', -1);
swobj = copy(testCase.swobj_tri);
spec1 = swobj.spinwave({[0 0 0] [1 1 1] 50}); % Run once for the JIT compiler
t0 = tic; spec1 = swobj.spinwave(hkl); t1 = toc(t0);
swpref.setpref('nthread', 1);
t0 = tic; spec1 = swobj.spinwave(hkl); t2 = toc(t0);
testCase.verifyTrue(t1 < t2);
swpref.setpref('usemex', mexpref.val);
swpref.setpref('nthread', nthrprf.val);
end
end
function test_mex_prefs_nspinlarge(testCase)
swobj = copy(testCase.swobj_tri);
swobj.addmatrix('label', 'K', 'value', 0.1);
swobj.addaniso('K');
mexpref = swpref.getpref('usemex');
nslpref = swpref.getpref('nspinlarge');
swpref.setpref('usemex', 1);
omega = [1e-5; -1e-5];
Sab = reshape([[0.5 0 -0.5j; 0 0 0; 0.5j 0 0.5] [0.5 0 0.5j; 0 0 0; -0.5j 0 0.5]], 3, 3, 2);
swloop_mock = sw_tests.utilities.mock_function('swloop', {omega, Sab, 1, 0});
spec1 = swobj.spinwave([1; 1; 1], 'hermit', false);
testCase.assertEqual(swloop_mock.n_calls, 1);
swobj.genmagstr('nExt', 0.01); % Convert to commensurate with 9 spins in unit cell (3x3)
swpref.setpref('nspinlarge', 5);
spec1 = swobj.spinwave([1; 1; 1], 'hermit', false);
testCase.assertEqual(swloop_mock.n_calls, 1); % Make sure swloop wasn't called this time
swpref.setpref('usemex', mexpref.val);
swpref.setpref('nspinlarge', nslpref.val);
end
end
methods (Test, TestTags = {'Symbolic'})
function test_sw_symbolic_no_qpts(testCase)
Expand Down
3 changes: 3 additions & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,3 +35,6 @@ includedSupportPackages.txt
mccExcludedFiles.log
requiredMCRProducts.txt
unresolvedSymbols.txt

coverage*xml
junit*xml
2 changes: 1 addition & 1 deletion external/eigorth.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@
nStack = size(M,3);

% Use OpenMP parallelised mex file if it exists
if nStack>1 && useMex
if nStack>1 && any(useMex)
% eigenvalues are already orthogonalised by eig_omp
[V, D] = eig_omp(M,'orth','sort','descend');
return
Expand Down
Loading

0 comments on commit f11cf75

Please sign in to comment.