[petsc-users] MKL BLAS interface inconsistency

[Alexander Grayver]

Hello,

The inconsistency is still there. I pulled petsc-dev today and came 
across with the same problem:

src/vec/vec/impls/seq/bvec2.c
239:    *z = PetscRealPart(BLASdot_(&bn,xx,&one,xx,&one));
240:    *z = PetscSqrtReal(*z);

This doesn't work with Intel MKL. However in one of the mailing list 
discussions it was said there is now check for that.
This check didn't work for me apparently. Any ideas why?

Regards,
Alexander

On 19.01.2012 17:55, Alexander Grayver wrote:
> 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