Templated GEMM references C when beta == 0.0.

[rileyjmurray]

The fully templated GEMM code (below) references the input value of C(i,j) even when beta == 0.0. As a result, if C isn't initialized and contains a NaN, then that NaN can still be propagated through to the output when beta == 0.0.

blaspp/include/blas/gemm.hh

Lines 166 to 265 in e954a9b

	// alpha != zero
	if (transA == Op::NoTrans) {
	if (transB == Op::NoTrans) {
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i)
	C(i, j) *= beta;
	for (int64_t l = 0; l < k; ++l) {
	scalar_t alpha_Blj = alpha*B(l, j);
	for (int64_t i = 0; i < m; ++i)
	C(i, j) += A(i, l)*alpha_Blj;
	}
	}
	}
	else if (transB == Op::Trans) {
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i)
	C(i, j) *= beta;
	for (int64_t l = 0; l < k; ++l) {
	scalar_t alpha_Bjl = alpha*B(j, l);
	for (int64_t i = 0; i < m; ++i)
	C(i, j) += A(i, l)*alpha_Bjl;
	}
	}
	}
	else { // transB == Op::ConjTrans
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i)
	C(i, j) *= beta;
	for (int64_t l = 0; l < k; ++l) {
	scalar_t alpha_Bjl = alpha*conj(B(j, l));
	for (int64_t i = 0; i < m; ++i)
	C(i, j) += A(i, l)*alpha_Bjl;
	}
	}
	}
	}
	else if (transA == Op::Trans) {
	if (transB == Op::NoTrans) {
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i) {
	scalar_t sum = zero;
	for (int64_t l = 0; l < k; ++l)
	sum += A(l, i)*B(l, j);
	C(i, j) = alpha*sum + beta*C(i, j);
	}
	}
	}
	else if (transB == Op::Trans) {
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i) {
	scalar_t sum = zero;
	for (int64_t l = 0; l < k; ++l)
	sum += A(l, i)*B(j, l);
	C(i, j) = alpha*sum + beta*C(i, j);
	}
	}
	}
	else { // transB == Op::ConjTrans
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i) {
	scalar_t sum = zero;
	for (int64_t l = 0; l < k; ++l)
	sum += A(l, i)*conj(B(j, l));
	C(i, j) = alpha*sum + beta*C(i, j);
	}
	}
	}
	}
	else { // transA == Op::ConjTrans
	if (transB == Op::NoTrans) {
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i) {
	scalar_t sum = zero;
	for (int64_t l = 0; l < k; ++l)
	sum += conj(A(l, i))*B(l, j);
	C(i, j) = alpha*sum + beta*C(i, j);
	}
	}
	}
	else if (transB == Op::Trans) {
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i) {
	scalar_t sum = zero;
	for (int64_t l = 0; l < k; ++l)
	sum += conj(A(l, i))*B(j, l);
	C(i, j) = alpha*sum + beta*C(i, j);
	}
	}
	}
	else { // transB == Op::ConjTrans
	for (int64_t j = 0; j < n; ++j) {
	for (int64_t i = 0; i < m; ++i) {
	scalar_t sum = zero;
	for (int64_t l = 0; l < k; ++l)
	sum += A(l, i)*B(j, l); // little improvement here
	C(i, j) = alpha*conj(sum) + beta*C(i, j);
	}
	}
	}
	}

@weslleyspereira, recently I suggested that tests for this situation be added to your BLAS stress tests. Once you add such tests, it would be good to run against BLAS++. It might be that other fully-templated functions in BLAS++ make similar mistakes.

[weslleyspereira]

Hi. Good catch, @rileyjmurray !

This test is already in testBLAS. See
https://github.com/tlapack/testBLAS/blob/b177c5b3ae3dbbaece254727e036f487d2271506/src/test_corner_cases.cpp#L993-L997
You can use testBLAS to test the templates from BLAS++. It is prepared for that.

That said, I am not sure what fix should be considered in BLAS++. In <T>LAPACK we adopted a different strategy. There, we have 2 interfaces for gemm:

1. C := alpha*A*B+beta*C
2. C := alpha*A*B

See https://github.com/tlapack/tlapack/blob/master/include/tlapack/blas/gemm.hpp

[rileyjmurray]

The fix in BLAS++ should just be to write C(i, j) = 0.0 for all (i, j) as a preliminary step if beta is zero. This is the only thing that would conform with the BLAS standard, as I understand it. It could even be performed at the start of the function and subsequent code could be left unchanged relative to its current form (note: that would mean $O(n^2)$ multiplications are wasted when beta == 0.0, but that's really not something we should care about).

[weslleyspereira]

I agree with you. Moreover, this fix should be applied to all lv3 BLAS routines.

[mgates3]

Depending on the case, reference BLAS either sets C = 0 before, or sets C = alpha*A*B after (without beta*C). Not really sure why they chose 2 different implementations. The later one has less memory traffic (one less write of C and read of C).

Also, if alpha = 0, it should just do beta*C as reference BLAS does, no need to do O(n^3) work in alpha*A*B. Though I argue that using gemm to scale or initialize a matrix is an abuse of gemm; use lascl or laset.

[weslleyspereira]

Also, if alpha = 0, it should just do beta*C as reference BLAS does, no need to do O(n^3) work in alpha*A*B. Though I argue that using gemm to scale or initialize a matrix is an abuse of gemm; use lascl or laset.

I agree with your argument. One should use lascl or laset if alpha is known to be zero. The current implementation of gemm in <T>LAPACK will do the O(n^3) operation alpha*A*B even if alpha is zero.

[weslleyspereira]

I noticed that I wasn't testing BLAS++ with the option BLAS_USE_TEMPLATE. Now the automatic tests in testBLAS are able to catch the issues related to this thread. See https://github.com/tlapack/testBLAS/actions/runs/4484673095/jobs/7885422971. [Edit:] The output of the tests is currently too verbose. It is easy to look at the raw log of the test.

One of the failures:

2023-03-21T23:14:00.5093998Z  83/172 Test  #83: gemm satisfies all corner cases - float ...................................................................***Failed    0.00 sec
2023-03-21T23:14:00.5094720Z Filters: gemm satisfies all corner cases - float
2023-03-21T23:14:00.5095406Z Randomness seeded to: 555933941
2023-03-21T23:14:00.5095731Z 
2023-03-21T23:14:00.5096228Z ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2023-03-21T23:14:00.5096660Z tester is a Catch2 v3.0.1 host application.
2023-03-21T23:14:00.5097312Z Run with -? for options
2023-03-21T23:14:00.5097616Z 
2023-03-21T23:14:00.5098247Z -------------------------------------------------------------------------------
2023-03-21T23:14:00.5098925Z gemm satisfies all corner cases - float
2023-03-21T23:14:00.5099551Z   C := beta C if m, n > 0 and k = 0
2023-03-21T23:14:00.5100133Z -------------------------------------------------------------------------------
2023-03-21T23:14:00.5100893Z /home/runner/work/testBLAS/testBLAS/src/test_corner_cases.cpp:987
2023-03-21T23:14:00.5101353Z ...............................................................................
2023-03-21T23:14:00.5101875Z 
2023-03-21T23:14:00.5102201Z /home/runner/work/testBLAS/testBLAS/src/test_corner_cases.cpp:990: FAILED:
2023-03-21T23:14:00.5102844Z   CHECK( C[0] == C11 )
2023-03-21T23:14:00.5103230Z with expansion:
2023-03-21T23:14:00.5104610Z   1.0f == 2.0f
2023-03-21T23:14:00.5104943Z 
2023-03-21T23:14:00.5105584Z -------------------------------------------------------------------------------
2023-03-21T23:14:00.5106147Z gemm satisfies all corner cases - float
2023-03-21T23:14:00.5106784Z   C does not need to be set if beta = 0
2023-03-21T23:14:00.5107345Z -------------------------------------------------------------------------------
2023-03-21T23:14:00.5108125Z /home/runner/work/testBLAS/testBLAS/src/test_corner_cases.cpp:993
2023-03-21T23:14:00.5108604Z ...............................................................................
2023-03-21T23:14:00.5109297Z 
2023-03-21T23:14:00.5109626Z /home/runner/work/testBLAS/testBLAS/src/test_corner_cases.cpp:996: FAILED:
2023-03-21T23:14:00.5110306Z   CHECK( (!isnan(C_[0]) && !isnan(C_[1]) && !isnan(C_[2]) && !isnan(C_[3])) )
2023-03-21T23:14:00.5110735Z with expansion:
2023-03-21T23:14:00.5111291Z   false
2023-03-21T23:14:00.5111559Z 
2023-03-21T23:14:00.5112026Z ===============================================================================
2023-03-21T23:14:00.5112426Z test cases:  1 |  0 passed | 1 failed
2023-03-21T23:14:00.5113036Z assertions: 12 | 10 passed | 2 failed

[mgates3]

@weslleyspereira Which option? Your post says: option ``.

[weslleyspereira]

@weslleyspereira Which option? Your post says: option ``.

I'm sorry. I meant BLAS_USE_TEMPLATE.

[mgates3]

Fixed in #99.