Re: [petsc-dev] Mac OSX complex BLAS problem

18 Feb 2012

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 <[email protected]>
	Subject: 	Re: [petsc-users] MKL BLAS interface inconsistency
	Date: 	January 19, 2012 10:55:21 AM CST
	To: 	PETSc users list <[email protected]>
	Reply-To: 	PETSc users list <[email protected]>

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 <[email protected]> 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-ve...
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 <[email protected]> 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

    

Barry Smith

tags

participants (1)