diff --git a/paddle/cuda/include/hl_cpu_matrix_kernel.cuh b/paddle/cuda/include/hl_cpu_matrix_kernel.cuh
index f35bfbc5c8253d632f8089f5037421f527633aad..9c49a4bd2083794e98b099b25944bedec3d5a2ff 100644
--- a/paddle/cuda/include/hl_cpu_matrix_kernel.cuh
+++ b/paddle/cuda/include/hl_cpu_matrix_kernel.cuh
@@ -17,7 +17,11 @@ limitations under the License. */
 
 #include <stdio.h>
 #include "hl_base.h"
+#if defined(__ARM_NEON__) || defined(__ARM_NEON)
+#include "hl_neon_matrix_kernel.cuh"
+#else
 #include "hl_sse_matrix_kernel.cuh"
+#endif
 
 /**
  * @brief   cpu element wise unary operator.
diff --git a/paddle/cuda/include/hl_matrix_base.cuh b/paddle/cuda/include/hl_matrix_base.cuh
index db35ee2037433163ebb3673edb350e3fab71fba9..8b755c1095c2c4fdb7e74d8cddc948e6a6af380b 100644
--- a/paddle/cuda/include/hl_matrix_base.cuh
+++ b/paddle/cuda/include/hl_matrix_base.cuh
@@ -66,6 +66,8 @@ typedef BaseOp SSESquaredDiff;
 typedef BaseOp SSEFirst;
 typedef BaseOp SSESecond;
 typedef BaseOp SSEClassificationError;
+#elif defined(__ARM__NEON__) || defined(__ARM_NEON)
+#include "hl_matrix_base_neon.cuh"
 #else
 #include "hl_matrix_base_sse.cuh"
 #endif
diff --git a/paddle/cuda/include/hl_matrix_base_neon.cuh b/paddle/cuda/include/hl_matrix_base_neon.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..e13019f5ee24ad600005c99678426ee3808b0e54
--- /dev/null
+++ b/paddle/cuda/include/hl_matrix_base_neon.cuh
@@ -0,0 +1,161 @@
+/* Copyright (c) 2016 PaddlePaddle Authors. All Rights Reserve.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License. */
+
+
+#ifndef HL_MATRIX_BASE_NEON_CUH_
+#define HL_MATRIX_BASE_NEON_CUH_
+
+namespace aggregate {
+class SSESum {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    return vaddq_f32(a, b);
+  }
+};
+
+class SSEMax {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    return vmaxq_f32(a, b);
+  }
+};
+
+class SSEMin {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    return vminq_f32(a, b);
+  }
+};
+}  // namespace aggregate
+
+namespace base {
+namespace unary {
+class SSEIdentity {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a) const {
+    return a;
+  }
+};
+}  // namespace unary
+
+namespace binary {
+class SSEAdd {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    return vaddq_f32(a, b);
+  }
+};
+
+class SSEAdd2 {
+public:
+  static const bool sse = true;
+  const real p1;
+  const real p2;
+  float32x4_t mp1;
+  float32x4_t mp2;
+
+public:
+  SSEAdd2(const real s1, const real s2) : p1(s1), p2(s2) {
+    mp1 = vdupq_n_f32(p1);
+    mp2 = vdupq_n_f32(p2);
+  }
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    float32x4_t tmp1, tmp2;
+    tmp1 = vmulq_f32(mp1, a);
+    tmp2 = vmulq_f32(mp2, b);
+    return vaddq_f32(tmp1, tmp2);
+  }
+};
+
+class SSESub {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    return vsubq_f32(a, b);
+  }
+};
+
+class SSEMul {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    return vmulq_f32(a, b);
+  }
+};
+
+class SSEDiv {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    float32x4_t tmp;
+    tmp = vrecpeq_f32(b);
+    return vmulq_f32(a, tmp);
+  }
+};
+
+class SSESquaredDiff {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    float32x4_t tmp;
+    tmp = vsubq_f32(a, b);
+    return vmulq_f32(tmp, tmp);
+  }
+};
+
+class SSEFirst {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    return a;
+  }
+};
+
+class SSESecond {
+public:
+  static const bool sse = true;
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    return b;
+  }
+};
+
+class SSEClassificationError {
+public:
+  static const bool sse = true;
+  const real p;
+  float32x4_t mp;
+  uint32x4_t result;
+
+public:
+  explicit SSEClassificationError(const real s) : p(s) {
+    mp = vdupq_n_f32(p);
+    result = vdupq_n_u32(1);
+  }
+  // TODO: to be check
+  INLINE float32x4_t vecOp(const float32x4_t a, const float32x4_t b) const {
+    uint32x4_t tmp1 = vcgtq_f32(a, mp);
+    uint32x4_t tmp2 = vcgtq_f32(b, mp);
+    uint32x4_t tmp3 = veorq_u32(tmp1, tmp2);
+    return vcvtq_f32_u32(vandq_u32(tmp3, result));
+  }
+};
+}  // namespace binary
+}  // namespace base
+
+#endif /* HL_MATRIX_BASE_NEON_CUH_ */
diff --git a/paddle/cuda/include/hl_neon_matrix_kernel.cuh b/paddle/cuda/include/hl_neon_matrix_kernel.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..7b4e5b00079b66d0a46a1344a43f41962cf50f10
--- /dev/null
+++ b/paddle/cuda/include/hl_neon_matrix_kernel.cuh
@@ -0,0 +1,299 @@
+/* Copyright (c) 2016 PaddlePaddle Authors. All Rights Reserve.
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License. */
+
+
+#ifndef HL_NEON_MATRIX_KERNEL_CUH_
+#define HL_NEON_MATRIX_KERNEL_CUH_
+
+#include "hl_matrix_type.cuh"
+
+#define VECTOR_SIZE     16
+
+/* number of float in vector */
+#define     VECTOR_LEN      4
+#define     VECTOR_SET      vdupq_n_f32
+
+inline bool hl_check_align(size_t size) {
+  return !(size & (VECTOR_SIZE - 1));
+}
+
+inline bool hl_check_align(void *ptr) {
+  return hl_check_align(reinterpret_cast<size_t>(ptr));
+}
+
+template <class Agg>
+inline real hl_agg_op(Agg agg, vecType mm) {
+  float32x4_t rev = vrev64q_f32(mm);
+  float32x4_t tmp1 = agg.vecOp(rev, rev);
+  float32x2_t lo = vget_high_f32(rev);
+  float32x2_t hi = vget_low_f32(rev);
+  float32x4_t tmp2 = vcombine_f32(hi, lo);
+  float32x4_t ret = agg.vecOp(tmp1, tmp2);
+
+  return vgetq_lane_f32(ret, 0);
+}
+
+template <class Agg, class Op, class Saver>
+void hl_sse_matrix_row_op(Agg agg, Op op, Saver sv,
+                          int dimM, int dimN,
+                          real *dst, int ld,
+                          real *A, int lda) {
+  for (int i = 0; i < dimM; i++, A += lda) {
+    vecType mm = VECTOR_SET(agg.init());
+    vecType *a = (vecType*)(A);
+    for (int j = 0; j < dimN / VECTOR_LEN; j++, a++) {
+      mm = agg.vecOp(mm, op.vecOp(*a));
+    }
+
+    int rem = dimN % VECTOR_LEN;
+    if (rem) {
+      real tmp = hl_agg_op(agg, mm);
+      real *a = A + (dimN / VECTOR_LEN) * VECTOR_LEN;
+      for (int j = 0; j < rem; j++) {
+        tmp = agg(tmp, op(a[j]));
+      }
+      dst[i*ld] = sv(dst[i*ld], tmp);
+    } else {
+      dst[i*ld] = sv(dst[i*ld], hl_agg_op(agg, mm));
+    }
+  }
+}
+
+template <class Agg, class Op, class Saver>
+void hl_sse_matrix_row_op(Agg agg, Op op, Saver sv,
+                          int dimM, int dimN,
+                          real *dst, int ld,
+                          real *A, int lda,
+                          real *B, int ldb) {
+  for (int i = 0; i < dimM; i++, A += lda, B += ldb) {
+    vecType mm = VECTOR_SET(agg.init());
+    vecType *a = (vecType*)(A);
+    vecType *b = (vecType*)(B);
+    for (int j = 0; j < dimN / VECTOR_LEN; j++, a++, b++) {
+        mm = agg.vecOp(mm, op.vecOp(*a, *b));
+    }
+
+    int rem = dimN % VECTOR_LEN;
+    if (rem) {
+      real tmp = hl_agg_op(agg, mm);
+      real *a = A + (dimN / VECTOR_LEN) * VECTOR_LEN;
+      real *b = B + (dimN / VECTOR_LEN) * VECTOR_LEN;
+      for (int j = 0; j < rem; j++) {
+          tmp = agg(tmp, op(a[j], b[j]));
+      }
+      dst[i*ld] = sv(dst[i*ld], tmp);
+    } else {
+        dst[i*ld] = sv(dst[i*ld], hl_agg_op(agg, mm));
+    }
+  }
+}
+
+template <class Agg, class Op, class Saver>
+void hl_matrix_column_op(Agg agg, Op op, Saver sv,
+                         int dimM, int dimN,
+                         real *dst,
+                         real *A, int lda) {
+  for (int j = 0; j < dimN; j++) {
+    real tmp = agg.init();
+    for (int i = 0; i < dimM; i++) {
+        tmp = agg(tmp, op(A[i * lda + j]));
+    }
+    dst[j] = sv(dst[j], tmp);
+  }
+}
+
+template <class Agg, class Op, class Saver>
+void hl_matrix_column_op(Agg agg, Op op, Saver sv,
+                         int dimM, int dimN,
+                         real *dst,
+                         real *A, int lda,
+                         real *B, int ldb) {
+  for (int j = 0; j < dimN; j++) {
+    real tmp = agg.init();
+    for (int i = 0; i < dimM; i++) {
+        tmp = agg(tmp, op(A[i * lda + j], B[i * ldb + j]));
+    }
+    dst[j] = sv(dst[j], tmp);
+  }
+}
+
+/*
+ * MaxRow greater than or equal dimN
+ * dimN is multiples of VECTOR_LEN
+ * so rem <= MaxRow / VECTOR_LEN
+ */
+template <int MaxRow, class Agg, class Op, class Saver>
+void hl_sse_column_op_with_rem(Agg agg, Op op, Saver sv,
+                               int dimM, int dimN,
+                               real *dst,
+                               real *A, int lda) {
+  vecType mm[MaxRow / VECTOR_LEN];
+  for (int n = 0; n < MaxRow / VECTOR_LEN; n++) {
+    mm[n] = VECTOR_SET(agg.init());
+  }
+
+  for (int i = 0; i < dimM; i++) {
+    vecType *a = (vecType*)(A + i * lda);
+    for (int n = 0; n < dimN / VECTOR_LEN; n++) {
+      mm[n] = agg.vecOp(mm[n], op.vecOp(a[n]));
+    }
+  }
+
+  vecType *result = (vecType*)(dst);
+  for (int n = 0; n < dimN / VECTOR_LEN; n++) {
+    result[n] = sv.vecOp(result[n], mm[n]);
+  }
+
+  int rem = dimN % VECTOR_LEN;
+  if (rem) {
+    A += (dimN / VECTOR_LEN) * VECTOR_LEN;
+    dst += (dimN / VECTOR_LEN) * VECTOR_LEN;
+    hl_matrix_column_op(agg, op, sv, dimM, rem, dst, A, lda);
+  }
+}
+
+/*
+ * dimN is multiples of VECTOR_LEN
+ * dimN greater than Step
+ */
+template <int Step, class Agg, class Op, class Saver>
+void hl_sse_matrix_column_op(Agg agg, Op op, Saver sv,
+                             int dimM, int dimN,
+                             real *dst,
+                             real *A, int lda) {
+  for (int j = 0; j < dimN / Step; j++, dst += Step, A += Step) {
+    vecType mm[Step / VECTOR_LEN];
+    for (int n = 0; n < Step / VECTOR_LEN; n++) {
+      mm[n] = VECTOR_SET(agg.init());
+    }
+
+    for (int i = 0; i < dimM; i++) {
+      vecType *a = (vecType*)(A + i * lda);
+      for (int n = 0; n < Step / VECTOR_LEN; n++) {
+        mm[n] = agg.vecOp(mm[n], op.vecOp(a[n]));
+      }
+    }
+
+    vecType *result = (vecType*)(dst);
+    for (int n = 0; n < Step / VECTOR_LEN; n++) {
+      result[n] = sv.vecOp(result[n], mm[n]);
+    }
+  }
+
+  int remRow = dimN % Step;
+  if (remRow) {
+    hl_sse_column_op_with_rem<Step>(agg, op, sv, dimM, remRow, dst, A, lda);
+  }
+}
+
+template <class Agg, class Op, class Saver>
+void hl_sse_matrix_column_op(Agg agg, Op op, Saver sv,
+                             int dimM, int dimN,
+                             real *dst,
+                             real *A, int lda) {
+  if (dimN <= 16) {
+    hl_sse_matrix_column_op<16>(agg, op, sv, dimM, dimN, dst, A, lda);
+  } else if (dimN <= 32) {
+    hl_sse_matrix_column_op<32>(agg, op, sv, dimM, dimN, dst, A, lda);
+  } else if (dimN <= 1024 || dimM <= 512) {
+    hl_sse_matrix_column_op<64>(agg, op, sv, dimM, dimN, dst, A, lda);
+  } else {
+    hl_sse_matrix_column_op<1024>(agg, op, sv, dimM, dimN, dst, A, lda);
+  }
+}
+
+template <int MaxRow, class Agg, class Op, class Saver>
+void hl_sse_column_op_with_rem(Agg agg, Op op, Saver sv,
+                               int dimM, int dimN,
+                               real *dst,
+                               real *A, int lda,
+                               real *B, int ldb) {
+  vecType mm[MaxRow / VECTOR_LEN];
+  for (int n = 0; n < MaxRow / VECTOR_LEN; n++) {
+    mm[n] = VECTOR_SET(agg.init());
+  }
+
+  for (int i = 0; i < dimM; i++) {
+    vecType *a = (vecType*)(A + i * lda);
+    vecType *b = (vecType*)(B + i * ldb);
+    for (int n = 0; n < dimN / VECTOR_LEN; n++) {
+      mm[n] = agg.vecOp(mm[n], op.vecOp(a[n], b[n]));
+    }
+  }
+
+  vecType *result = (vecType*)(dst);
+  for (int n = 0; n < dimN / VECTOR_LEN; n++) {
+    result[n] = sv.vecOp(result[n], mm[n]);
+  }
+
+  int rem = dimN % VECTOR_LEN;
+  if (rem) {
+    A += (dimN / VECTOR_LEN) * VECTOR_LEN;
+    B += (dimN / VECTOR_LEN) * VECTOR_LEN;
+    dst += (dimN / VECTOR_LEN) * VECTOR_LEN;
+    hl_matrix_column_op(agg, op, sv, dimM, rem, dst, A, lda, B, ldb);
+  }
+}
+
+template <int Step, class Agg, class Op, class Saver>
+void hl_sse_matrix_column_op(Agg agg, Op op, Saver sv,
+                             int dimM, int dimN,
+                             real *dst,
+                             real *A, int lda,
+                             real *B, int ldb) {
+  for (int j = 0; j < dimN / Step; j++, dst += Step, A += Step, B += Step) {
+    vecType mm[Step / VECTOR_LEN];
+    for (int n = 0; n < Step / VECTOR_LEN; n++) {
+      mm[n] = VECTOR_SET(agg.init());
+    }
+
+    for (int i = 0; i < dimM; i++) {
+      vecType *a = (vecType*)(A + i * lda);
+      vecType *b = (vecType*)(B + i * ldb);
+      for (int n = 0; n < Step / VECTOR_LEN; n++) {
+        mm[n] = agg.vecOp(mm[n], op.vecOp(a[n], b[n]));
+      }
+    }
+
+    vecType *result = (vecType*)(dst);
+    for (int n = 0; n < Step / VECTOR_LEN; n++) {
+      result[n] = sv.vecOp(result[n], mm[n]);
+    }
+  }
+
+  int remRow = dimN % Step;
+  if (remRow) {
+    hl_sse_column_op_with_rem<Step>(
+        agg, op, sv, dimM, remRow, dst, A, lda, B, ldb);
+  }
+}
+
+template <class Agg, class Op, class Saver>
+void hl_sse_matrix_column_op(Agg agg, Op op, Saver sv,
+                             int dimM, int dimN,
+                             real *dst,
+                             real *A, int lda,
+                             real *B, int ldb) {
+  if (dimN <= 16) {
+    hl_sse_matrix_column_op<16>(agg, op, sv, dimM, dimN, dst, A, lda, B, ldb);
+  } else if (dimN <= 32) {
+    hl_sse_matrix_column_op<32>(agg, op, sv, dimM, dimN, dst, A, lda, B, ldb);
+  } else if (dimN <= 1024 || dimM <= 512) {
+    hl_sse_matrix_column_op<64>(agg, op, sv, dimM, dimN, dst, A, lda, B, ldb);
+  } else {
+    hl_sse_matrix_column_op<1024>(agg, op, sv, dimM, dimN, dst, A, lda, B, ldb);
+  }
+}
+
+#endif /* HL_NEON_MATRIX_KERNEL_CUH_ */
diff --git a/paddle/math/SIMDFunctions.cpp b/paddle/math/SIMDFunctions.cpp
index 95219debf50e57407b668d315b91141d259fc779..60d37cef4000e9bd9ea78e5d72f065dadb4b4b61 100644
--- a/paddle/math/SIMDFunctions.cpp
+++ b/paddle/math/SIMDFunctions.cpp
@@ -13,10 +13,12 @@ See the License for the specific language governing permissions and
 limitations under the License. */
 
 #include "SIMDFunctions.h"
+#ifdef __SSE__
 #include <immintrin.h>
+#endif
 #include <algorithm>
 
-#ifndef __AVX__
+#ifdef __SSE__
 static void addto_sse(float* a, const float* b, size_t len) {
   int offset = len % 16;
   __m128 ma0, ma1, ma2, ma3;
@@ -125,7 +127,8 @@ static void col_max_sse(float* result,
   }
 }
 
-#else
+#elif defined(__AVX__)
+
 static void addto_avx(float* a, const float* b, size_t len) {
   int offset = len % 32;
 
@@ -357,15 +360,16 @@ static void decayL1_avx(
 
 #endif
 
-#ifndef __AVX__
+#ifdef __SSE__
 #define SIMD_INVOKE(func, ...) func##_sse(__VA_ARGS__)
-#else
+#elif __AVX__
 #define SIMD_INVOKE(func, ...) func##_avx(__VA_ARGS__)
 #endif
 
 namespace paddle {
 namespace simd {
 namespace internal {
+#ifdef __SSE__
 void addToImpl(float* a, const float* b, size_t len) {
   SIMD_INVOKE(addto, a, b, len);
 }
@@ -376,6 +380,7 @@ void batchAddToImpl(float* a, const float* b[], int batch, size_t len) {
 void colMaxImpl(float* result, const float* data, int dim, int numSamples) {
   SIMD_INVOKE(col_max, result, data, dim, numSamples);
 }
+#endif
 
 #ifdef __AVX__
 void decayL1AvxImpl(float* dst, float* src, float lambda, size_t len) {
diff --git a/paddle/math/SIMDFunctions.h b/paddle/math/SIMDFunctions.h
index 9b0a8719b287a2b88e966484090974586d64521f..9df02faa9027dae7f00b08d0f8cea2d459b0324e 100644
--- a/paddle/math/SIMDFunctions.h
+++ b/paddle/math/SIMDFunctions.h
@@ -128,17 +128,29 @@ void decayL1AvxImpl(
 
 template <>
 inline void addTo(float* a, const float* b, size_t len) {
+#ifdef __SSE__
   internal::addToImpl(a, b, len);
+#else
+  naive::addTo(a, b, len);
+#endif
 }
 
 template <>
 inline void batchAddTo(float* a, const float* b[], int batch, size_t len) {
+#ifdef __SSE__
   internal::batchAddToImpl(a, b, batch, len);
+#else
+  naive::batchAddTo(a, b, batch, len);
+#endif
 }
 
 template <>
 inline void colMax(float* result, const float* data, int dim, int numSamples) {
+#ifdef __SSE__
   internal::colMaxImpl(result, data, dim, numSamples);
+#else
+  naive::colMax(result, data, dim, numSamples);
+#endif
 }
 
 template <>
diff --git a/paddle/utils/tests/test_SIMDFlags.cpp b/paddle/utils/tests/test_SIMDFlags.cpp
index 8200a24ce7b7df75b48a89fbb7af15f304c5957f..185789c927be19385d6ddc7a1889b6cc56109d38 100644
--- a/paddle/utils/tests/test_SIMDFlags.cpp
+++ b/paddle/utils/tests/test_SIMDFlags.cpp
@@ -18,7 +18,8 @@ limitations under the License. */
 using namespace paddle;  // NOLINT
 
 TEST(SIMDFlags, gccTest) {
-#if (defined(__GNUC__) || defined(__GNUG__)) && !(defined(__clang__))
+#if (defined(__GNUC__) || defined(__GNUG__)) && !(defined(__clang__)) && \
+    !defined(__arm__)
   // clang-format off
   CHECK(!__builtin_cpu_supports("sse")    != HAS_SSE);
   CHECK(!__builtin_cpu_supports("sse2")   != HAS_SSE2);
@@ -43,4 +44,5 @@ TEST(SIMDFlags, normalPrint) {
   LOG(INFO) << "Has AVX:     " << std::boolalpha << HAS_AVX;
   LOG(INFO) << "Has AVX2:    " << std::boolalpha << HAS_AVX2;
   LOG(INFO) << "Has AVX512:  " << std::boolalpha << HAS_AVX512;
+  LOG(INFO) << "Has NEON:    " << std::boolalpha << HAS_NEON;
 }