diff --git a/utest/CMakeLists.txt b/utest/CMakeLists.txt
index 4771d8a27..6a61899da 100644
--- a/utest/CMakeLists.txt
+++ b/utest/CMakeLists.txt
@@ -18,6 +18,7 @@ else ()
     test_zscal.c
     test_amin.c
     test_axpby.c
+    test_gemv.c
   )
 endif ()
 
diff --git a/utest/Makefile b/utest/Makefile
index 36acf96cd..b82937093 100644
--- a/utest/Makefile
+++ b/utest/Makefile
@@ -14,7 +14,7 @@ UTESTEXTBIN=openblas_utest_ext
 include $(TOPDIR)/Makefile.system
 
 OBJS=utest_main.o test_min.o test_amax.o test_ismin.o test_rotmg.o test_axpy.o test_dotu.o test_dsdot.o test_swap.o test_rot.o test_dnrm2.o test_zscal.o \
-     test_amin.o test_axpby.o
+     test_amin.o test_axpby.o test_gemv.o
 #test_rot.o test_swap.o test_axpy.o test_dotu.o test_dsdot.o test_fork.o
 OBJS_EXT=utest_main.o $(DIR_EXT)/xerbla.o $(DIR_EXT)/common.o 
 OBJS_EXT+=$(DIR_EXT)/test_isamin.o $(DIR_EXT)/test_idamin.o $(DIR_EXT)/test_icamin.o $(DIR_EXT)/test_izamin.o 
diff --git a/utest/test_gemv.c b/utest/test_gemv.c
new file mode 100644
index 000000000..c85ef3f38
--- /dev/null
+++ b/utest/test_gemv.c
@@ -0,0 +1,126 @@
+#include "openblas_utest.h"
+#include <cblas.h>
+
+#ifndef NAN
+#define NAN 0.0/0.0
+#endif
+#ifndef INFINITY
+#define INFINITY 1.0/0.0
+#endif
+
+#ifdef BUILD_SINGLE
+
+CTEST(sgemv, 0_nan_inf)
+{
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 1;
+    float alpha = 0.0;
+    float beta = 0.0;
+    char  trans = 'N';
+    float A[N * N];
+    float X[N];
+    float Y[N];
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    for (int i = 0; i < (N - 1); i += 2)
+    {
+        Y[i]     = NAN;
+        Y[i + 1] = INFINITY;
+    }
+    Y[N - 1] = NAN;
+    BLASFUNC(sgemv)(&trans, &N, &N, &alpha, A, &N, X, &incX, &beta, Y, &incY);
+    for (int i = 0; i < N; i ++)
+        ASSERT_TRUE(Y[i] == 0.0);
+}
+
+CTEST(sgemv, 0_nan_inf_incy_2)
+{
+    blasint N  = 17;
+    blasint Ny = 33;
+    blasint incX = 1;
+    blasint incY = 2;
+    float alpha = 0.0;
+    float beta = 0.0;
+    char  trans = 'N';
+    float A[N * N];
+    float X[N];
+    float Y[Ny];
+    float *ay = Y;
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    memset(Y, 0, sizeof(Y));
+    for (int i = 0; i < (N - 1); i += 2)
+    {
+        ay[0]   = NAN;
+        ay     += 2;
+        ay[0]   = INFINITY;
+        ay     += 2;
+    }
+    Y[Ny - 1] = NAN;
+    BLASFUNC(sgemv)(&trans, &N, &N, &alpha, A, &N, X, &incX, &beta, Y, &incY);
+    for (int i = 0; i < Ny; i ++)
+        ASSERT_TRUE(Y[i] == 0.0);
+}
+
+#endif
+
+#ifdef BUILD_DOUBLE
+CTEST(dgemv, 0_nan_inf)
+{
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 1;
+    double alpha = 0.0;
+    double beta = 0.0;
+    char  trans = 'N';
+    double A[N * N];
+    double X[N];
+    double Y[N];
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    for (int i = 0; i < (N - 1); i += 2)
+    {
+        Y[i]     = NAN;
+        Y[i + 1] = INFINITY;
+    }
+    Y[N - 1] = NAN;
+    BLASFUNC(dgemv)(&trans, &N, &N, &alpha, A, &N, X, &incX, &beta, Y, &incY);
+    for (int i = 0; i < N; i ++)
+        ASSERT_TRUE(Y[i] == 0.0);
+}
+
+CTEST(dgemv, 0_nan_inf_incy_2)
+{
+    blasint N  = 17;
+    blasint Ny = 33;
+    blasint incX = 1;
+    blasint incY = 2;
+    double alpha = 0.0;
+    double beta = 0.0;
+    char  trans = 'N';
+    double A[N * N];
+    double X[N];
+    double Y[Ny];
+    double *ay = Y;
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    memset(Y, 0, sizeof(Y));
+    for (int i = 0; i < (N - 1); i += 2)
+    {
+        ay[0]   = NAN;
+        ay     += 2;
+        ay[0]   = INFINITY;
+        ay     += 2;
+    }
+    Y[Ny - 1] = NAN;
+    BLASFUNC(dgemv)(&trans, &N, &N, &alpha, A, &N, X, &incX, &beta, Y, &incY);
+    for (int i = 0; i < Ny; i ++)
+        ASSERT_TRUE(Y[i] == 0.0);
+}
+
+#endif