diff --git a/interface/zscal.c b/interface/zscal.c
index 498377343c..0e52d113be 100644
--- a/interface/zscal.c
+++ b/interface/zscal.c
@@ -98,7 +98,7 @@ void CNAME(blasint n, FLOAT alpha_r, void *vx, blasint incx){
   if (nthreads == 1) {
 #endif
 
-  SCAL_K(n, 0, 0, alpha[0], alpha[1], x, incx, NULL, 0, NULL, 0);
+  SCAL_K(n, 0, 0, alpha[0], alpha[1], x, incx, NULL, 0, NULL, 1);
 
 #ifdef SMP
   } else {
@@ -108,7 +108,7 @@ void CNAME(blasint n, FLOAT alpha_r, void *vx, blasint incx){
     mode  =  BLAS_SINGLE | BLAS_COMPLEX;
 #endif
 
-    blas_level1_thread(mode, n, 0, 0,  alpha, x, incx, NULL, 0, NULL, 0, (int (*)(void))SCAL_K, nthreads);
+    blas_level1_thread(mode, n, 0, 0,  alpha, x, incx, NULL, 0, NULL, 1, (int (*)(void))SCAL_K, nthreads);
 
   }
 #endif
diff --git a/kernel/zarch/cscal.c b/kernel/zarch/cscal.c
index e623f306b9..1c9f2cda73 100644
--- a/kernel/zarch/cscal.c
+++ b/kernel/zarch/cscal.c
@@ -210,7 +210,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
       BLASLONG n1 = n & -2;
 
       if (da_i == 0.0) {
-
+       if (dummy2 == 0) {
         while (j < n1) {
 
           x[i] = 0.0;
@@ -230,11 +230,43 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
           j++;
 
         }
+       } else {
+	  while (j < n1) {
+	  if (isnan(x[i]) || isinf(x[i]) || isnan(x[i+1])) {
+	  x[i] = NAN;
+	  x[i+1] = NAN;
+	  }else{
+          x[i] = 0.0;
+          x[i + 1] = 0.0;
+	  }
+	  if (isnan(x[i+inc_x]) || isinf(x[i+inc_x]) || isnan(x[i+1+inc_x])) {
+          x[i + inc_x] = NAN;
+          x[i + 1 + inc_x] = NAN;
+	  } else {
+          x[i + inc_x] = 0.0;
+          x[i + 1 + inc_x] = 0.0;
+	  }
+          i += 2 * inc_x;
+          j += 2;
+
+        }
 
+        while (j < n) {
+	  if (isnan(x[i]) || isinf(x[i]) || isnan(x[i+1])) {
+	  x[i] = NAN;
+	  x[i+1] = NAN;
+	  }else{
+          x[i] = 0.0;
+          x[i + 1] = 0.0;
+	  }
+          i += inc_x;
+          j++;
+	 }
+        }
       } else {
 
         while (j < n1) {
-          if (isnan(x[i]) || isinf(x[i]))
+          if (isnan(x[i]) || isinf(x[i])) 
             temp0   = NAN;
 	  else
 	    temp0 = -da_i * x[i + 1];
@@ -276,7 +308,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 
     } else {
 
-      if (da_i == 0.0) {
+      if (da_i == 0.0 && dummy2) {
         BLASLONG n1 = n & -2;
 
         while (j < n1) {
@@ -335,12 +367,16 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
     alpha[1] = da_i;
 
     if (da_r == 0.0)
-      if (da_i == 0)
+      if (da_i == 0 && dummy2 == 0)
         cscal_kernel_16_zero(n1, x);
-      else
+      else {
+/*	if (dummy2 == 0)
         cscal_kernel_16_zero_r(n1, alpha, x);
-    else if (da_i == 0)
-      cscal_kernel_16_zero_i(n1, alpha, x);
+	else*/
+	cscal_kernel_16(n1, da_r, da_i, x);
+      }
+/*    else if (da_i == 0 && !isnan(da_r))
+      cscal_kernel_16/*_zero_i(n1, alpha, x);*/
     else
       cscal_kernel_16(n1, da_r, da_i, x);
 
@@ -354,7 +390,8 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
       float res = 0.0;
       if (isnan(da_r)) res = da_r;
       while (j < n) {
-
+	if (dummy2)
+		if (isnan(x[i])|| isnan(x[i+1])) res=NAN;
         x[i] = res;
         x[i + 1] = res;
         i += 2;
@@ -382,7 +419,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
           x[i + 1] = da_i * x[i];
 	else
 	  x[i + 1] = NAN;
-	if (x[i] == x[i])
+	if (!isnan(x[i]))
           x[i] = temp0;
         i += 2;
         j++;
@@ -398,7 +435,18 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
       while (j < n) {
 
         temp0 = da_r * x[i];
-        x[i + 1] = da_r * x[i + 1];
+	if (dummy2) {
+	if (isnan(x[i])||isinf(x[i]))temp0 = NAN;
+	if (isnan(x[i+1])||isinf(x[i+1])) 
+		x[i+1] = NAN;
+	else
+		x[i+1] = da_r * x[i + 1];
+	} else {
+		if (isnan(x[i]))
+			x[i + 1] = NAN;
+		else
+        		x[i + 1] = da_r * x[i + 1];
+	}
         x[i] = temp0;
         i += 2;
         j++;
@@ -411,7 +459,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 
         temp0 = da_r * x[i] - da_i * x[i + 1];
         x[i + 1] = da_r * x[i + 1] + da_i * x[i];
-        x[i] = temp0;
+        if (!isnan(x[i])) x[i] = temp0;
         i += 2;
         j++;
 
diff --git a/kernel/zarch/zscal.c b/kernel/zarch/zscal.c
index 36466a6e0a..5111bc4551 100644
--- a/kernel/zarch/zscal.c
+++ b/kernel/zarch/zscal.c
@@ -208,7 +208,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
       BLASLONG n1 = n & -2;
 
       if (da_i == 0.0) {
-
+       if (dummy2 == 0) {
         while (j < n1) {
 
           x[i] = 0.0;
@@ -228,7 +228,38 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
           j++;
 
         }
-
+       } else {
+        while (j < n1) {
+	if (isnan(x[i]) || isinf(x[i]) || isnan(x[i+1])) {
+		x[i] = NAN;
+		x[i+1] = NAN;
+	} else {
+		x[i] = 0.0;
+		x[i+1] = 0.0;
+	}
+	if (isnan(x[i+inc_x]) || isinf(x[i+inc_x]) || isnan(x[i+inc_x+1])) {
+		x[i + inc_x] = NAN;
+		x[i + inc_x + 1] = NAN;
+	} else {
+		x[i + inc_x] = 0.;
+		x[i + inc_x + 1] = 0.;
+	}
+	i += 2 * inc_x;
+	j += 2;
+	}
+	       
+	while (j < n) {
+	if (isnan(x[i]) || isinf(x[i]) || isnan(x[i+1])) {
+		x[i] = NAN;
+		x[i+1] = NAN;
+	} else {
+		x[i] = 0.;
+		x[i+1] = 0.;
+	}
+		i += inc_x;
+		j++;
+        }
+       }
       } else {
 
         while (j < n1) {
@@ -276,7 +307,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 
     } else {
 
-      if (da_i == 0.0) {
+      if (da_i == 0.0 && dummy2) {
         BLASLONG n1 = n & -2;
 
         while (j < n1) {
@@ -335,12 +366,10 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
     alpha[1] = da_i;
 
     if (da_r == 0.0)
-      if (da_i == 0)
+      if (da_i == 0 && dummy2 == 0)
         zscal_kernel_8_zero(n1, x);
       else
         zscal_kernel_8(n1, da_r, da_i, x);
-    else if (da_i == 0 && da_r == da_r)
-      zscal_kernel_8_zero_i(n1, alpha, x);
     else
       zscal_kernel_8(n1, da_r, da_i, x);
 
@@ -354,7 +383,8 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
       double res= 0.0;
       if (isnan(da_r)) res = da_r;
       while (j < n) {
-
+	if (dummy2)
+		if (isnan(x[i]) || isnan(x[i+1])) res = NAN;
         x[i] = res;
         x[i + 1] = res;
         i += 2;
@@ -381,7 +411,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
         	x[i + 1] = da_i * x[i];
 	else
 		x[i + 1] = NAN;
-	if (x[i]==x[i])
+	if (!isnan(x[i]))
           x[i] = temp0;
         i += 2;
         j++;
@@ -397,8 +427,19 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
       while (j < n) {
 
         temp0 = da_r * x[i];
-        x[i + 1] = da_r * x[i + 1];
-        x[i] = temp0;
+	if (dummy2) {
+		if (isnan(x[i]) || isinf(x[i])) temp0 = NAN;
+		if (isnan(x[i + 1]) || isinf(x[i + 1]))
+			x[i + 1] = NAN;
+		else
+			x[i + 1] = da_r * x[i + 1];
+	} else {
+		if (isnan(x[i]))
+			x[i + 1] = NAN;
+		else
+	        x[i + 1] = da_r * x[i + 1];
+	}
+	x[i] = temp0;
         i += 2;
         j++;
 
@@ -410,7 +451,7 @@ int CNAME(BLASLONG n, BLASLONG dummy0, BLASLONG dummy1, FLOAT da_r, FLOAT da_i,
 
         temp0 = da_r * x[i] - da_i * x[i + 1];
         x[i + 1] = da_r * x[i + 1] + da_i * x[i];
-        x[i] = temp0;
+        if (!isnan(x[i])) x[i] = temp0;
         i += 2;
         j++;
 
diff --git a/utest/test_gemv.c b/utest/test_gemv.c
index dab6d2f11f..a0ac3e69e1 100644
--- a/utest/test_gemv.c
+++ b/utest/test_gemv.c
@@ -128,3 +128,487 @@ CTEST(dgemv, 0_nan_inf_incy_2)
 }
 
 #endif
+
+#ifdef BUILD_COMPLEX
+
+CTEST(cgemv, 0_nan_inf)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 1;
+    float alpha[2] = {0.0, 0.0};
+    float beta[2] = {0.0, 0.0};
+    char  trans = 'N';
+    float A[17 * 17 * 4];
+    float X[17 * 2];
+    float Y[17 * 2];
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        Y[i]     = NAN;
+        Y[i + 1] = NAN;
+
+        Y[i + 2] = INFINITY;
+        Y[i + 3] = INFINITY;
+    }
+    Y[2 * N - 1] = NAN;
+    Y[2 * N - 2] = NAN;
+    BLASFUNC(cgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 2 * N; i ++)fprintf(stderr,"should be zero %g\n",Y[i]);
+    for (i = 0; i < 2 * N; i ++)
+        ASSERT_TRUE(Y[i] == 0.0);
+}
+
+CTEST(cgemv, 0_nan_inf_incy_2)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 2;
+    float alpha[2] = {0.0, 0.0};
+    float beta[2] = {0.0, 0.0};
+    char  trans = 'N';
+    float A[17 * 17 * 4];
+    float X[17];
+    float Y[17 * 2 * 2];
+    float *ay = Y;
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    memset(Y, 0, sizeof(Y));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        ay[0]   = NAN;
+        ay[1]   = NAN;
+        ay     += 4;
+        ay[0]   = INFINITY;
+        ay[1]   = INFINITY;
+        ay     += 4;
+    }
+    Y[4 * N - 4] = NAN;
+    Y[4 * N - 3] = NAN;
+    BLASFUNC(cgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 4 * N; i ++)
+        ASSERT_TRUE(Y[i] == 0.0);
+}
+
+CTEST(cgemv, 0_2_nan_1_inf_1)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 1;
+    float alpha[2] = {0.0, 0.0};
+    float beta[2] = {0.0, 2.0};
+    char  trans = 'N';
+    float A[17 * 17 * 4];
+    float X[17 * 2];
+    float Y[17 * 2];
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        Y[i]     = NAN;
+        Y[i + 1] = 1.0;
+
+        Y[i + 2] = INFINITY;
+        Y[i + 3] = 1.0;
+    }
+    Y[2 * N - 2] = NAN;
+    Y[2 * N - 1] = 1.0;
+    BLASFUNC(cgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 2 * N; i += 2) {
+    fprintf(stderr,"should be NAN %g\n",Y[i]);
+    if ((i>>1)%2)fprintf(stderr,"should be INF %g\n",Y[i+1]);
+    else fprintf(stderr,"should be NAN too %g\n",Y[i+1]);
+        if ((i >> 1) % 2){
+            ASSERT_TRUE(isnan(Y[i]));
+            ASSERT_TRUE(isinf(Y[i + 1]));
+        }
+        else {
+            ASSERT_TRUE(isnan(Y[i]));
+            ASSERT_TRUE(isnan(Y[i + 1]));
+        }
+    }
+}
+
+CTEST(cgemv, 0_2_nan_1_inf_1_incy_2)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 2;
+    float alpha[2] = {0.0, 0.0};
+    float beta[2] = {0.0, 2.0};
+    char  trans = 'N';
+    float A[17 * 17 * 4];
+    float X[17];
+    float Y[17 * 2 * 2];
+    float *ay = Y;
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    memset(Y, 0, sizeof(Y));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        ay[0]   = NAN;
+        ay[1]   = 1.0;
+        ay     += 4;
+        ay[0]   = INFINITY;
+        ay[1]   = 1.0;
+        ay     += 4;
+    }
+    Y[4 * N - 4] = NAN;
+    Y[4 * N - 3] = 1.0;
+    BLASFUNC(cgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 4 * N; i += 2) {
+        if ((i >> 1) % 2) {
+            ASSERT_TRUE(Y[i] == 0.0);
+            ASSERT_TRUE(Y[i + 1] == 0.0);
+        }
+        else {
+            if ((i >> 2) % 2) {
+                ASSERT_TRUE(isnan(Y[i]));
+                ASSERT_TRUE(isinf(Y[i + 1]));
+            }
+            else {
+                ASSERT_TRUE(isnan(Y[i]));
+                ASSERT_TRUE(isnan(Y[i + 1]));
+            }
+        }
+    }
+}
+
+CTEST(cgemv, 2_0_nan_1_inf_1)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 1;
+    float alpha[2] = {0.0, 0.0};
+    float beta[2] = {2.0, 0.0};
+    char  trans = 'N';
+    float A[17 * 17 * 4];
+    float X[17 * 2];
+    float Y[17 * 2];
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        Y[i]     = NAN;
+        Y[i + 1] = 1.0;
+
+        Y[i + 2] = INFINITY;
+        Y[i + 3] = 1.0;
+    }
+    Y[2 * N - 2] = NAN;
+    Y[2 * N - 1] = 1.0;
+    BLASFUNC(cgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 2 * N; i += 2) {
+    if ((i>>1)%2)fprintf(stderr,"Yi should be INF %g\n",Y[i]);
+    else fprintf(stderr,"Yi should be NAN %g\n",Y[i]);
+    fprintf(stderr,"Yi+1 should be NAN %g\n",Y[i+1]);
+        if ((i >> 1) % 2){
+            ASSERT_TRUE(isinf(Y[i]));
+            ASSERT_TRUE(isnan(Y[i + 1]));
+        }
+        else {
+            ASSERT_TRUE(isnan(Y[i]));
+            ASSERT_TRUE(isnan(Y[i + 1]));
+        }
+    }
+}
+
+CTEST(cgemv, 2_0_nan_1_inf_1_incy_2)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 2;
+    float alpha[2] = {0.0, 0.0};
+    float beta[2] = {2.0, 0.0};
+    char  trans = 'N';
+    float A[17 * 17 * 4];
+    float X[17];
+    float Y[17 * 2 * 2];
+    float *ay = Y;
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    memset(Y, 0, sizeof(Y));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        ay[0]   = NAN;
+        ay[1]   = 1.0;
+        ay     += 4;
+        ay[0]   = INFINITY;
+        ay[1]   = 1.0;
+        ay     += 4;
+    }
+    Y[4 * N - 4] = NAN;
+    Y[4 * N - 3] = 1.0;
+    BLASFUNC(cgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 4 * N; i += 2) {
+        if ((i >> 1) % 2) {
+            ASSERT_TRUE(Y[i] == 0.0);
+            ASSERT_TRUE(Y[i + 1] == 0.0);
+        }
+        else {
+            if ((i >> 2) % 2) {
+                ASSERT_TRUE(isinf(Y[i]));
+                ASSERT_TRUE(isnan(Y[i + 1]));
+            }
+            else {
+                ASSERT_TRUE(isnan(Y[i]));
+                ASSERT_TRUE(isnan(Y[i + 1]));
+            }
+        }
+    }
+}
+
+#endif
+
+#ifdef BUILD_COMPLEX16
+
+CTEST(zgemv, 0_nan_inf)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 1;
+    double alpha[2] = {0.0, 0.0};
+    double beta[2] = {0.0, 0.0};
+    char  trans = 'N';
+    double A[17 * 17 * 4];
+    double X[17 * 2];
+    double Y[17 * 2];
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        Y[i]     = NAN;
+        Y[i + 1] = NAN;
+
+        Y[i + 2] = INFINITY;
+        Y[i + 3] = INFINITY;
+    }
+    Y[2 * N - 1] = NAN;
+    Y[2 * N - 2] = NAN;
+    BLASFUNC(zgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 2 * N; i ++)
+        ASSERT_TRUE(Y[i] == 0.0);
+}
+
+CTEST(zgemv, 0_nan_inf_incy_2)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 2;
+    double alpha[2] = {0.0, 0.0};
+    double beta[2] = {0.0, 0.0};
+    char  trans = 'N';
+    double A[17 * 17 * 4];
+    double X[17];
+    double Y[17 * 2 * 2];
+    double *ay = Y;
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    memset(Y, 0, sizeof(Y));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        ay[0]   = NAN;
+        ay[1]   = NAN;
+        ay     += 4;
+        ay[0]   = INFINITY;
+        ay[1]   = INFINITY;
+        ay     += 4;
+    }
+    Y[4 * N - 4] = NAN;
+    Y[4 * N - 3] = NAN;
+    BLASFUNC(zgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 4 * N; i ++)
+        ASSERT_TRUE(Y[i] == 0.0);
+}
+
+CTEST(zgemv, 0_2_nan_1_inf_1)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 1;
+    double alpha[2] = {0.0, 0.0};
+    double beta[2] = {0.0, 2.0};
+    char  trans = 'N';
+    double A[17 * 17 * 4];
+    double X[17 * 2];
+    double Y[17 * 2];
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        Y[i]     = NAN;
+        Y[i + 1] = 1.0;
+
+        Y[i + 2] = INFINITY;
+        Y[i + 3] = 1.0;
+    }
+    Y[2 * N - 2] = NAN;
+    Y[2 * N - 1] = 1.0;
+    BLASFUNC(zgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 2 * N; i += 2) {
+        if ((i >> 1) % 2){
+            ASSERT_TRUE(isnan(Y[i]));
+            ASSERT_TRUE(isinf(Y[i + 1]));
+        }
+        else {
+            ASSERT_TRUE(isnan(Y[i]));
+            ASSERT_TRUE(isnan(Y[i + 1]));
+        }
+    }
+}
+
+CTEST(zgemv, 0_2_nan_1_inf_1_incy_2)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 2;
+    double alpha[2] = {0.0, 0.0};
+    double beta[2] = {0.0, 2.0};
+    char  trans = 'N';
+    double A[17 * 17 * 4];
+    double X[17];
+    double Y[17 * 2 * 2];
+    double *ay = Y;
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    memset(Y, 0, sizeof(Y));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        ay[0]   = NAN;
+        ay[1]   = 1.0;
+        ay     += 4;
+        ay[0]   = INFINITY;
+        ay[1]   = 1.0;
+        ay     += 4;
+    }
+    Y[4 * N - 4] = NAN;
+    Y[4 * N - 3] = 1.0;
+    BLASFUNC(zgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 4 * N; i += 2) {
+        if ((i >> 1) % 2) {
+            ASSERT_TRUE(Y[i] == 0.0);
+            ASSERT_TRUE(Y[i + 1] == 0.0);
+        }
+        else {
+            if ((i >> 2) % 2) {
+                ASSERT_TRUE(isnan(Y[i]));
+                ASSERT_TRUE(isinf(Y[i + 1]));
+            }
+            else {
+                ASSERT_TRUE(isnan(Y[i]));
+                ASSERT_TRUE(isnan(Y[i + 1]));
+            }
+        }
+    }
+}
+
+CTEST(zgemv, 2_0_nan_1_inf_1)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 1;
+    double alpha[2] = {0.0, 0.0};
+    double beta[2] = {2.0, 0.0};
+    char  trans = 'N';
+    double A[17 * 17 * 4];
+    double X[17 * 2];
+    double Y[17 * 2];
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        Y[i]     = NAN;
+        Y[i + 1] = 1.0;
+
+        Y[i + 2] = INFINITY;
+        Y[i + 3] = 1.0;
+    }
+    Y[2 * N - 2] = NAN;
+    Y[2 * N - 1] = 1.0;
+    BLASFUNC(zgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 2 * N; i += 2) {
+        if ((i >> 1) % 2)fprintf(stderr,"Yi sollte inf sein %g\n",Y[i]);
+	else fprintf(stderr,"Yi sollte nan sein %g\n",Y[i]);
+	fprintf(stderr,"Yi+1 sollte nan sein %g\n",Y[i+1]);
+        if ((i >> 1) % 2){
+            ASSERT_TRUE(isinf(Y[i]));
+            ASSERT_TRUE(isnan(Y[i + 1]));
+        }
+        else {
+            ASSERT_TRUE(isnan(Y[i]));
+            ASSERT_TRUE(isnan(Y[i + 1]));
+        }
+    }
+}
+
+CTEST(zgemv, 2_0_nan_1_inf_1_incy_2)
+{
+    int i;
+    blasint N = 17;
+    blasint incX = 1;
+    blasint incY = 2;
+    double alpha[2] = {0.0, 0.0};
+    double beta[2] = {2.0, 0.0};
+    char  trans = 'N';
+    double A[17 * 17 * 4];
+    double X[17];
+    double Y[17 * 2 * 2];
+    double *ay = Y;
+
+    memset(A, 0, sizeof(A));
+    memset(X, 0, sizeof(X));
+    memset(Y, 0, sizeof(Y));
+    for (i = 0; i < (2 * N - 2); i += 4)
+    {
+        ay[0]   = NAN;
+        ay[1]   = 1.0;
+        ay     += 4;
+        ay[0]   = INFINITY;
+        ay[1]   = 1.0;
+        ay     += 4;
+    }
+    Y[4 * N - 4] = NAN;
+    Y[4 * N - 3] = 1.0;
+    BLASFUNC(zgemv)(&trans, &N, &N, alpha, A, &N, X, &incX, beta, Y, &incY);
+    for (i = 0; i < 4 * N; i += 2) {
+        if ((i >> 1) % 2) {
+            ASSERT_TRUE(Y[i] == 0.0);
+            ASSERT_TRUE(Y[i + 1] == 0.0);
+        }
+        else {
+            if ((i >> 2) % 2) {
+                ASSERT_TRUE(isinf(Y[i]));
+                ASSERT_TRUE(isnan(Y[i + 1]));
+            }
+            else {
+                ASSERT_TRUE(isnan(Y[i]));
+                ASSERT_TRUE(isnan(Y[i + 1]));
+            }
+        }
+    }
+}
+
+#endif
diff --git a/utest/test_zscal.c b/utest/test_zscal.c
index 09e63752c2..57d78b6901 100644
--- a/utest/test_zscal.c
+++ b/utest/test_zscal.c
@@ -442,6 +442,33 @@ CTEST(cscal, i_0inf_inc_2)
     ASSERT_TRUE(isnan(inf[17]));
 }
 
+CTEST(cscal, i00_NAN)
+{
+    blasint N=9;
+    blasint incX=1;
+    float i[] = {0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0 };
+    float nan[] = {NAN, 0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
+    BLASFUNC(cscal)(&N, i, nan, &incX);
+    ASSERT_TRUE(isnan(nan[0]));
+    ASSERT_TRUE(isnan(nan[1]));
+    ASSERT_TRUE(isnan(nan[16]));
+    ASSERT_TRUE(isnan(nan[17]));
+}
+
+CTEST(cscal, i00_NAN_incx_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    float i[] = {0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0 };
+    float nan[] = {0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN,
+                   0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN};
+    BLASFUNC(cscal)(&N, i, nan, &incX);
+    ASSERT_TRUE(isnan(nan[0]));
+    ASSERT_TRUE(isnan(nan[1]));
+    ASSERT_TRUE(isnan(nan[16]));
+    ASSERT_TRUE(isnan(nan[17]));
+}
+
 #endif
 
 #ifdef BUILD_COMPLEX16
@@ -588,4 +615,31 @@ CTEST(zscal, i_0inf_inc_2)
     ASSERT_TRUE(isnan(inf[17]));
 }
 
+CTEST(zscal, i00_NAN)
+{
+    blasint N=9;
+    blasint incX=1;
+    double i[] = {0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0 };
+    double nan[] = {NAN, 0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0, NAN,0};
+    BLASFUNC(zscal)(&N, i, nan, &incX);
+    ASSERT_TRUE(isnan(nan[0]));
+    ASSERT_TRUE(isnan(nan[1]));
+    ASSERT_TRUE(isnan(nan[16]));
+    ASSERT_TRUE(isnan(nan[17]));
+}
+
+CTEST(zscal, i00_NAN_incx_2)
+{
+    blasint N=9;
+    blasint incX=2;
+    double i[] = {0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0 };
+    double nan[] = {0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN,
+                    0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN, 0,NAN};
+    BLASFUNC(zscal)(&N, i, nan, &incX);
+    ASSERT_TRUE(isnan(nan[0]));
+    ASSERT_TRUE(isnan(nan[1]));
+    ASSERT_TRUE(isnan(nan[16]));
+    ASSERT_TRUE(isnan(nan[17]));
+}
+
 #endif