Hi Michael,
I think I've found the source of the complex matrix multiplication problem I'm having. The xgemm function in CMatrix.h tests a couple cases before executing a LAPACK function to do the multiplication. For most cases, ZGEMM is called and this works correctly in all builds. However, if a 1xn matrix A is multiplied by a nx1 matrix B, as in c = A*B, XZDOTU is called instead and in the VS2010 build the result is not always (ever?) correct. The octfile below should reproduce the problem.
It seems that this is a problem
with the BLAS/LAPACK/ATLAS package being used in the VS2010 build. Does this sound correct? (I've never understood the difference between these three things, so maybe now is a good time to learn.) There are some posts from 2007/08 that touch on this issue with XZDOTU. Any information you might have about the packages used in the build (or how to find that information) would be greatly appreciated.
MTP
#include <octave/oct.h>
DEFUN_DLD (complex_mult_test_oct, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} address@hidden, @var{b}, @var{c}, @var{d}] =} complex_mult_test_oct ()\n\
Test complex matrix multiplication.\n\
@end deftypefn")
{
int nargin = args.length ();
octave_value_list retval;
if (nargin != 0)
{
print_usage();
}
else
{
ComplexMatrix a(2,2);
a(0,0) = Complex(0.0,0.0);
a(0,1) = Complex(1.0,0.0);
a(1,0) = Complex(0.0,0.0);
a(1,1) = Complex(1.0,0.0);
ComplexMatrix b(2,2);
b(0,0) = Complex(0.0,0.0);
b(0,1) = Complex(0.0,0.0);
b(1,0) = Complex(1.0,0.0);
b(1,1) = Complex(1.0,0.0);
/******** Test 1 *********/
//The following compiles and executes correctly on both the MINGW and VS2010 builds.
//The multiplication operator in this case calls xgemm->ZGEMM.
ComplexMatrix c = a*b;
/******** Test 2 *********/
//The following compiles on VS2010 build, but does not execute correctly.
//The multiplication operator in this case calls xgemm->XZDOTU, and XZDOTU seems to have a problem in the VS2010 build.
ComplexMatrix d = a.extract(0,0,0,1)*b.extract(0,0,1,0);
retval.append(a);
retval.append(b);
retval.append(c);
retval.append(d);
}
return retval;
}
From: Mike Puglia <address@hidden>
To: Michael Goffioul <address@hidden>
Cc: "address@hidden" <address@hidden>
Sent: Saturday, September 21, 2013 11:00 AM
Subject: Re: Octave 3.6.4 VS2010 and the C++ API
Hi,
I have found another seemingly major difference between between the VS2010 and MINGW builds, and this would seem to be a bug, rather than a difference in the API definition. I just wanted to confirm before looking further into this.
The oct file below is a simple demonstration of a
complex matrix multiplication. In the MINGW build, it correctly results in c = 1; In the VS2010 build it incorrectly results in c = 5.2998e-315 + 2.3146e-245i. I can compile and link in both platforms with no errors or warnings coming from mkoctfile.
Looking at the header definition of the complex matrix multiplication, it would seem that it is exposed to the DLL, since it is decorated with extern
OCTAVE_API. Also, other members of the ComplexMatrix class (like transpose() and conj()) are exposed and do work correctly in the VS2010 build.
Before I look into this further, do you have any thoughts on what might be happening with the complex matrix multiplication?
And following up on the original issue, the work-around for
the xleftdiv() issue below does exist. The solve() method is exposed to the API through the Matrix classes and can be used to write an xleftdiv() function locally with out too much effort (though it is of course messier than having it in the API).
Thanks is advance for your help.
MTP
#include <octave/oct.h>
DEFUN_DLD (complex_mult_test_oct, args, nargout,
"-*- texinfo -*-\n\
@deftypefn {Function File} address@hidden, @var{b}, @var{c}, @var{d}]
=} complex_mult_test_oct ()\n\
Test complex matrix multiplication. c = a*b.\n\
@end deftypefn")
{
int nargin = args.length ();
octave_value_list retval;
if (nargin != 0)
{
print_usage();
}
else
{
ComplexMatrix a(1,2);
a(0,0) = Complex(0.0);
a(0,1) = Complex(1.0);
ComplexMatrix b(2,1);
b(0,0) = Complex(0.0);
b(1,0) = Complex(1.0);
ComplexMatrix c = a*b;
Complex d = Complex(0.0)*Complex(0.0) + Complex(1.0)*Complex(1.0);
retval.append(a);
retval.append(b);
retval.append(c);
retval.append(d);
}
return retval;
}
From: Mike Puglia <address@hidden>
To: Michael Goffioul <address@hidden>
Cc: "address@hidden" <address@hidden>
Sent: Thursday, September 19, 2013 11:01 PM
Subject: Re: Octave 3.6.4 VS2010 and the C++ API
Thank you for the help. You have saved me a great deal.
Off-hand, can you think of anything in the API that is exposed and can achieve the same functionality? Specifically I'll need to somehow replicate:
extern ComplexMatrix xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b, MatrixType &typ, blas_trans_type transt = blas_no_trans);
I'm guessing a complete recompile is going to be a bit beyond my
abilities/resources at the moment. I
saw the scripts and pre-compiled VS2010 packages you've posted at Octave for Windows though. Will those work for v3.6.4? Generally, how
would I proceed with a compilation? I assume I'll need Cygwin to run the sh scripts? Do the scripts work in conjunction with the pre-compiled libraries, so that a complete recompile would not be necessary if only the functions in xdiv.h are decorated/exposed?
Thanks again for the help.
From: Michael Goffioul <address@hidden>
To: Mike Puglia <address@hidden>
Cc: "address@hidden" <address@hidden>
Sent: Thursday, September 19, 2013 10:01 PM
Subject: Re: Octave 3.6.4 VS2010 and the C++ API
_______________________________________________
Help-octave mailing list
address@hiddenhttps://mailman.cae.wisc.edu/listinfo/help-octave