bugGNU Octave - Bugs: bug #32120, complex times complex failure

 
 

bug #32120: complex times complex failure

Submitter:  Pascal Dupuis <cdemills>
Submitted:  Thu 13 Jan 2011 01:53:09 PM UTC
   
 
Category:  Libraries Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Works For Me Assigned to:  None
Originator Name:  Pascal Dupuis Open/Closed:  * Closed
Release:  * 3.2.3 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 27 Mar 2014 05:18:48 AM UTC, comment #37: 

Ok, closing bug report as requested in comment #36.

Mike Miller <mtmiller>
Group Member
Sat 22 Mar 2014 03:31:00 PM UTC, comment #36: 

Folks, I would now recommend that this bug be closed all over again.

It turns out that the problems we're seeing are largely caused by an inconsistent treatment of BLAS across the various libraries that Octave depends upon: QRUpdate, ARPACK.

For instance, if you simply compile QRUpdate and ARPACK with -ff2c, and Octave with -ff2c as well, everything works fine. And since ARPACK and QRUpdate export only subroutines, compiling them with -ff2c does not break any GNU Fortran code that links with them. They work both sides of the proverbial street.

Problems arise when using "shims" to fix the BLAS issues on the Mac, such as "blaswrap.c" here in Octave or "dotwrp" out on GitHub. If Octave is configured not to use dotwrp, but QRUpdate includes it statically (which is currently the case on Homebrew), Octave will break in exactly the fashion described here.

Bottom line: this isn't really Octave's fault. It should be documented heavily, and the package managers may need fixing. I'm working hard on the Homebrew side of things. Hopefully we'll get this nailed down soon.

I've led you on a couple of wild goose chases here. I apologize. If it is any consolation, I've created a bunch of brand new code to fix BLAS/LAPACK incompatibilities that I'm now recommending nobody use.

Michael C. Grant <mcgrant>
Thu 20 Mar 2014 06:18:32 AM UTC, comment #35: 

[Damnit, I keep forgetting the bug tracker has no email gateway. I sure wish that could be changed.]

Yes, it does. It needs the first ax_blas pass to fail, and then it should try "blaswrap.o -framework vecLib".

I've actually learned a lot about the BLAS issue, and about dynamic library handling in OSX, studying this. Armed with this info I've basically written an entirely new blas wrapper (I'm calling it vecLibFort). It fixes all of the offending BLAS and LAPACK calls (hopefully all of them), using a combination of static and dynamic replacement. It also provides the ability to patch compiled applications using the Mac equivalent of LD_PRELOAD (DYLD_INSERT_LIBRARIES). It's also a lot more compact thanks to some preprocessing tricks.

I think it's something that would be of great benefit to the wider Mac community; I know there are other Homebrew projects that have latent bugs due to the gfortran/vecLib incompatibility.

If you'll hang just a bit longer, I'll let you know when that code is ready to check out (likely tomorrow).

Michael C. Grant <mcgrant>
Thu 20 Mar 2014 05:18:50 AM UTC, comment #34: 

At first glance this change looks OK to me.  Does it also need the earlier change in the configure script that I proposed (to reorder the test for OS X veclib and the check for BLAS with -ff2c)?

John W. Eaton <jwe>
Group administrator
Wed 19 Mar 2014 05:10:08 PM UTC, comment #33: 

I've been toying around with blaswrap.c and I think I have a better solution, which I've  attached. What it does is make each function responsible for making its own call to dlsym(), and storing its function pointer in a local static variable.

I came up with this based on some reading on the interwebs about interposing. I tried replacing the direct vecLib reference with RTLD_NEXT; alas, that does not seem to work, for reasons I cannot determine (because according to my reading of the documentation, it should). That's too bad, because it would enable us to remove the construct and destructor calls altogether.

(file #30947)

Michael C. Grant <mcgrant>
Tue 18 Mar 2014 03:13:59 AM UTC, comment #32: 

I've attached a replacement blaswrap.c. I removed the abort() commands and added sdot_. Obviously you may not want to remove the abort()s permanently without understanding why they are being triggered, but at least this replacement blaswrap.c now passes all the tests where these functions are involved.

I answered my own question re: ARPACK I think; this doesn't seem to help that.

(file #30935)

Michael C. Grant <mcgrant>
Tue 18 Mar 2014 01:24:19 AM UTC, comment #31: 

Even if I bypass the "abort" errors, "sdot" still fails the configure test. To fix this, you need to add "sdot_" to LAPACK_LIST(_). (Yes,  I know that it's not LAPACK, but that list does what we want :-))

I'm coming around to the idea that this dynamic wrapper is the best approach, because I've noticed that ARPACK calls suffer as well from the BLAS issue. It seems to me this fixes it all around. Does this wrapper handle the needs of these other libraries as well, or just liboctave?

Michael C. Grant <mcgrant>
Mon 17 Mar 2014 11:40:11 PM UTC, comment #30: 

I did attempt to use USE_BLASWRAP. I got an "abort()" error, which I traced to a dlsym() error for "slamc3_", which is odd, because that function does exist. I bypassed the abort() call and was able to get Octave to load, but I've not done further testing yet. I will continue.

Michael C. Grant <mcgrant>
Mon 17 Mar 2014 10:54:50 PM UTC, comment #29: 

I'd still like to know whether the existing blaswrap.c file works properly if -ff2c is not used and USE_BLASWRAP is defined.  If so, then I'd probably just fix the configure script on the stable branch and then work on better ways to interface with BLAS and LAPACK libraries on the default branch.

John W. Eaton <jwe>
Group administrator
Mon 17 Mar 2014 10:52:27 PM UTC, comment #28: 

Yes, the r1mach problem looks like a cut and paste error.  I checked in the following changeset for that problem:

http://hg.savannah.gnu.org/hgweb/octave/rev/c08776badd3d

John W. Eaton <jwe>
Group administrator
Mon 17 Mar 2014 08:19:28 PM UTC, comment #27: 

I'm making some solid progress here.

As an experiment, I ran configure with --with-blas and --with-lapack set to "-ldotwrp -framework vecLib". "-ldotwrp" is the wrapper library that I linked to below, that rewrites the 5 dot product functions that cause problems for BLAS.

With this setting, the existing configure passes all of the BLAS tests without -ff2c, and the resulting executable passes all of the dot-product related Octave tests.

However, new failures in several SLATEC functions arise. As far as I can tell, they are caused specifically by "r1mach.f" and its call to the LAPACK function "slamch". I was able to fix these errors, however, by editing r1mach.f and changing it from a "double precision" function to a "real" function. Looking at the SLATEC calls that use "r1mach", that seems to be the correct signature for R1MACH anyway.

For reference, I'm attaching the "clapack.h" header file supplied by Apple. You can see that it does follow the convention of returning "double" even if the underlying function is single precision.

(file #30930)

Michael C. Grant <mcgrant>
Mon 17 Mar 2014 05:57:44 PM UTC, comment #26: 

Perhaps the best bet is this: create a new "blaswrap.c" that basically does the same thing, but instead of rewriting the functions dynamically, it rewrites them statically---using the CBLAS calls to avoid the inevitable name conflicts.

If you want I can work on this.

Michael C. Grant <mcgrant>
Mon 17 Mar 2014 02:53:59 PM UTC, comment #25: 

OK, I didn't realize that the cblas interface was included in the apple framework thing.

Given that info, let me think about how I'd prefer to clean this up so that it might be more reliable and simpler.

John W. Eaton <jwe>
Group administrator
Mon 17 Mar 2014 02:41:29 PM UTC, comment #24: 

I will give this a shot when I get home in a 2-3 hours. But from what I'm reading in blaswrap.c, that code is broken for 64-bit Matlab anyway. I am glad you explained it though, because when I hard-coded USE_BLASWRAP I did not disable -ff2c, so that was probably why it didn't work for me.

I'm honestly not sure what you're asking re: Apple's framework BLAS and f2c. The CBLAS calls are part of the framework, so there's no need to compile it. I'm not sure how it is implemented internally---that is, if the Fortran interface calls CBLAS or vice versa. But I am not sure we should care.

I think one reason why the configure script was never right is because the folks doing Mac ports, like Homebrew, etc., were either using OpenBLAS or the "dotwrp" workaround like I linked to below. I know for a fact that Homebrew used that previously, for instance; but it was configured to work only with a particular version of OSX. So the configure script may not have been exercised this way before.

Michael C. Grant <mcgrant>
Mon 17 Mar 2014 02:35:48 PM UTC, comment #23: 

Sorry, I forgot to attach the patch before I hit submit for my last comment.

(file #30926)

John W. Eaton <jwe>
Group administrator
Mon 17 Mar 2014 02:34:19 PM UTC, comment #22: 

I agree the blaswrap solution is not great, but I'm thinking that the configure script has never been right.

Could you please try a fresh build with the attached patch?  After running bootstrap, you should have USE_BLASWRAP defined and FFLAGS should NOT contain -ff2c.  If that's not what happens, you can skip the build and just let me know.  Or edit config.h and all the Makefiles to remove -ff2c from FFLAGS and try building.

Again, what does cblas do to handle this situation if the apple framework blas is compiled with f2c conventions and it is just using it, not replacing it?  I think we should be able to do the same and avoid adding cblas as a dependency.

John W. Eaton <jwe>
Group administrator
Mon 17 Mar 2014 12:13:35 PM UTC, comment #21: 

Without the proposed changes, USE_BLASWRAP is not defined. That's because, for some reason, the -ff2c test passes, even for ZDOTU. I think the way configure is built, it expects the -ff2c test to fail. And you would think it would, given the experience I am having with it, but it doesn't.

The dynamic linker tricks employed by USE_BLASWRAP seems like a horribly complex and fragile solution. I tried forcibly adding USE_BLASWRAP to config.h, and octave immediately aborted. I tried forcibly running the USE_BLASWRAP tests in the configure script, and they aborted as well.

I don't think the dynamic linker trick is necessary. I did some searching, and it looks like previous Homebrew recipes for Octave solved the problem by providing alternative versions of SDOT, ZDOTU, ZDOTC, CDOTU, and CDOTC that wrap the CBLAS versions. This gives you the freedom to stick with standard calls to BLAS. Here's the wrapper code they use:

https://github.com/tenomoto/dotwrp/blob/master/dotwrp.f90

Perhaps I can try and figure out how to add this into cruft in a clean way.

Michael C. Grant <mcgrant>
Mon 17 Mar 2014 04:48:01 AM UTC, comment #20: 

Without any of the proposed changes, is USE_BLASWRAP defined in config.h?  Does FFLAGS include -ff2c?

On OS X, I think the intent is to have USE_BLASWRAP defined and not add -ff2c to FFLAGS, but I'm guessing that the ordering of the configure tests is wrong and we are getting just the opposite, with -ff2c in FFLAGS and USE_BLASWRAP not defined.

John W. Eaton <jwe>
Group administrator
Mon 17 Mar 2014 04:35:53 AM UTC, comment #19: 

Sincd people often replace the blas library that Octave is linked to with some accelerated version, maybe the configure test is not all we need here.  We can get that test right and then someone can use dynamic linker tricks or just replace the system library with some other version and still be screwed.  Maybe we need to also perform the configure checks when Octave starts and refuse to run if the blas library looks like it is not working properly?

Is Octave being compiled with -ff2c?  If so, then I think all Fortran libraries will need to be compiled this way.  Are they?

Looking at the configure tests, it seems that the wrapper that we are using on OS X systems is supposed to work without using -ff2c.

Is Apple's veclib compiled with f2c conventions or not?  the comments in blaswrap.c say that it is, but if we are using that code, we should not be using -ff2c to compile the Fortran bits of Octave.

If we use cblas, is it still possible for people to replace the blas library by using dynamic linker tricks or by replacing the blas and lapack shared libraries, or do they then have to also replace the clas library?  I guess this question is really asking whether cblas is just a wrapper library or if it actually provide the implementation of the blas library.

John W. Eaton <jwe>
Group administrator
Mon 17 Mar 2014 01:40:15 AM UTC, comment #18: 

OK, this is only a partial fix. It turns out that you need to do something similar in fMatrix.cc, in order to handle SDOT properly; otherwise you get a FAIL result there, too...

...and in libinterp/corefcn/dot.c as well. That's a much bigger one. I wouldn't be surprised if every complex and single-precision dot product function in that file is broken here.

ARGH! There has got to be a better way.

Michael C. Grant <mcgrant>
Mon 17 Mar 2014 01:19:42 AM UTC, comment #17: 

OK, here's a patch that seems to do the trick.

Here's what it does: in ax_blas.m4, it creates a new define USE_CBLAS_COMPLEX that is set to 1 when the "-framework vecBlas" version of the BLAS is detected. This define is then used in CMatrix.cc and fCMatrix.cc to replace the calls to xzdotu, etc. to corresponding CBLAS calls.

I've tested this and found it does indeed produce the correct results for both single and double precision.

It seems to me that it might be a good idea to generalize this: instead of defining USE_CBLAS_COMPLEX for the specific case of vecLib, you could test for the existence of these functions in whatever BLAS library was successfully called, and then define USE_CBLAS_COMPLEX whenever they are found.

By the way, USE_BLASWRAP and the corresponding code in blaswrap.c should probably be dropped. It seems to me this is much less fragile solution.

(file #30922)

Michael C. Grant <mcgrant>
Sun 16 Mar 2014 08:52:02 PM UTC, comment #16: 

Sorry for spamming y'all, but the rabbit hole is getting deeper. Just to make sure there isn't a subtle difference in compilation causing the problem, I linked my last source code below to the exact object file created for xzdotu.f by the Octave build script---and everything works fine.

So for whatever reason, linking my test code below to Octave's xzdotu.o and -framework vecLib works just fine. But xzdotu.o and -framework vecLib in Octave does not work.

I think the best thing to do here is not question why, and just find a workaround, perhaps by calling vecLib's own cblas_* code.

Michael C. Grant <mcgrant>
Sun 16 Mar 2014 07:46:42 PM UTC, comment #15: 

I have to admit I am stumped. Looking at config.log, I see signs that configure is performing the right tests to see what to do. Here is the test that it runs:

       program main

       double complex zdotu,a(1),b(1),w
       external zdotu
       a(1) = dcmplx(1d0,1d0)
       b(1) = dcmplx(1d0,2d0)
       w = zdotu(1,a,1,b,1)
       if (w .ne. a(1)*b(1)) stop 1

       end

This program fails when compiled with "gfortran -framework vecLib", but succeeds when compiled with "gfortran -ff2c -framework vecLib". Just to make sure that this test just "got lucky" because it use n=1, I created this:

       program main

       double complex zdotu,a(3),b(3),w
       external zdotu
       a(1) = dcmplx(1d0,1d0)
       a(2) = dcmplx(2d0,1d0)
       a(3) = dcmplx(3d0,1d0)
       b(1) = dcmplx(4d0,1d0)
       b(2) = dcmplx(5d0,1d0)
       b(3) = dcmplx(6d0,1d0)
       w = zdotu(3,a,1,b,1)
       write(*,*) a(1), b(1), w

       end

Sure enough, this test passes with -ff2c, and fails without it. Then I tried this:

       program main

       double complex a(3),b(3),w
       external xzdotu
       a(1) = dcmplx(1d0,1d0)
       a(2) = dcmplx(2d0,1d0)
       a(3) = dcmplx(3d0,1d0)
       b(1) = dcmplx(4d0,1d0)
       b(2) = dcmplx(5d0,1d0)
       b(3) = dcmplx(6d0,1d0)
       call xzdotu(3,a,1,b,1,w)
       write(*,*) a(1), b(1), w

       end

And I compiled this with the exact code used in liboctave/cruft/blas-xtra. And again, with the -ff2c flag, it works. And this is exactly how it is compiled in Octave:

libtool: compile:  /usr/local/bin/gfortran -ff2c -O -c blas-xtra/xzdotu.f  -fno-common -o blas-xtra/.libs/libcruft_la-xzdotu.o

And yet, I verified, by inserting WRITE statements in xzdotu.f, that it is the call to zdotu that is failing within Octave.

The wrapper is working in isolation, when compiled using the same flags Octave uses, but it's not working within Octave itself. I don't get it.

Michael C. Grant <mcgrant>
Sun 16 Mar 2014 02:01:50 PM UTC, comment #14: 

I will try. But it looks like the existing wrappers ought to work if only the right flags are set. Do you agree? Is it worth trying to use them, or should I figure out how to swap in the cblas_* routines instead?

Michael C. Grant <mcgrant>
Thu 13 Mar 2014 02:38:38 PM UTC, comment #13: 

Looks like we need an autoconf check to avoid using the broken Apple Accelerate. Can you write it?

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Wed 12 Mar 2014 11:07:53 PM UTC, comment #12: 

I was able to achieve correct operation by forcibly replacing the ZDOTU call with a call to "cblas_zdotu_sub" from vecLib.

Now, I don't know the right way to patch this up, but at least I've narrowed it down... I can do some work on this if someone will give me direction. Perhaps this is just a configure error?

Michael C. Grant <mcgrant>
Wed 12 Mar 2014 10:02:55 PM UTC, comment #11: 

The current version of the Octave Homebrew recipe has this:


  # When building 64-bit binaries on Snow Leopard, there are naming issues with
  # the dot product functions in the BLAS library provided by Apple's
  # Accelerate framework. See the following thread for the gory details:
  #   http://www.macresearch.org/lapackblas-fortran-106
  # We can work around the issues using dotwrp.
  depends_on 'dotwrp' if MacOS.version == :snow_leopard and MacOS.prefer_64_bit?


I am not currently using this library, and I suspect if I were, this problem would go away. But it seems to me this is exactly the same issue that XZDOTU, etc. are tasked to handle. So perhaps the configure script isn't doing the right thing on the mac?

vecLib actually comes with subroutine versions of these four routines; so in theory, they could be used in place of the XZDOTU, etc. I think they might be drop-in replacements to what you have now.

Michael C. Grant <mcgrant>
Wed 12 Mar 2014 09:52:03 PM UTC, comment #10: 

I'm afraid I can't---valgrind on Mac OS 10.8/10.9 is "experimental and mostly broken". I tried anyway, but it died before it ever began the test.

Michael C. Grant <mcgrant>
Wed 12 Mar 2014 09:45:17 PM UTC, comment #9: 

I will try. But it's well known that the BLAS routines that return a DOUBLE COMPLEX are problematic. I think we have good reason to believe that the issue is with the XZDOTU wrapper approach and/or the configure script.

Michael C. Grant <mcgrant>
Wed 12 Mar 2014 09:42:46 PM UTC, comment #8: 

Those numbers look like uninitialised memory. Can you run Octave under valgrind?

./run-octave -valgrind samplescript.m

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Wed 12 Mar 2014 09:40:41 PM UTC, comment #7: 

I'm not surprised. This seems to be specific to certain BLAS implementations. I'm using vecLib on the Mac at the moment.

There's quite a bit of configure cruft devoted to testing the correctness of ZDOTU, which is the offending function here, and other BLAS functions.

Michael C. Grant <mcgrant>
Wed 12 Mar 2014 09:31:31 PM UTC, comment #6: 

ubuntu 12.04 works correctly.
x = [1 + i, 2 + i, 3 + i]
x =

   1 + 1i   2 + 1i   3 + 1i

y = [4 + i; 5 + i; 6 + i]
y =

   4 + 1i
   5 + 1i
   6 + 1i

x * y
ans =  29 + 21i

Doug Stewart <dastew>
Wed 12 Mar 2014 09:14:41 PM UTC, comment #5: 

Looks like it's the ZDOTU issue, actually, that a search of your mailing lists reveals has been a thorn in your side for years.

Michael C. Grant <mcgrant>
Wed 12 Mar 2014 09:01:35 PM UTC, comment #4: 

Can we consider re-opening this? I'm seeing a bug on this very test of CMatrix.cc when building gui-release on the Mac using the native BLAS. Now, it could very well be that the bug lies in the BLAS. But I would suggest that this particular build case is important enough that we should figure out why it is happening and attempt to work around it.

>> (1+i)*(4+i)+(2+i)*(5+i)+(3+i)*(6+i)
ans =  29 + 21i
>> x = [1+i 2+i 3+i]
x =

   1 + 1i   2 + 1i   3 + 1i

>> y = [ 4+i ; 5+i ; 6+i]
y =

   4 + 1i
   5 + 1i
   6 + 1i

>> x*y
ans =  6.9263e-310 + 6.9263e-310i
>> sum(x.'.*y)
ans =  29 + 21i


Michael C. Grant <mcgrant>
Fri 14 Jan 2011 12:28:41 AM UTC, comment #3: 

If, as I think you posted on the help list, you can make this problem go away by linking with a different BLAS library, then I don't think this is a bug in Octave.  So I'm closing this report.

John W. Eaton <jwe>
Group administrator
Thu 13 Jan 2011 04:01:32 PM UTC, comment #2: 

and running
doug@doug-desktop:~/octave_dev/main/octave$ hg tip
changeset:   11514:a1386fab2a02
tag:         tip
user:        John W. Eaton <jwe@octave.org>
date:        Thu Jan 13 08:04:23 2011 -0500
files:       ChangeLog configure.ac

I get

octave:1>  [1+i 2+i ; 3+i 4+i ] * [5+i ; 6+i]
ans =

   15 + 14i
   37 + 18i

octave:2>

Doug Stewart

Anonymous
Thu 13 Jan 2011 02:19:10 PM UTC, comment #1: 

in 3.2.3 I get

>>> [1+i 2+i ; 3+i 4+i ] * [5+i ; 6+i]

ans =

   15 + 14i
   37 + 18i

>>>


Doug Stewart

Anonymous
Thu 13 Jan 2011 01:53:09 PM UTC, original submission:  

This bug is present in both 3.2.3 and 3.3.54+. I tried to execute the test found in liboctave/CMatrix.cc

assert([1+i 2+i ; 3+i 4+i ] * [5+i ; 6+i], [15 + 14i ; 37 + 18i], 1e-14)

It fails,  and the funny thing is:
[1+i 2+i] * [5+i ; 6+i]
ans =  15 + 14i %#  OK
[3+i 4+i ] * [5+i ; 6+i]
ans =  37 + 18i %# OK

[1+i 2+i ; 3+i 4+i ] * [5+i ; 6+i]
ans =

   17 + 11i
   39 + 11i

Just plain wrong.

Regards

Pascal

Pascal Dupuis <cdemills>
Group Member

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #30947:  blaswrap.c added by mcgrant (5KiB - application/octet-stream)
file #30935:  blaswrap.c added by mcgrant (10KiB - application/octet-stream)
file #30930:  clapack.h added by mcgrant (429KiB - application/octet-stream)
file #30926:  diffs.txt added by jwe (1KiB - text/plain)
file #30922:  cblas.diff added by mcgrant (5KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by dastew (Posted a comment)
  • -email is unavailable- added by mcgrant (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by cdemills (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-03-27 mtmiller StatusConfirmed Works For Me
        Open/ClosedOpen Closed
    2014-03-19 mcgrant Attached File- Added blaswrap.c, #30947
    2014-03-18 mcgrant Attached File- Added blaswrap.c, #30935
    2014-03-17 mcgrant Attached File- Added clapack.h, #30930
    2014-03-17 jwe Attached File- Added diffs.txt, #30926
    2014-03-17 mcgrant Attached File- Added cblas.diff, #30922
    2014-03-12 jordigh StatusWorks For Me Confirmed
        Open/ClosedClosed Open
    2011-01-14 jwe StatusNone Works For Me
        Open/ClosedOpen Closed

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code