[petsc-dev] Mac OSX complex BLAS problem

[Barry Smith]

    This is a known problem. Intel changed their complex dot a while ago causing grief.  Perhaps Apple follows the Intel blas libraries since they are in bed with Intel. 

     Anyways, yes a proper set of configure tests is needed for different versions of complex (double and single) dot and then the PETSc source code needs adjusting to handle the situations. It is likely this means ugly ifdef crap in our source. Thanks Jack we really appreciate it.

   Barry



	From: 	Alexander Grayver <agrayver at gfz-potsdam.de>
	Subject: 	Re: [petsc-users] MKL BLAS interface inconsistency
	Date: 	January 19, 2012 10:55:21 AM CST
	To: 	PETSc users list <petsc-users at mcs.anl.gov>
	Reply-To: 	PETSc users list <petsc-users at mcs.anl.gov>

This (dirty) patch solves the problem:

diff -r a3e9ca59ab58 src/vec/vec/impls/seq/bvec2.c
--- a/src/vec/vec/impls/seq/bvec2.c     Tue Jan 17 22:04:05 2012 -0600
+++ b/src/vec/vec/impls/seq/bvec2.c     Thu Jan 19 17:28:39 2012 +0100
@@ -232,12 +232,13 @@
  PetscErrorCode    ierr;
  PetscInt          n = xin->map->n;
  PetscBLASInt      one = 1, bn = PetscBLASIntCast(n);
+  PetscScalar cnorm;

  PetscFunctionBegin;
  if (type == NORM_2 || type == NORM_FROBENIUS) {
    ierr = VecGetArrayRead(xin,&xx);CHKERRQ(ierr);
-    *z = BLASdot_(&bn,xx,&one,xx,&one);
-    *z = PetscSqrtReal(*z);
+    zdotc(&cnorm,&bn,xx,&one,xx,&one);
+    *z = PetscSqrtReal(PetscAbsScalar(cnorm));
    ierr = VecRestoreArrayRead(xin,&xx);CHKERRQ(ierr);
    ierr = PetscLogFlops(PetscMax(2.0*n-1,0.0));CHKERRQ(ierr);
  } else if (type == NORM_INFINITY) {

The same is applied to mpi vector implementation from /petsc-dev/src/vec/vec/impls/mpi/pvec2.c
Of course it works only if one uses Intel MKL BLAS/LAPACK and double complex arithmetics. Which is my case.


On Feb 18, 2012, at 5:47 PM, Matthew Knepley wrote:

> On Sat, Feb 18, 2012 at 5:34 PM, Jack Poulson <jack.poulson at gmail.com> wrote:
> I think I mentioned this to several PETSc folks two years ago. It is a bad idea to return complex numbers from Fortran, as there is no portable solution as far as I know; some compilers need you to declare it as 
> 
> complex zdotc_( <regular args> )
> 
> some as
> 
> void zdotc_( <regular args>, complex* z )
> 
> and I think I ran into problems with both on OSX (but don't quote me on it). I ended up settling on writing my own complex dotc and dotu routines.
> 
> That is exactly it. Scipy has a patch for it
> 
>   http://projects.scipy.org/scipy/attachment/ticket/1321/scipy-r6856-isolve-veclib.patch
> 
> So what are we going to do in PETSc folks? It sounds like we need a configure test for complex dot, and then
> and then an #define wrapper for ones that return rather than taking an argument.
> 
>     MAtt
>  
> 
> Jack
> 
> 
> On Sat, Feb 18, 2012 at 5:28 PM, Matthew Knepley <knepley at gmail.com> wrote:
> We have a problem with complex dot product on OSX. Here is the easiest way to see it (I think):
> 
> knepley:/PETSc3/petsc/petsc-dev$ ./arch-complex-fftw-debug/lib/ex5-obj/ex5 -snes_monitor
> ./arch-complex-fftw-debug/lib/ex5-obj/ex5 -snes_monitor
>   0 SNES Function norm 0.000000000000e+00 
> 
> With the debugger:
> 
> Breakpoint 1, VecNorm_Seq (xin=0x1022fc170, type=NORM_2, z=0x7fff5fbfea10) at bvec2.c:238
> 238         ierr = VecGetArrayRead(xin,&xx);CHKERRQ(ierr);
> (gdb) n
> Current language:  auto; currently c++
> 239         *z = BLASdot_(&bn,xx,&one,xx,&one);
> (gdb) p bn
> $1 = 16
> (gdb) p xx[0]@16
> $2 = {{
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>    _M_value = -0.10378182131158753 + 0 * I
>   }, {
>     _M_value = -0.10378182131158753 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = -0.10378182131158753 + 0 * I
>   }, {
>     _M_value = -0.10378182131158753 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }, {
>     _M_value = 0 + 0 * I
>   }}
> (gdb) p one
> $3 = 1
> (gdb) n
> 240         *z = PetscSqrtReal(*z);
> (gdb) p *z
> $4 = 0
> (gdb) p zdotc
> $5 = {<text variable, no debug info>} 0x7fff863e71e2 <zdotc_>
> 
> It appears that the dot product is producing nothing. Didn't we have a similar problem a little while ago?
> 
>    Matt
> 
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener
> 
> 
> 
> 
> -- 
> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.
> -- Norbert Wiener