@jwe: I checked the patch in under your credentials here (http://hg.savannah.gnu.org/hgweb/octave/rev/e6b037514f56).

I amended it slightly.  There seems to be a bit of a convention to use do_XXX for the internal function which implements the user visible XXX function.  Hence, I renamed blkmm_ov to do_blkmm.  I also rolled the templated blkmm and blkmm_ov in to the one function do_blkmm because they were both short and it didn't seem necessary to have so many templated functions.  I added a BIST test for this bug as well.

Marking as fixed and closing report.

@jwe: Maybe you can answer this question.  The fix for this bug is conceptually simple, check for empty matrix inputs and return an empty matrix output if that is the case.  However, in order to get the correct output type, double; double complex; single; single complex, it appears like the code to check for an empty matrix would need to be repeated four times.  Simple, but ugly.

My proof of concept working only with double arrays is below

  // Special case of empty matrices
  if (m == 0 || n == 0)
    {
      retval = NDArray (dimz);
      return retval;
    }

  if (argx.iscomplex () || argy.iscomplex ())
    {
      if (argx.is_single_type () || argy.is_single_type ())
        {
          FloatComplexNDArray x = argx.float_complex_array_value ();
          FloatComplexNDArray y = argy.float_complex_array_value ();
          FloatComplexNDArray z (dimz);

          F77_XFCN (cmatm3, CMATM3, (m, n, k, np,
                                     F77_CONST_CMPLX_ARG (x.data ()), F77_CONST_CMPLX_ARG (y.data ()),
                                     F77_CMPLX_ARG (z.fortran_vec ())));
          retval = z;
        }
      else
        {
          ComplexNDArray x = argx.complex_array_value ();
          ComplexNDArray y = argy.complex_array_value ();
          ComplexNDArray z (dimz);

          F77_XFCN (zmatm3, ZMATM3, (m, n, k, np,
                                     F77_CONST_DBLE_CMPLX_ARG (x.data ()), F77_CONST_DBLE_CMPLX_ARG (y.data ()),
                                     F77_DBLE_CMPLX_ARG (z.fortran_vec ())));
          retval = z;
        }
    }

Is there any good way of avoiding the repetition of "if (m ==0 || n == 0)" within each branch of the if/else tree which decodes type?

If not, I'll just duplicate the tests.