Testsuite: improve and integrate unit test for A·(B⊗C) MEX files

In particular, add a test for the sparse A·(B⊗C) where B≠C.

tests/.gitignore

@@ -95,7 +95,6 @@ wsOct
 !/kalman_filter_smoother/testsmoother.m
 !/kalman_steady_state/test1.m
 !/kronecker/nash_matrices.mat
-!/kronecker/small_matrices.mat
 !/kronecker/test_kron.m
 !/load_octave_packages.m
 !/ls2003/data_ca1.m
@@ -152,6 +151,7 @@ wsOct
 !/run_all_unitary_tests.m
 !/run_block_byte_tests_matlab.m
 !/run_block_byte_tests_octave.m
+!/run_kronecker_tests.m
 !/run_m_script.m
 !/run_o_script.m
 !/run_reporting_test_matlab.m

tests/Makefile.am

@@ -1075,7 +1075,8 @@ M_TRS_FILES += 	run_block_byte_tests_matlab.m.trs \
 		run_all_unitary_tests.m.trs \
 		histval_initval_file_unit_tests.m.trs \
 		test_aggregate_routine_1_2.m.trs \
-		test_aggregate_routine_1_2_3.m.trs
+		test_aggregate_routine_1_2_3.m.trs \
+		run_kronecker_tests.m.trs
 
 M_XFAIL_TRS_FILES = $(patsubst %.mod, %.m.trs, $(XFAIL_MODFILES))
 
@@ -1085,7 +1086,8 @@ O_TRS_FILES += $(patsubst %.mod, %.o.trs, $(PARTICLEFILES))
 O_TRS_FILES += 	run_block_byte_tests_octave.o.trs \
 		run_reporting_test_octave.o.trs \
 		run_all_unitary_tests.o.trs \
-	 	histval_initval_file_unit_tests.o.trs
+		histval_initval_file_unit_tests.o.trs \
+		run_kronecker_tests.o.trs
 
 O_XFAIL_TRS_FILES = $(patsubst %.mod, %.o.trs, $(XFAIL_MODFILES))
 
@@ -1255,7 +1257,10 @@ EXTRA_DIST = \
 	histval_initval_file/my_assert.m \
 	histval_initval_file/ramst_data.xls \
 	histval_initval_file/ramst_data.xlsx \
-	histval_initval_file/ramst_initval_file_data.m
+	histval_initval_file/ramst_initval_file_data.m \
+	run_kronecker_tests.m \
+	kronecker/test_kron.m \
+	kronecker/nash_matrices.mat
 
 if ENABLE_MATLAB
 check-local: check-matlab

tests/kronecker/test_kron.m

 function info = test_kron(test,number_of_threads)
-% Copyright (C) 2007-2020 Dynare Team
+% Copyright (C) 2007-2021 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -28,10 +28,11 @@ if test == 1
     disp('')
     disp('I''m building the test problem...')
     tic
-    percentage_of_non_zero_elements = 10e-4;
-    NumberOfVariables = 549;%100;
-    NumberOfEquations = 256;%50
-    NumberOfColsInB   = 50 ; 
+    percentage_of_non_zero_elements = 1e-3;
+    NumberOfVariables = 120;
+    NumberOfEquations = 70;
+    NumberOfColsInB   = 50;
+    NumberOfColsInC   = 60;
     A = zeros(NumberOfEquations,NumberOfVariables^2);
     for i = 1:NumberOfEquations
         for j = 1:NumberOfVariables
@@ -47,38 +48,50 @@ if test == 1
     end
     A = sparse(A);
     B = randn(NumberOfVariables,NumberOfColsInB);
+    C = randn(NumberOfVariables,NumberOfColsInC);
     disp('Done!')
     toc
     disp('')
-    disp('Computation of A*kron(B,B) with the mex file (v1):')
+    disp('Computation of A*kron(B,B) with the mex file:')
     tic 
     D1 = sparse_hessian_times_B_kronecker_C(A,B,number_of_threads);
     toc
-    disp('')
-    disp('Computation of A*kron(B,B) with the mex file (v2):')
+    disp('Computation of A*kron(B,B) with two nested loops:')
     tic
-    D2 = sparse_hessian_times_B_kronecker_C(A,B,B,number_of_threads);
+    D2 = zeros(NumberOfEquations,NumberOfColsInB*NumberOfColsInB);
+    k = 1;
+    for i1 = 1:NumberOfColsInB
+        for i2 = 1:NumberOfColsInB
+            D2(:,k) = A*kron(B(:,i1),B(:,i2));
+            k = k+1;
+        end
+    end
     toc
     disp('');
     disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
-    if max(max(abs(D1-D2)))>1e-10
+    if max(max(abs(D1-D2)))>1e-13
         test1_1=0; 
     end
-    disp(' ')
-    disp('Computation of A*kron(B,B) with two nested loops:')
+    disp('')
+    disp('Computation of A*kron(B,C) with the mex file:')
     tic
-    D3 = zeros(NumberOfEquations,NumberOfColsInB*NumberOfColsInB);
+    D3 = sparse_hessian_times_B_kronecker_C(A,B,C,number_of_threads);
+    toc
+    disp('Computation of A*kron(B,C) with two nested loops:')
+    tic
+    D4 = zeros(NumberOfEquations,NumberOfColsInB*NumberOfColsInC);
     k = 1;
     for i1 = 1:NumberOfColsInB
-        for i2 = 1:NumberOfColsInB
-            D3(:,k) = A*kron(B(:,i1),B(:,i2)); 
+        for i2 = 1:NumberOfColsInC
+            D4(:,k) = A*kron(B(:,i1),C(:,i2));
             k = k+1;
         end
     end
     toc
+    disp(' ')
     disp('');
-    disp(['Difference between D1 and D3 = ' num2str(max(max(abs(D1-D3))))]);
-    if max(max(abs(D1-D3)))>1e-10
+    disp(['Difference between D3 and D4 = ' num2str(max(max(abs(D3-D4))))]);
+    if max(max(abs(D3-D4)))>1e-13
         test1_2=0; 
     end
 % $$$ FOR THE DIMENSIONS CONSIDERED HERE THIS PART WILL RESULT IN A OUT OF MEMORY ERROR.   
@@ -163,7 +176,7 @@ if test==3
     disp('Test with full format matrix -- 1(a)')
     D1 = A*kron(B,C);
     tic
-    D2 = A_times_B_kronecker_C(A,B,C,number_of_threads);
+    D2 = A_times_B_kronecker_C(A,B,C);
     toc
     disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
     if max(max(abs(D1-D2)))>1e-10
@@ -172,7 +185,7 @@ if test==3
     disp('Test with full format matrix -- 1(b)')
     D1 = A*kron(B,B);
     tic
-    D2 = A_times_B_kronecker_C(A,B,number_of_threads);
+    D2 = A_times_B_kronecker_C(A,B);
     toc
     disp(['Difference between D1 and D2 = ' num2str(max(max(abs(D1-D2))))]);
     if max(max(abs(D1-D2)))>1e-10

tests/run_kronecker_tests.m

+% Copyright (C) 2021 Dynare Team
+%
+% This file is part of Dynare.
+%
+% Dynare is free software: you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation, either version 3 of the License, or
+% (at your option) any later version.
+%
+% Dynare is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
+
+top_test_dir = getenv('TOP_TEST_DIR');
+addpath([top_test_dir filesep '..' filesep 'matlab']);
+dynare_config
+
+cd kronecker
+
+disp('**** Testing sparse_hessian_times_B_kronecker_C MEX...')
+disp('')
+disp('** Test 1')
+info(1) = test_kron(1, num_procs);
+disp('')
+disp('** Test 2')
+info(2) = test_kron(2, num_procs);
+disp('')
+disp('**** Testing A_times_B_kronecker_C MEX...')
+info(3) = test_kron(3, num_procs);
+
+cd ..
+
+if isoctave
+    ext = '.o.trs';
+else
+    ext = '.m.trs';
+end
+fid = fopen([ 'run_kronecker_tests' ext ], 'w+');
+num_failed_tests = sum(~info);
+tests = { 'sparse1', 'sparse2', 'dense' };
+failed_tests = tests(find(~info));
+if num_failed_tests > 0
+  fprintf(fid,':test-result: FAIL\n');
+  fprintf(fid,':number-tests: 3\n');
+  fprintf(fid,':number-failed-tests: %d\n', num_failed_tests);
+  fprintf(fid,':list-of-failed-tests: %s\n', failed_tests{:});
+else
+  fprintf(fid,':test-result: PASS\n');
+  fprintf(fid,':number-tests: 3\n');
+  fprintf(fid,':number-failed-tests: 0\n');
+end
+fclose(fid);