ZGEMV / Blas returns NAN

[bose]

Hi,

I have been working on some simple code to determine some eigenvalues and I needed to compute a matrix vector product for an Arnoldi iteration (using ARPACK). I have manually verified (by parsing each element in large sparse matrices ) that elements in both the matrix and vector of order ~1000, are initialized, finite, and most importantly not NAN.

I make a call to zgemv as follows :

CALL ZGEMV('N', n, n, DCMPLX(1d0, 0d0), A, n,
& t1, 1, DCMPLX(0d0, 0d0), t2, 1)


As mentioned before, A, t1, and t2 are all properly initialized. But on return, t2 is a vector of NAN's. I can post more code if that would be helpful. Does anyone know what the trouble might be?

Thank you in advance for all your help ,

Sanjeeb T. Bose

[Julien Langou]

At a glance your call looks correct to me.

Which BLAS are you using?

You can try to use another BLAS to give a change, here is one:
http://www.netlib.org/blas/zgemv.f
this is the reference BLAS implementation. This one is correct for sure.

-julien

[bose]

I have installed BLAS from an ubuntu distribution; the shared library I linke at runtime is libblas.so.3. Although, I have since programmed a small test function just to try ZGEMV and it returns an answer that is not NAN.

Do you know of anything that would cause ZGEMV to return NAN? Thanks for all your help,

Sanjeeb

[Julien Langou]

Do you know of anything that would cause ZGEMV to return NAN? Thanks

A buggy ZGEMV or a buggy call.
Otherwise no there is no reason that comes to my mind.
Can you try to ZGEMV I gave you?
This will have the same interface as the one from the Ubuntu distribution.
Julien

[bose]

I tried the call with your zgemv from the netlib repository. Unfortunately, that did not seem to work. Would it be possible for me to send you the code and have you take a quick glance at it? I understand that this is a large inconvience for you -- thank you so much for all your help.

Thanks again,
Sanjeeb.

[Julien Langou]

Although looking at a simple line of code is not enough, I'd rather have a small main that you can post on this forum to look at. Please try to clean around and to restrict your codes to basically initialize a matrix (could be random I guess), call zgemv, then check for you NaN. If you have NaN then I'll have a look. ok?
-j

[bose]

I have played around with the code and have noticed something strange. I use ZGEMV early in the code (for testing purposes) and it returns an answer : (1000154, 0). Some other lines are executed that do not involve any of A, t1, or t2. Upon recomputation, ZGEMV returns (NAN, NAN), and I'm not sure why.

I am posting the relevant lines of code :

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
DO i = 1,n
t1(i) = DCMPLX(i)
t2(i) = DCMPLX(0d0, 0d0)
ENDDO

CALL ZGEMV('N', n, n, DCMPLX(1d0, 0d0), A, n,
& t1, 1, DCMPLX(0d0, 0d0), t2, 1)

WRITE(*, *) t2(1)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

c Compute the LU factorization of the banded matrix B
c B has 2 sub- and 2 super- diagonals

CALL ZGBTRF(n, n, kl, ku, BAND, (2*kl + ku + 1), ipiv, info)


IF (info .ne. 0) THEN

WRITE(*, *) 'LU factorization failed ...'
END IF

WRITE(*, *) 'Computing eigenvalues...'



c Compute selected eigenvalues / eigenvectors using
c appropriate krylov subspace method (Arnoldi method)
c
c Call ARPACK fortran subroutines, following ARPACK calls
c below models function call in the ARPACK documentation

c initialize call of znaupd, arpack doc p. 18

ido = 0
iparam(1) = 1
iparam(3) = 2000
iparam(7) = 3
info = 0
iter = 0


10 CONTINUE

iter = iter + 1
WRITE(*, *) 'Iteration = ', iter, '...'


CALL ZNAUPD(ido, 'G', n, 'LM', nev, tol, resid,
& ncv, v, ldv, iparam, ipntr, workd, workl, lworkl,
& rwork, info )

IF (ido .eq. -1 .OR. ido .eq. 1) THEN


c First, compute the product z <- Az

DO i = 1, n

c t1(i) = DCMPLX(workd(ipntr(1) + i - 1))
t1(i) = DCMPLX(i)
t2(i) = DCMPLX(0d0, 0d0)
ENDDO


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

CALL ZGEMV('N', n, n, DCMPLX(1d0, 0d0), A, n,
& t1, 1, DCMPLX(0d0, 0d0), t2, 1)

c CALL ZGBMV('N', n, n, kl, ku, DCMPLX(1d0, 0d0), A_BAND,
c & kl + ku + 1, t1, 1, DCMPLX(0d0, 0d0), t2, 1)

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

WRITE(*, *) t2(1)


The calls to ZGEMV are contained within comment lines.