topic Re: Bug in GEMV ? in Intel® oneAPI Math Kernel Library

Bug in GEMV ?

OP1 — Thu, 04 Sep 2008 21:39:03 GMT

The following simple program, which uses MKL's GEMV for complex array sections, produce erroneous results.

The program is made of one module (containing one complex array LX), and two subroutines (including the main program). The program links mkl_blas95.lib (which has been created following the MKL documentation). The file mkl_blas.f90 is included in the project.

I am using VS 2005 - Intel Fortran 10.1 and MKL 10.1 Beta 2. The code is compiled as a 32-bit executable and runs on 64-bit XP.

Here is the module:

MODULE
MOD
COMPLEX(8),ALLOCATABLE :: LX(:,:)
END MODULE MOD

Here is the subroutine (which accesses the module, and calls GEMV):

SUBROUTINE
SUB(L,X,Y,N,M)
USE MOD
USE MKL95_BLAS
IMPLICIT NONE
INTEGER(4) N,M
COMPLEX(8),INTENT(IN) :: X(N,M)
COMPLEX(8),INTENT(OUT) :: Y(N,M)
COMPLEX(8),INTENT(IN) :: L(N,N,M)
CALL GEMV(L(1:N,1:N,1),X(1:N,1),LX(1:N,1))
Y = LX
END SUBROUTINE SUB

Here is the main program:

PROGRAM
TEST
USE MOD
IMPLICIT NONE
COMPLEX(8),ALLOCATABLE :: L(:,:,:),X(:,:),Y(:,:)
INTEGER(4) N,M
N=2
M=3
ALLOCATE(L(N,N,M),X(N,M),Y(N,M),LX(N,M))
L = DCMPLX(1.0D+0,0.0D+0)
X = DCMPLX(1.0D+0,0.0D+0)
Y = DCMPLX(0.0D+0,0.0D+0)
CALL SUB(L,X,Y,N,M)
WRITE(*,*) Y
END PROGRAM TEST

The output of TEST is of the 6.277E+66 kind...

Now, if I declare and initializeLX in SUB the program works as intended. Similarly, the program works as is if I use ZGEMV instead of GEMV (with the necessary modifications to SUB).

Any explanation? I am sure that somethin gs escapes my understanding - but could this be a bug as well?

Olivier

Re: Bug in GEMV ?

TimP — Sat, 06 Sep 2008 03:20:35 GMT

I don't have the same compilers installed, as I have little interest in 32-bit compilation on 64-bit OS, but this example seems to work OK here. Interesting combination of legacy non-standard intrinsic with f90.
If you have a problem with the beta compilers offered at the top of the Fortran forum, which should include MKL, you should file a report on premier.intel.com.

Re: Bug in GEMV ?

OP1 — Mon, 08 Sep 2008 14:45:58 GMT

Hi Tim18,

I don't think it's a matter of 32-bit code running on 64-bit machine. I just indicated this for the sake of completeness. I believe there must be something fishy going on with the compiler or GEMV in this very particular case (since it works with other constructs). The problem is with the access in SUB to the module variable LX.

For instance, if I modify the last lines of SUB as follows:

CALL
GEMV(L(1:N,1:N,1),X(1:N,1),LX(1:N,1))
Y = LX

transformed into

CALL GEMV(L(1:N,1:N,1),X(1:N,1),Y(1:N,1))

then the code works flawlessly (the module variables is not accessed anymore). And, as I said, the subroutine ZGEMV works with the code as is (except it has more parameters, of course) - in other words, ZGEMV has no issues accessing properly the desired section of the complex module variable LX, whereas GEMV's behavior is troublesome.

If anybody has an idea about what's going on here, this would be great (especially if somebody can reproduce this behavior). The sample code above is very representative of what my real code does. Right now everything works fine with the ZGEMV subroutine - so this is what I use. But I would appreciate being able to use GEMV... it's cleaner... faster to code... and it's supposed to work.

Tim18, what do you mean when you say it's "a combination of legacy non-standard intrinsic with F90?"

Thanks,

Olivier

Re: Bug in GEMV ?

TimP — Mon, 08 Sep 2008 15:38:16 GMT

dcmplx is non-standard Fortran, such as would have been used 20 years ago. There have been a few discussions about the lack of a generic f77 styleequivalent. I guess one can rely on dcmplx(x,y) to work the same as cmplx(x,y,kind(1d0)), if it works at all. In your style, it would be cmplx(x,y,8).

Your use of kinds 4 and 8of course is non-portable, although it will work with most compilers for the same platforms.

Re: Bug in GEMV ?

Gennady_F_Intel — Tue, 09 Sep 2008 06:12:32 GMT

Oliver.

Id like to add my five cents did you try to initialize array LX before calculation?

--Gennady

Re: Bug in GEMV ?

Vladimir_Koldakov__I — Tue, 09 Sep 2008 06:22:29 GMT

It's really, according to GEMV description,y = Ax + y.

ALLOCATE statement allocates memory only. It does not nullify elements.

Often we have zero elements after ALLOCATE, but in common case undefined values.

Thanks,

Vladimir

Re: Bug in GEMV ?

OP1 — Tue, 09 Sep 2008 12:07:14 GMT

Thanks Gennady and Vladimir... you both were right... the solution to this *problem*is that the default beta value in GEMV is 1... and therefore I should have initialized LX... Arrggh... I hate hidden default values :)

Thanks again!!

Olivier

Re: Bug in GEMV ?

Vladimir_Koldakov__I — Thu, 11 Sep 2008 03:47:32 GMT

Oliver,

Yes, I see,default values:alpha=1, beta=0are more reasonable.

Thanks,

Vladimir