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
