[petsc-users] MKL BLAS interface inconsistency

[Alexander Grayver]

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 19.01.2012 15:30, Alexander Grayver wrote:
> Dear petsc team,
>
> I've been struggling with this problem for a long time and seeking for 
> your advice now.
> Let's take simple petsc program:
>
> #include <petscksp.h>
> int main(int argc,char **args)
> {
>   Vec                b;
>   PetscReal      norm;
>
>   PetscInitialize(&argc,&args,(char *)0,NULL);
>
>   VecCreateSeq(PETSC_COMM_SELF,100,&b);
>   VecSet(b,2.0);
>   VecNorm(b,NORM_1_AND_2,&norm);
>
>   return 0;
> }
>
> This program works well if I compile petsc with non-mkl blas/lapack.
> However, if I compile petsc with mkl blas/lapack this program crashes:
>
> zdotc,                   FP=7fff8f121e30
> VecNorm_Seq,    FP=7fff8f121f90
> VecNorm_Seq,    FP=7fff8f1220f0
> VecNorm,            FP=7fff8f122230
> main,                   FP=7fff8f122290
>
> It crashes in zdotc. You call this routine as following:
>
> *z = BLASdot_(&bn,xx,&one,xx,&one);
>
> When I look at the zdotc interface in mkl.h I see:
> void    ZDOTC(MKL_Complex16 *pres, const MKL_INT *n, const 
> MKL_Complex16 *x, const MKL_INT *incx, const  MKL_Complex16 *y, const 
> MKL_INT *incy);
>
> I also found example here: 
> http://software.intel.com/en-us/articles/intel-math-kernel-library-intel-mkl-blas-cblas-and-lapack-compilinglinking-functions-fortran-and-cc-calls/
> There are 6 input parameters. But in "classical" BLAS implementation 5 
> input and 1 output.
>
> For example, NORM_1 works well since interface to dzasum is the same.
>
> Any ideas?
>


-- 
Regards,
Alexander