diff --git a/.gitignore b/.gitignore
index e257658..f0bbcbb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,34 +1,3 @@
-# ---> C++
-# Prerequisites
-*.d
-
-# Compiled Object files
-*.slo
-*.lo
-*.o
-*.obj
-
-# Precompiled Headers
-*.gch
-*.pch
-
-# Compiled Dynamic libraries
-*.so
-*.dylib
-*.dll
-
-# Fortran module files
-*.mod
-*.smod
-
-# Compiled Static libraries
-*.lai
-*.la
-*.a
-*.lib
-
-# Executables
-*.exe
-*.out
-*.app
-
+.DS_Store
+build/
+.vscode/
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000..2999c7f
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,40 @@
+cmake_minimum_required(VERSION 3.15.2)
+# 设置项目名称与语言
+project(GCTL_AI VERSION 1.0)
+# 添加配置配件编写的函数
+include(CMakePackageConfigHelpers)
+
+# 添加编译选项
+option(GCTL_AI_OPENMP "Use the OpenMP library" OFF)
+
+message(STATUS "Platform: " ${CMAKE_HOST_SYSTEM_NAME})
+message(STATUS "Install prefix: " ${CMAKE_INSTALL_PREFIX})
+message(STATUS "Processor: " ${CMAKE_HOST_SYSTEM_PROCESSOR})
+message(STATUS "[GCTL_AI] Use the OpenMP library: " ${GCTL_AI_OPENMP})
+
+find_package(GCTL REQUIRED)
+include_directories(${GCTL_INC_DIR})
+
+find_package(GCTL_OPTIMIZATION REQUIRED)
+include_directories(${GCTL_OPTIMIZATION_INC_DIR})
+
+if(GCTL_AI_OPENMP)
+	if(NOT OpenMP_CXX_FOUND)
+		find_package(OpenMP REQUIRED)
+		set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+		set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+		include_directories(${OpenMP_CXX_INCLUDE_DIRS})
+	endif()
+endif()
+
+# 加入一个头文件配置，让cmake对源码进行操作
+configure_file(
+	"${PROJECT_SOURCE_DIR}/config.h.in"
+	"${PROJECT_SOURCE_DIR}/lib/dnn/config.h"
+	)
+
+# 添加库源文件地址
+add_subdirectory(lib)
+
+# 去掉注释编译示例
+add_subdirectory(examples)
\ No newline at end of file
diff --git a/GCTL_AIConfig.cmake.in b/GCTL_AIConfig.cmake.in
new file mode 100644
index 0000000..db9a325
--- /dev/null
+++ b/GCTL_AIConfig.cmake.in
@@ -0,0 +1,38 @@
+@PACKAGE_INIT@
+
+set(@PROJECT_NAME@_VERSION "@PROJECT_VERSION@")
+set_and_check(@PROJECT_NAME@_INSTALL_PREFIX "${PACKAGE_PREFIX_DIR}")
+set_and_check(@PROJECT_NAME@_INC_DIR "${PACKAGE_PREFIX_DIR}/include")
+set_and_check(@PROJECT_NAME@_INCLUDE_DIR "${PACKAGE_PREFIX_DIR}/include")
+set_and_check(@PROJECT_NAME@_LIB_DIR "${PACKAGE_PREFIX_DIR}/lib")
+set_and_check(@PROJECT_NAME@_LIBRARY_DIR "${PACKAGE_PREFIX_DIR}/lib")
+
+set(@PROJECT_NAME@_LIB gctl_ai)
+set(@PROJECT_NAME@_LIBRARY gctl_ai)
+
+set(@PROJECT_NAME@_OPENMP @GCTL_AI_OPENMP@)
+message(STATUS "[GCTL_AI] Use the OpenMP library:" @GCTL_AI_OPENMP@)
+
+if(NOT GCTL_FOUND)
+    find_package(GCTL REQUIRED)
+    include_directories(${GCTL_INC_DIR})
+endif()
+
+if(NOT GCTL_OPTIMIZATION_FOUND)
+    find_package(GCTL_OPTIMIZATION REQUIRED)
+    include_directories(${GCTL_OPTIMIZATION_INC_DIR})
+endif()
+
+if(@PROJECT_NAME@_OPENMP)
+    if(NOT OpenMP_CXX_FOUND)
+        find_package(OpenMP REQUIRED)
+        set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+        set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+        include_directories(${OpenMP_CXX_INCLUDE_DIRS})
+    endif()
+endif()
+
+# include target information
+include("${CMAKE_CURRENT_LIST_DIR}/@PROJECT_NAME@Targets.cmake")
+
+check_required_components(@PROJECT_NAME@)
\ No newline at end of file
diff --git a/README.md b/README.md
index 3399982..a442297 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,4 @@
-# gctl_ai
-
+# gctl_ai
+
+#### 介绍
+此库包含人工智能算法代码。
\ No newline at end of file
diff --git a/config.h.in b/config.h.in
new file mode 100644
index 0000000..58a7bce
--- /dev/null
+++ b/config.h.in
@@ -0,0 +1 @@
+#cmakedefine GCTL_AI_OPENMP
\ No newline at end of file
diff --git a/data/MNIST/mnist_database.h b/data/MNIST/mnist_database.h
new file mode 100644
index 0000000..ed1c1e6
--- /dev/null
+++ b/data/MNIST/mnist_database.h
@@ -0,0 +1,153 @@
+#ifndef _MNIST_DATABASE_H
+#define _MNIST_DATABASE_H
+
+#include "string"
+#include "iostream"
+#include "fstream"
+#include "vector"
+
+int ReverseInt(int i)  
+{  
+    unsigned char ch1, ch2, ch3, ch4;  
+    ch1 = i & 255;  
+    ch2 = (i >> 8) & 255;  
+    ch3 = (i >> 16) & 255;  
+    ch4 = (i >> 24) & 255;  
+    return((int)ch1 << 24) + ((int)ch2 << 16) + ((int)ch3 << 8) + ch4;  
+}
+
+class mnist_database
+{
+public:
+    mnist_database(std::string dir);
+    virtual ~mnist_database(){}
+
+    const std::vector<std::vector<double> > &train_images();
+    const std::vector<std::vector<double> > &test_images();
+    const std::vector<double> &train_labels();
+    const std::vector<double> &test_labels();
+    void image_dimension(int &rows, int &cols);
+
+private:
+    void read_mnist_images(std::ifstream &fs, std::vector<std::vector<double> > &images);
+    void read_mnist_labels(std::ifstream &fs, std::vector<double> &labels);
+
+private:
+    std::vector<std::vector<double> > train_images_, test_images_;
+    std::vector<double> train_labels_, test_labels_;
+};
+
+mnist_database::mnist_database(std::string dir)
+{
+    std::string file = dir + "/t10k-images.idx3-ubyte";
+    std::ifstream infile(file, std::ios::binary);
+    if (!infile) throw std::runtime_error("[mnist_database] Database is not found.");
+    read_mnist_images(infile, test_images_);
+    infile.close();
+
+    file = dir + "/t10k-labels.idx1-ubyte";
+    infile.open(file, std::ios::binary);
+    if (!infile) throw std::runtime_error("[mnist_database] Database is not found.");
+    read_mnist_labels(infile, test_labels_);
+    infile.close();
+
+    file = dir + "/train-images.idx3-ubyte";
+    infile.open(file, std::ios::binary);
+    if (!infile) throw std::runtime_error("[mnist_database] Database is not found.");
+    read_mnist_images(infile, train_images_);
+    infile.close();
+
+    file = dir + "/train-labels.idx1-ubyte";
+    infile.open(file, std::ios::binary);
+    if (!infile) throw std::runtime_error("[mnist_database] Database is not found.");
+    read_mnist_labels(infile, train_labels_);
+    infile.close();
+}
+
+const std::vector<std::vector<double> > &mnist_database::train_images()
+{
+    return train_images_;
+}
+
+const std::vector<std::vector<double> > &mnist_database::test_images()
+{
+    return test_images_;
+}
+
+const std::vector<double> &mnist_database::train_labels()
+{
+    return train_labels_;
+}
+
+const std::vector<double> &mnist_database::test_labels()
+{
+    return test_labels_;
+}
+
+void mnist_database::image_dimension(int &rows, int &cols)
+{
+    rows = cols = 28;
+    return;
+}
+
+void mnist_database::read_mnist_images(std::ifstream &fs, std::vector<std::vector<double> > &images)
+{
+    int magic_number = 0;  
+    int number_of_images = 0;  
+    int n_rows = 0;  
+    int n_cols = 0;  
+    unsigned char label;  
+    fs.read((char*)&magic_number, sizeof(magic_number));  
+    fs.read((char*)&number_of_images, sizeof(number_of_images));  
+    fs.read((char*)&n_rows, sizeof(n_rows));  
+    fs.read((char*)&n_cols, sizeof(n_cols));  
+    magic_number = ReverseInt(magic_number);  
+    number_of_images = ReverseInt(number_of_images);  
+    n_rows = ReverseInt(n_rows);  
+    n_cols = ReverseInt(n_cols);  
+
+    //std::cout << "magic number = " << magic_number << std::endl;  
+    //std::cout << "number of images = " << number_of_images << std::endl;  
+    //std::cout << "rows = " << n_rows << std::endl;  
+    //std::cout << "cols = " << n_cols << std::endl;  
+
+    std::vector<double> tp;
+    for (int i = 0; i < number_of_images; i++)  
+    {
+        tp.clear();
+        for (int r = 0; r < n_rows; r++)  
+        {  
+            for (int c = 0; c < n_cols; c++)  
+            {  
+                unsigned char image = 0;  
+                fs.read((char*)&image, sizeof(image));  
+                tp.push_back(image);  
+            }  
+        }  
+        images.push_back(tp);  
+    }
+    return;
+}
+
+void mnist_database::read_mnist_labels(std::ifstream &fs, std::vector<double> &labels)
+{
+    int magic_number = 0;  
+    int number_of_images = 0;  
+    fs.read((char*)&magic_number, sizeof(magic_number));  
+    fs.read((char*)&number_of_images, sizeof(number_of_images));  
+    magic_number = ReverseInt(magic_number);  
+    number_of_images = ReverseInt(number_of_images);
+
+    //std::cout << "magic number = " << magic_number << std::endl;  
+    //std::cout << "number of images = " << number_of_images << std::endl;  
+        
+    for (int i = 0; i < number_of_images; i++)  
+    {  
+        unsigned char label = 0;  
+        fs.read((char*)&label, sizeof(label));  
+        labels.push_back((double)label);  
+    }  
+    return;
+}
+
+#endif // _MNIST_DATABASE_H
\ No newline at end of file
diff --git a/data/MNIST/t10k-images.idx3-ubyte b/data/MNIST/t10k-images.idx3-ubyte
new file mode 100644
index 0000000..1170b2c
Binary files /dev/null and b/data/MNIST/t10k-images.idx3-ubyte differ
diff --git a/data/MNIST/t10k-labels.idx1-ubyte b/data/MNIST/t10k-labels.idx1-ubyte
new file mode 100644
index 0000000..d1c3a97
Binary files /dev/null and b/data/MNIST/t10k-labels.idx1-ubyte differ
diff --git a/data/MNIST/train-images.idx3-ubyte b/data/MNIST/train-images.idx3-ubyte
new file mode 100644
index 0000000..bbce276
Binary files /dev/null and b/data/MNIST/train-images.idx3-ubyte differ
diff --git a/data/MNIST/train-labels.idx1-ubyte b/data/MNIST/train-labels.idx1-ubyte
new file mode 100644
index 0000000..d6b4c5d
Binary files /dev/null and b/data/MNIST/train-labels.idx1-ubyte differ
diff --git a/data/saved_networks/mnist_m1.gctl.dnn b/data/saved_networks/mnist_m1.gctl.dnn
new file mode 100644
index 0000000..5b17d16
Binary files /dev/null and b/data/saved_networks/mnist_m1.gctl.dnn differ
diff --git a/data/saved_networks/mnist_m2.gctl.dnn b/data/saved_networks/mnist_m2.gctl.dnn
new file mode 100644
index 0000000..a44ef58
Binary files /dev/null and b/data/saved_networks/mnist_m2.gctl.dnn differ
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
new file mode 100644
index 0000000..3ff7f9d
--- /dev/null
+++ b/examples/CMakeLists.txt
@@ -0,0 +1,13 @@
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --std=c++11 -O3")
+set(EXECUTABLE_OUTPUT_PATH ${PROJECT_BINARY_DIR}/bin/examples)
+
+macro(add_example name)
+    add_executable(${name} ${name}.cpp)
+    target_link_libraries(${name} PUBLIC gctl_ai)
+endmacro()
+
+add_example(ex1)
+add_example(ex2)
+add_example(ex_mnist)
+add_example(ex_mnist2)
+add_example(ex_mnist3)
\ No newline at end of file
diff --git a/examples/ex1.cpp b/examples/ex1.cpp
new file mode 100644
index 0000000..1bed7e0
--- /dev/null
+++ b/examples/ex1.cpp
@@ -0,0 +1,120 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "../lib/dnn.h"
+using namespace gctl;
+
+void data_generator(const matrix<double> &train_obs, matrix<double> &train_tar)
+{
+    for (int j = 0; j < train_obs.col_size(); j++)
+    {
+        train_tar[0][j] = sqrt(train_obs[0][j]*train_obs[0][j] + train_obs[1][j]*train_obs[1][j]);
+    }
+    return;
+}
+
+int main(int argc, char const *argv[]) try
+{
+    // Prepare the data. In this example, we try to learn the sin() function.
+    matrix<double> train_obs(2, 1000), train_tar(1, 1000), pre_obs(2, 10), pre_tar(1, 10), predicts(1, 10);
+
+    unsigned int seed = 101;
+    srand(seed);
+    for (int j = 0; j < 1000; j++)
+    {
+        for (int i = 0; i < 2; i++)
+        {
+            train_obs[i][j] = random(0.0, 1.0);
+        }
+    }
+
+    for (int j = 0; j < 10; j++)
+    {
+        for (int i = 0; i < 2; i++)
+        {
+            pre_obs[i][j] = random(0.0, 1.0);
+        }
+    }
+
+    data_generator(train_obs, train_tar);
+    data_generator(pre_obs, pre_tar);
+
+    dnn my_nn("Ex-1");
+    my_nn.add_hind_layer(2, 100, FullyConnected, Identity);
+    my_nn.add_hind_layer(100, 100, FullyConnected, PReLU);
+    my_nn.add_hind_layer(100, 100, FullyConnected, PReLU);
+    my_nn.add_hind_layer(100, 1, FullyConnected, Identity);
+    my_nn.add_output_layer(RegressionMSE);
+    my_nn.add_train_set(train_obs, train_tar, 200);
+
+    my_nn.show_network();
+    my_nn.init_network(0.0, 0.1, seed);
+
+    sgd_para my_para = my_nn.default_sgd_para();
+    my_nn.train_network(my_para, gctl::ADAM);
+    //lgd_para my_para = my_nn.default_lgd_para();
+    //my_para.flight_times = 5000;
+    //my_para.lambda = 5e-5;
+    //my_para.epsilon = 1e-5;
+    //my_para.batch = 10;
+    //my_nn.train_network(my_para, gctl::LGD);
+
+    my_nn.predict(pre_obs, predicts);
+
+    double diff = 0;
+    for (int i = 0; i < 1; i++)
+    {
+        for (int j = 0; j < 10; j++)
+        {
+            diff = std::max(fabs(predicts[i][j] - pre_tar[i][j]), diff);
+        }
+    }
+    std::clog << "Max difference = " << diff << "\n";
+/*
+    my_nn.save_network("ex1");
+
+    dnn file_nn("File NN");
+    file_nn.load_network("ex1");
+    file_nn.show_network();
+    file_nn.predict(pre_obs, predicts);
+
+    diff = 0;
+    for (int i = 0; i < 1; i++)
+    {
+        for (int j = 0; j < 10; j++)
+        {
+            diff = std::max(fabs(predicts[i][j] - pre_tar[i][j]), diff);
+        }
+    }
+    std::clog << "Max difference = " << diff << "\n";
+*/
+    return 0;
+}
+catch (std::exception &e)
+{
+    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}
\ No newline at end of file
diff --git a/examples/ex2.cpp b/examples/ex2.cpp
new file mode 100644
index 0000000..91e27bd
--- /dev/null
+++ b/examples/ex2.cpp
@@ -0,0 +1,66 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "../lib/dnn.h"
+using namespace gctl;
+
+int main(int argc, char const *argv[]) try
+{
+    convolution cnl_layer(0, 3, 10, 10, 3, 3, 2, 2, Valid, ReLU);
+    std::clog << cnl_layer.layer_info() << std::endl;
+
+    matrix<double> t(100, 3);
+    std::clog << "T = \n";
+    for (size_t c = 0; c < 3; c++)
+    {
+        for (size_t i = 0; i < 10; i++)
+        {
+            for (size_t j = 0; j < 10; j++)
+            {
+                t[i*10+j][c] = i*10 + j + 1;
+                std::clog << t[i*10+j][c] << " ";
+            }
+            std::clog << "\n";
+        }
+    }
+
+    array<double> weights(28);
+    for (size_t i = 0; i < 28; i++)
+    {
+        weights[i] = 1.0;
+    }
+    
+    cnl_layer.forward_propagation(weights, t);
+
+    const matrix<double> &d = cnl_layer.forward_propagation_data();
+    d.show(std::clog);
+    return 0;
+}
+catch (std::exception &e)
+{
+    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}
\ No newline at end of file
diff --git a/examples/ex_mnist.cpp b/examples/ex_mnist.cpp
new file mode 100644
index 0000000..b587502
--- /dev/null
+++ b/examples/ex_mnist.cpp
@@ -0,0 +1,116 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "../data/MNIST/mnist_database.h"
+#include "../lib/dnn.h"
+
+using namespace gctl;
+
+int main(int argc, char const *argv[]) try
+{
+    mnist_database data("data/MNIST");
+
+    matrix<double> train_obs(784, 60000), train_lab(10, 60000, 0.0);
+    matrix<double> test_obs(784, 10000), test_lab(10, 10000, 0.0), predicts(10, 10000);
+
+    const std::vector<std::vector<double> > &dt_obs = data.train_images();
+    for (size_t i = 0; i < 60000; i++)
+    {
+        for (size_t j = 0; j < 784; j++)
+        {
+            train_obs[j][i] = dt_obs[i][j]/255.0;
+        }
+    }
+
+    const std::vector<std::vector<double> > &dt_obs2 = data.test_images();
+    for (size_t i = 0; i < 10000; i++)
+    {
+        for (size_t j = 0; j < 784; j++)
+        {
+            test_obs[j][i] = dt_obs2[i][j]/255.0;
+        }
+    }
+
+    const std::vector<double> &dt_lab = data.train_labels();
+    for (size_t i = 0; i < 60000; i++)
+    {
+        train_lab[dt_lab[i]][i] = 1.0;
+    }
+
+    const std::vector<double> &dt_lab2 = data.test_labels();
+    for (size_t i = 0; i < 10000; i++)
+    {
+        test_lab[dt_lab2[i]][i] = 1.0;
+    }
+
+    dnn my_nn("Ex-MNIST");
+    my_nn.add_hind_layer(784, 800, FullyConnected, PReLU);
+    my_nn.add_hind_layer(800, 10, FullyConnected, SoftMax);
+    my_nn.add_output_layer(MultiClassEntropy);
+    my_nn.add_train_set(train_obs, train_lab, 1000);
+    my_nn.init_network(0.0, 0.1);
+
+    //sgd_para my_para = my_nn.default_sgd_para();
+    //my_para.alpha = 0.01;
+    //my_para.epsilon = 1e-5;
+    //my_nn.train_network(my_para, gctl::ADAM);
+    lgd_para my_para = my_nn.default_lgd_para();
+    my_para.flight_times = 1000;
+    my_para.alpha = 0.08;
+    my_para.beta = 1.8;
+    my_nn.train_network(my_para);
+
+    my_nn.predict(test_obs, predicts);
+    my_nn.save_network("data/saved_networks/mnist_m1");
+
+    int wrong_predicts = 0;
+    int test_id, pre_id;
+    double test_scr, pre_scr;
+    for (size_t j = 0; j < 10000; j++)
+    {
+        test_id = pre_id = 0;
+        test_scr = pre_scr = 0;
+        for (size_t i = 0; i < 10; i++)
+        {
+            if (test_lab[i][j] > test_scr) {test_scr = test_lab[i][j]; test_id = i;}
+            if (predicts[i][j] > pre_scr) {pre_scr = predicts[i][j]; pre_id = i;}
+        }
+
+        if (test_id != pre_id)
+        {
+            wrong_predicts++;   
+        }
+    }
+
+    std::cout << "Correct Rate = " << (10000 - wrong_predicts)/100.0 << "%\n";
+
+    return 0;
+}
+catch (std::exception &e)
+{
+    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}
\ No newline at end of file
diff --git a/examples/ex_mnist2.cpp b/examples/ex_mnist2.cpp
new file mode 100644
index 0000000..f5b24f9
--- /dev/null
+++ b/examples/ex_mnist2.cpp
@@ -0,0 +1,116 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "../data/MNIST/mnist_database.h"
+#include "../lib/dnn.h"
+
+using namespace gctl;
+
+int main(int argc, char const *argv[]) try
+{
+    mnist_database data("data/MNIST");
+
+    matrix<double> train_obs(784, 60000), train_lab(10, 60000, 0.0);
+    matrix<double> test_obs(784, 10000), test_lab(10, 10000, 0.0), predicts(10, 10000);
+
+    const std::vector<std::vector<double> > &dt_obs = data.train_images();
+    for (size_t i = 0; i < 60000; i++)
+    {
+        for (size_t j = 0; j < 784; j++)
+        {
+            train_obs[j][i] = dt_obs[i][j]/255.0;
+        }
+    }
+
+    const std::vector<std::vector<double> > &dt_obs2 = data.test_images();
+    for (size_t i = 0; i < 10000; i++)
+    {
+        for (size_t j = 0; j < 784; j++)
+        {
+            test_obs[j][i] = dt_obs2[i][j]/255.0;
+        }
+    }
+
+    const std::vector<double> &dt_lab = data.train_labels();
+    for (size_t i = 0; i < 60000; i++)
+    {
+        train_lab[dt_lab[i]][i] = 1.0;
+    }
+
+    const std::vector<double> &dt_lab2 = data.test_labels();
+    for (size_t i = 0; i < 10000; i++)
+    {
+        test_lab[dt_lab2[i]][i] = 1.0;
+    }
+
+    dnn my_nn("Ex-MNIST");
+    my_nn.add_hind_layer(1, 28, 28, 4, 4, 2, 2, Convolution, Same, PReLU);
+    my_nn.add_hind_layer(169, 256, FullyConnected, PReLU);
+    my_nn.add_hind_layer(256, 10, FullyConnected, SoftMax);
+    my_nn.add_output_layer(MultiClassEntropy);
+    my_nn.add_train_set(train_obs, train_lab, 5000);
+    my_nn.init_network(0.0, 0.1);
+
+    sgd_para my_para = my_nn.default_sgd_para();
+    my_para.alpha = 0.01;
+    my_para.epsilon = 1e-5;
+    my_nn.train_network(my_para, gctl::ADAM);
+    //lgd_para my_para = my_nn.default_lgd_para();
+    //my_para.flight_times = 2000;
+    //my_para.alpha = 0.1;
+    //my_nn.train_network(my_para);
+
+    my_nn.predict(test_obs, predicts);
+    my_nn.save_network("data/saved_networks/mnist_m2");
+
+    int wrong_predicts = 0;
+    int test_id, pre_id;
+    double test_scr, pre_scr;
+    for (size_t j = 0; j < 10000; j++)
+    {
+        test_id = pre_id = 0;
+        test_scr = pre_scr = 0;
+        for (size_t i = 0; i < 10; i++)
+        {
+            if (test_lab[i][j] > test_scr) {test_scr = test_lab[i][j]; test_id = i;}
+            if (predicts[i][j] > pre_scr) {pre_scr = predicts[i][j]; pre_id = i;}
+        }
+
+        if (test_id != pre_id)
+        {
+            wrong_predicts++;   
+        }
+    }
+
+    std::cout << "Correct Rate = " << (10000 - wrong_predicts)/100.0 << "%\n";
+
+    return 0;
+}
+catch (std::exception &e)
+{
+    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}
\ No newline at end of file
diff --git a/examples/ex_mnist3.cpp b/examples/ex_mnist3.cpp
new file mode 100644
index 0000000..d7fae16
--- /dev/null
+++ b/examples/ex_mnist3.cpp
@@ -0,0 +1,87 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "../data/MNIST/mnist_database.h"
+#include "../lib/dnn.h"
+
+using namespace gctl;
+
+int main(int argc, char const *argv[]) try
+{
+    mnist_database data("data/MNIST");
+    dnn my_nn("Ex-MNIST");
+    my_nn.load_network("data/saved_networks/mnist_m1.gctl.dnn");
+    my_nn.show_network();
+    my_nn.save_layer2text(0, "data/saved_networks/mnist_m1_layer1");
+    my_nn.save_layer2text(1, "data/saved_networks/mnist_m1_layer2");
+/*    
+    matrix<double> test_obs(784, 10000), test_lab(10, 10000, 0.0), predicts(10, 10000);
+
+    const std::vector<std::vector<double> > &dt_obs2 = data.test_images();
+    for (size_t i = 0; i < 10000; i++)
+    {
+        for (size_t j = 0; j < 784; j++)
+        {
+            test_obs[j][i] = dt_obs2[i][j]/255.0;
+        }
+    }
+
+    const std::vector<double> &dt_lab2 = data.test_labels();
+    for (size_t i = 0; i < 10000; i++)
+    {
+        test_lab[dt_lab2[i]][i] = 1.0;
+    }
+
+    my_nn.predict(test_obs, predicts);
+
+    int wrong_predicts = 0;
+    int test_id, pre_id;
+    double test_scr, pre_scr;
+    for (size_t j = 0; j < 10000; j++)
+    {
+        test_id = pre_id = 0;
+        test_scr = pre_scr = 0;
+        for (size_t i = 0; i < 10; i++)
+        {
+            if (test_lab[i][j] > test_scr) {test_scr = test_lab[i][j]; test_id = i;}
+            if (predicts[i][j] > pre_scr) {pre_scr = predicts[i][j]; pre_id = i;}
+        }
+
+        if (test_id != pre_id)
+        {
+            wrong_predicts++;   
+        }
+    }
+
+    std::cout << "Correct Rate = " << (10000 - wrong_predicts)/100.0 << "%\n";
+*/
+    return 0;
+}
+catch (std::exception &e)
+{
+    GCTL_ShowWhatError(e.what(), GCTL_ERROR_ERROR, 0, 0, 0);
+}
\ No newline at end of file
diff --git a/installer b/installer
new file mode 100755
index 0000000..c49ddcc
--- /dev/null
+++ b/installer
@@ -0,0 +1,54 @@
+#!/bin/bash
+
+if [[ $# == 0 || ${1} == "help" ]]; then
+	echo "Compiles executables/libraries and maintains installed files. Two tools 'Cmake' and 'stow' are empolyed here. For more information, see https://cmake.org and https://www.gnu.org/software/stow/."
+	echo ""
+	echo "School of Earth Sciences, Zhejiang University"
+	echo "Yi Zhang (yizhang-geo@zju.edu.cn)"
+	echo ""
+	echo "Usage: ./config.sh [option] [Cmake options]"
+	echo ""
+	echo "Options:"
+	echo "(1) configure: Configure Cmake project(s). This option could take extra Cmake options as in <option>=<value>."
+	echo "(2)     build: Build executables/libraries."
+	echo "(3)   install: Install executables/libraries to the directory of CMAKE_INSTALL_PREFIX and sym-links them to the target address. This offers a quick and clean remove of the installed files."
+	echo "(4)     clean: Clean build/ folder(s)."
+	echo "(5) uninstall: Delete the installed files and sym-links."
+	echo "(6)      info: Print out current setups."
+	echo "(7)      help: Show help information."
+	exit 0
+fi
+
+package=gctl_potential
+address=/opt/stow
+taress=/usr/local
+option="-DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=${address}/${package}"
+
+if [[ $# -gt 1 ]]; then
+	for opt in "$@"; do
+		if [[ ${opt} != "configure" ]]; then
+			option="${option} -D${opt}"
+		fi
+	done
+fi
+
+if [[ ${1} == "configure" && ! -d "build/" ]]; then
+	mkdir build && cd build && cmake .. ${option}
+elif [[ ${1} == "configure" ]]; then
+	cd build && rm -rf * && cmake .. ${option}
+elif [[ ${1} == "build" ]]; then
+	cd build && make
+elif [[ ${1} == "install" ]]; then
+	cd build && sudo make install
+	sudo stow --dir=${address} --target=${taress} -S ${package}
+elif [[ ${1} == "clean" ]]; then
+	rm -rf build/
+elif [[ ${1} == "uninstall" ]]; then
+	sudo stow --dir=${address} --target=${taress} -D ${package}
+	sudo rm -rf ${address}/${package}
+elif [[ ${1} == "info" ]]; then
+	echo "package name:" ${package}
+	echo "stow address:" ${address}
+	echo "target address:" ${taress}
+	echo "Cmake options:" ${option}
+fi
diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
new file mode 100644
index 0000000..cc4c23d
--- /dev/null
+++ b/lib/CMakeLists.txt
@@ -0,0 +1,69 @@
+# 设置编译选项
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} --std=c++11 -O3")
+# 设置库文件的输出地址
+set(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib)
+# 设定库源文件文件夹
+aux_source_directory(dnn/ GCTL_DNN_SRC)
+
+# 以下部分为库的编译
+# 注意目标名必须唯一 所以不能直接生成相同名称的动态库与静态库
+# 注意此处不必为目标名称添加lib前缀和相应后缀，cmake会自行添加
+add_library(gctl_ai SHARED ${GCTL_DNN_SRC})
+# 首先添加静态库的生成命令
+add_library(gctl_ai_static STATIC ${GCTL_DNN_SRC})
+# 设置静态库的输出名称从而获得与动态库名称相同的静态库
+set_target_properties(gctl_ai_static PROPERTIES OUTPUT_NAME "gctl_ai")
+# 设置输出目标属性以同时输出动态库与静态库
+set_target_properties(gctl_ai PROPERTIES CLEAN_DIRECT_OUTPUT 1)
+set_target_properties(gctl_ai_static PROPERTIES CLEAN_DIRECT_OUTPUT 1)
+# 设置动态库的版本号
+set_target_properties(gctl_ai PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR}.${PROJECT_VERSION_MINOR})
+# 设置动态库的运行搜索地址
+set_target_properties(gctl_ai PROPERTIES INSTALL_RPATH /usr/local/lib)
+set_target_properties(gctl_ai_static PROPERTIES INSTALL_RPATH /usr/local/lib)
+
+#连接动态库
+target_link_libraries(gctl_ai ${GCTL_LIB})
+target_link_libraries(gctl_ai_static ${GCTL_LIB})
+target_link_libraries(gctl_ai ${GCTL_OPTIMIZATION_LIB})
+target_link_libraries(gctl_ai_static ${GCTL_OPTIMIZATION_LIB})
+
+if(GCTL_AI_OPENMP)
+	target_link_libraries(gctl_ai OpenMP::OpenMP_CXX)
+	target_link_libraries(gctl_ai_static OpenMP::OpenMP_CXX)
+endif()
+
+set(CONFIG_FILE_PATH lib/cmake/${PROJECT_NAME})
+
+configure_package_config_file(${PROJECT_SOURCE_DIR}/${PROJECT_NAME}Config.cmake.in 
+	${CMAKE_BINARY_DIR}/${PROJECT_NAME}Config.cmake
+	INSTALL_DESTINATION ${CONFIG_FILE_PATH})
+
+write_basic_package_version_file(${CMAKE_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake
+	VERSION ${PROJECT_VERSION}
+	COMPATIBILITY SameMajorVersion)
+
+# 库的安装命令
+if(WIN32)
+	install(TARGETS gctl_ai DESTINATION lib)
+	install(TARGETS gctl_ai_static DESTINATION lib)
+else()
+	install(TARGETS gctl_ai gctl_ai_static
+		EXPORT ${PROJECT_NAME}Targets
+		LIBRARY DESTINATION lib
+		ARCHIVE DESTINATION lib)
+	install(EXPORT ${PROJECT_NAME}Targets
+		DESTINATION ${CONFIG_FILE_PATH})
+	install(FILES
+		${CMAKE_BINARY_DIR}/${PROJECT_NAME}Config.cmake
+		${CMAKE_BINARY_DIR}/${PROJECT_NAME}ConfigVersion.cmake
+		DESTINATION ${CONFIG_FILE_PATH})
+endif()
+
+# 头文件安装命令
+file(GLOB GCTL_HEAD *.h)
+file(GLOB GCTL_AI_HEAD dnn/*.h)
+
+install(FILES ${GCTL_HEAD} DESTINATION include/gctl)
+install(FILES ${GCTL_AI_HEAD} DESTINATION include/gctl/dnn)
\ No newline at end of file
diff --git a/lib/dnn.h b/lib/dnn.h
new file mode 100644
index 0000000..a7ad151
--- /dev/null
+++ b/lib/dnn.h
@@ -0,0 +1,36 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_H
+#define _GCTL_DNN_H
+
+#include "dnn/activation.h"
+#include "dnn/hlayer.h"
+#include "dnn/olayer.h"
+#include "dnn/dnn.h"
+
+#endif // _GCTL_DNN_H
\ No newline at end of file
diff --git a/lib/dnn/activation.cpp b/lib/dnn/activation.cpp
new file mode 100644
index 0000000..ce7b98c
--- /dev/null
+++ b/lib/dnn/activation.cpp
@@ -0,0 +1,42 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "activation.h"
+
+gctl::dnn_activation::dnn_activation() {}
+
+gctl::dnn_activation::~dnn_activation() {}
+
+std::string gctl::dnn_activation::activation_name() const
+{
+    return "Unknown";
+}
+
+gctl::activation_type_e gctl::dnn_activation::activation_type() const
+{
+    return Unknown;
+}
\ No newline at end of file
diff --git a/lib/dnn/activation.h b/lib/dnn/activation.h
new file mode 100644
index 0000000..319027d
--- /dev/null
+++ b/lib/dnn/activation.h
@@ -0,0 +1,70 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_ACTIVATION_H
+#define _GCTL_DNN_ACTIVATION_H
+
+#include "gctl/core.h"
+#include "gctl/io.h"
+#include "gctl/algorithm.h"
+
+#include "config.h"
+#ifdef GCTL_AI_OPENMP
+#include "omp.h"
+#endif // GCTL_AI_OPENMP
+
+namespace gctl
+{
+    enum activation_type_e
+    {
+        Identity,
+        Mish,
+        ReLU,
+        PReLU,
+        Sigmoid,
+        SoftMax,
+        Tanh,
+        Unknown,
+    };
+
+    class dnn_activation
+    {
+    public:
+        dnn_activation();
+        virtual ~dnn_activation();
+
+        virtual void activate(const matrix<double> &z, matrix<double> &a) = 0;
+        virtual void apply_jacobian(const matrix<double> &z, const matrix<double> &a, 
+            const matrix<double> &f, matrix<double> &g) = 0;
+        virtual std::string activation_name() const;
+        virtual activation_type_e activation_type() const;
+    };
+}
+
+
+
+#endif // _GCTL_DNN_ACTIVATION_H
\ No newline at end of file
diff --git a/lib/dnn/activation_identity.cpp b/lib/dnn/activation_identity.cpp
new file mode 100644
index 0000000..5c92339
--- /dev/null
+++ b/lib/dnn/activation_identity.cpp
@@ -0,0 +1,64 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "activation_identity.h"
+
+gctl::identity::identity() {}
+
+gctl::identity::~identity() {}
+
+void gctl::identity::activate(const matrix<double> &z, matrix<double> &a)
+{
+    // a = activation(z) = z
+    // Z = [z1, ..., zn], A = [a1, ..., an], n observations
+    a.resize(z.row_size(), z.col_size());
+    a = z;
+    return;
+}
+
+void gctl::identity::apply_jacobian(const matrix<double> &z, 
+    const matrix<double> &a, const matrix<double> &f, matrix<double> &g)
+{
+    // Apply the Jacobian matrix J to a vector f
+    // J = d_a / d_z = I
+    // g = J * f = f
+    // Z = [z1, ..., zn], G = [g1, ..., gn], F = [f1, ..., fn]
+    // Note: When entering this function, Z and G may point to the same matrix
+    g.resize(a.row_size(), a.col_size());
+    g = f;
+    return;
+}
+
+std::string gctl::identity::activation_name() const
+{
+    return "Identity";
+}
+
+gctl::activation_type_e gctl::identity::activation_type() const
+{
+    return Identity;
+}
\ No newline at end of file
diff --git a/lib/dnn/activation_identity.h b/lib/dnn/activation_identity.h
new file mode 100644
index 0000000..3ff4181
--- /dev/null
+++ b/lib/dnn/activation_identity.h
@@ -0,0 +1,49 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_ACTIVATION_IDENTITY_H
+#define _GCTL_ACTIVATION_IDENTITY_H
+
+#include "activation.h"
+
+namespace gctl
+{
+    class identity : public dnn_activation
+    {
+    public:
+        identity();
+        virtual ~identity();
+
+        void activate(const matrix<double> &z, matrix<double> &a);
+        void apply_jacobian(const matrix<double> &z, const matrix<double> &a, 
+            const matrix<double> &f, matrix<double> &g);
+        std::string activation_name() const;
+        activation_type_e activation_type() const;
+    };
+}
+
+#endif // _GCTL_ACTIVATION_IDENTITY_H
\ No newline at end of file
diff --git a/lib/dnn/activation_mish.cpp b/lib/dnn/activation_mish.cpp
new file mode 100644
index 0000000..db77c4e
--- /dev/null
+++ b/lib/dnn/activation_mish.cpp
@@ -0,0 +1,102 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "activation_mish.h"
+
+gctl::mish::mish() {}
+
+gctl::mish::~mish() {}
+
+void gctl::mish::activate(const matrix<double> &z, matrix<double> &a)
+{
+    // Mish(x) = x * tanh(softplus(x))
+    // softplus(x) = log(1 + exp(x))
+    // a = activation(z) = Mish(z)
+    // Z = [z1, ..., zn], A = [a1, ..., an], n observations
+    // h(x) = tanh(softplus(x)) = (1 + exp(x))^2 - 1
+    //                            ------------------
+    //                            (1 + exp(x))^2 + 1
+    // Let s = exp(-abs(x)), t = 1 + s
+    // If x >= 0, then h(x) = (t^2 - s^2) / (t^2 + s^2)
+    // If x <= 0, then h(x) = (t^2 - 1) / (t^2 + 1)
+    a.resize(z.row_size(), z.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < a.row_size(); i++)
+    {
+        for (j = 0; j < a.col_size(); j++)
+        {
+            a[i][j] = z[i][j]*std::tanh(log(1.0 + exp(z[i][j])));
+        }
+    }
+    return;
+}
+
+void gctl::mish::apply_jacobian(const matrix<double> &z, 
+    const matrix<double> &a, const matrix<double> &f, matrix<double> &g)
+{
+    // Apply the Jacobian matrix J to a vector f
+    // J = d_a / d_z = diag(Mish'(z))
+    // g = J * f = Mish'(z) .* f
+    // Z = [z1, ..., zn], G = [g1, ..., gn], F = [f1, ..., fn]
+    // Note: When entering this function, Z and G may point to the same matrix
+    // Let h(x) = tanh(softplus(x))
+    // Mish'(x) = h(x) + x * h'(x)
+    // h'(x) = tanh'(softplus(x)) * softplus'(x)
+    //       = [1 - h(x)^2] * exp(x) / (1 + exp(x))
+    //       = [1 - h(x)^2] / (1 + exp(-x))
+    // Mish'(x) = h(x) + [x - Mish(x) * h(x)] / (1 + exp(-x))
+    // A = Mish(Z) = Z .* h(Z) => h(Z) = A ./ Z, h(0) = 0.6
+    g.resize(a.row_size(), a.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < g.row_size(); i++)
+    {
+        for (j = 0; j < g.col_size(); j++)
+        {
+            g[i][j] = std::tanh(log(1.0 + exp(z[i][j])));   
+        }
+
+        for (j = 0; j < g.col_size(); j++)
+        {
+            g[i][j] = f[i][j]*(g[i][j] + (z[i][j] - a[i][j]*g[i][j])/(1.0 + exp(-1.0*z[i][j])));   
+        }
+    }
+    return;
+}
+
+std::string gctl::mish::activation_name() const
+{
+    return "Mish";
+}
+
+gctl::activation_type_e gctl::mish::activation_type() const
+{
+    return Mish;
+}
\ No newline at end of file
diff --git a/lib/dnn/activation_mish.h b/lib/dnn/activation_mish.h
new file mode 100644
index 0000000..91c7dcd
--- /dev/null
+++ b/lib/dnn/activation_mish.h
@@ -0,0 +1,49 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_ACTIVATION_MISH_H
+#define _GCTL_ACTIVATION_MISH_H
+
+#include "activation.h"
+
+namespace gctl
+{
+    class mish : public dnn_activation
+    {
+    public:
+        mish();
+        virtual ~mish();
+
+        void activate(const matrix<double> &z, matrix<double> &a);
+        void apply_jacobian(const matrix<double> &z, const matrix<double> &a, 
+            const matrix<double> &f, matrix<double> &g);
+        std::string activation_name() const;
+        activation_type_e activation_type() const;
+    };
+}
+
+#endif // _GCTL_ACTIVATION_MISH_H
\ No newline at end of file
diff --git a/lib/dnn/activation_prelu.cpp b/lib/dnn/activation_prelu.cpp
new file mode 100644
index 0000000..fd89f29
--- /dev/null
+++ b/lib/dnn/activation_prelu.cpp
@@ -0,0 +1,84 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "activation_prelu.h"
+
+gctl::prelu::prelu() {}
+
+gctl::prelu::~prelu() {}
+
+void gctl::prelu::activate(const matrix<double> &z, matrix<double> &a)
+{
+    // a = activation(z) = max(z, 0)
+    // Z = [z1, ..., zn], A = [a1, ..., an], n observations
+    a.resize(z.row_size(), z.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < a.row_size(); i++)
+    {
+        for (j = 0; j < a.col_size(); j++)
+        {
+            if (z[i][j] > 0.0) a[i][j] = z[i][j];
+            else a[i][j] = 0.01*z[i][j];
+        }
+    }
+    return;
+}
+
+void gctl::prelu::apply_jacobian(const matrix<double> &z, 
+    const matrix<double> &a, const matrix<double> &f, matrix<double> &g)
+{
+    // Apply the Jacobian matrix J to a vector f
+    // J = d_a / d_z = diag(sign(a)) = diag(a > 0)
+    // g = J * f = (a > 0) .* f
+    // Z = [z1, ..., zn], G = [g1, ..., gn], F = [f1, ..., fn]
+    // Note: When entering this function, Z and G may point to the same matrix
+    g.resize(a.row_size(), a.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < g.row_size(); i++)
+    {
+        for (j = 0; j < g.col_size(); j++)
+        {
+            if (a[i][j] > 0.0) g[i][j] = f[i][j];
+            else g[i][j] = 0.01*f[i][j];
+        }
+    }
+    return;
+}
+
+std::string gctl::prelu::activation_name() const
+{
+    return "PReLU";
+}
+
+gctl::activation_type_e gctl::prelu::activation_type() const
+{
+    return PReLU;
+}
\ No newline at end of file
diff --git a/lib/dnn/activation_prelu.h b/lib/dnn/activation_prelu.h
new file mode 100644
index 0000000..d6143cc
--- /dev/null
+++ b/lib/dnn/activation_prelu.h
@@ -0,0 +1,49 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_ACTIVATION_PRELU_H
+#define _GCTL_ACTIVATION_PRELU_H
+
+#include "activation.h"
+
+namespace gctl
+{
+    class prelu : public dnn_activation
+    {
+    public:
+        prelu();
+        virtual ~prelu();
+
+        void activate(const matrix<double> &z, matrix<double> &a);
+        void apply_jacobian(const matrix<double> &z, const matrix<double> &a, 
+            const matrix<double> &f, matrix<double> &g);
+        std::string activation_name() const;
+        activation_type_e activation_type() const;
+    };
+}
+
+#endif // _GCTL_ACTIVATION_PRELU_H
diff --git a/lib/dnn/activation_relu.cpp b/lib/dnn/activation_relu.cpp
new file mode 100644
index 0000000..2a86909
--- /dev/null
+++ b/lib/dnn/activation_relu.cpp
@@ -0,0 +1,84 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "activation_relu.h"
+
+gctl::relu::relu() {}
+
+gctl::relu::~relu() {}
+
+void gctl::relu::activate(const matrix<double> &z, matrix<double> &a)
+{
+    // a = activation(z) = max(z, 0)
+    // Z = [z1, ..., zn], A = [a1, ..., an], n observations
+    a.resize(z.row_size(), z.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < a.row_size(); i++)
+    {
+        for (j = 0; j < a.col_size(); j++)
+        {
+            a[i][j] = std::max(z[i][j], 0.0);   
+        }
+    }
+    return;
+}
+
+void gctl::relu::apply_jacobian(const matrix<double> &z, 
+    const matrix<double> &a, const matrix<double> &f, matrix<double> &g)
+{
+    // Apply the Jacobian matrix J to a vector f
+    // J = d_a / d_z = diag(sign(a)) = diag(a > 0)
+    // g = J * f = (a > 0) .* f
+    // Z = [z1, ..., zn], G = [g1, ..., gn], F = [f1, ..., fn]
+    // Note: When entering this function, Z and G may point to the same matrix
+    g.resize(a.row_size(), a.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < g.row_size(); i++)
+    {
+        for (j = 0; j < g.col_size(); j++)
+        {
+            if (a[i][j] > 0.0) g[i][j] = f[i][j];
+            else g[i][j] = 0.0;
+            //g[i][j] = std::max(a[i][j], 0.0)*f[i][j];   
+        }
+    }
+    return;
+}
+
+std::string gctl::relu::activation_name() const
+{
+    return "ReLU";
+}
+
+gctl::activation_type_e gctl::relu::activation_type() const
+{
+    return ReLU;
+}
\ No newline at end of file
diff --git a/lib/dnn/activation_relu.h b/lib/dnn/activation_relu.h
new file mode 100644
index 0000000..ec124df
--- /dev/null
+++ b/lib/dnn/activation_relu.h
@@ -0,0 +1,49 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_ACTIVATION_RELU_H
+#define _GCTL_ACTIVATION_RELU_H
+
+#include "activation.h"
+
+namespace gctl
+{
+    class relu : public dnn_activation
+    {
+    public:
+        relu();
+        virtual ~relu();
+
+        void activate(const matrix<double> &z, matrix<double> &a);
+        void apply_jacobian(const matrix<double> &z, const matrix<double> &a, 
+            const matrix<double> &f, matrix<double> &g);
+        std::string activation_name() const;
+        activation_type_e activation_type() const;
+    };
+}
+
+#endif // _GCTL_ACTIVATION_RELU_H
diff --git a/lib/dnn/activation_sigmoid.cpp b/lib/dnn/activation_sigmoid.cpp
new file mode 100644
index 0000000..5731159
--- /dev/null
+++ b/lib/dnn/activation_sigmoid.cpp
@@ -0,0 +1,82 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "activation_sigmoid.h"
+
+gctl::sigmoid::sigmoid() {}
+
+gctl::sigmoid::~sigmoid() {}
+
+void gctl::sigmoid::activate(const matrix<double> &z, matrix<double> &a)
+{
+    // a = activation(z) = 1 / (1 + exp(-z))
+    // Z = [z1, ..., zn], A = [a1, ..., an], n observations
+    a.resize(z.row_size(), z.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < a.row_size(); i++)
+    {
+        for (j = 0; j < a.col_size(); j++)
+        {
+            a[i][j] = 1.0/(1.0 + exp(-1.0*z[i][j]));   
+        }
+    }
+    return;
+}
+
+void gctl::sigmoid::apply_jacobian(const matrix<double> &z, 
+    const matrix<double> &a, const matrix<double> &f, matrix<double> &g)
+{
+    // Apply the Jacobian matrix J to a vector f
+    // J = d_a / d_z = diag(a .* (1 - a))
+    // g = J * f = a .* (1 - a) .* f
+    // Z = [z1, ..., zn], G = [g1, ..., gn], F = [f1, ..., fn]
+    // Note: When entering this function, Z and G may point to the same matrix
+    g.resize(a.row_size(), a.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < g.row_size(); i++)
+    {
+        for (j = 0; j < g.col_size(); j++)
+        {
+            g[i][j] = a[i][j]*(1.0 - a[i][j])*f[i][j];   
+        }
+    }
+    return;
+}
+
+std::string gctl::sigmoid::activation_name() const
+{
+    return "Sigmoid";
+}
+
+gctl::activation_type_e gctl::sigmoid::activation_type() const
+{
+    return Sigmoid;
+}
\ No newline at end of file
diff --git a/lib/dnn/activation_sigmoid.h b/lib/dnn/activation_sigmoid.h
new file mode 100644
index 0000000..fc10235
--- /dev/null
+++ b/lib/dnn/activation_sigmoid.h
@@ -0,0 +1,49 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_ACTIVATION_SIGMOID_H
+#define _GCTL_ACTIVATION_SIGMOID_H
+
+#include "activation.h"
+
+namespace gctl
+{
+    class sigmoid : public dnn_activation
+    {
+    public:
+        sigmoid();
+        virtual ~sigmoid();
+
+        void activate(const matrix<double> &z, matrix<double> &a);
+        void apply_jacobian(const matrix<double> &z, const matrix<double> &a, 
+            const matrix<double> &f, matrix<double> &g);
+        std::string activation_name() const;
+        activation_type_e activation_type() const;
+    };
+}
+
+#endif // _GCTL_ACTIVATION_SIGMOID_H
\ No newline at end of file
diff --git a/lib/dnn/activation_softmax.cpp b/lib/dnn/activation_softmax.cpp
new file mode 100644
index 0000000..ee29bf1
--- /dev/null
+++ b/lib/dnn/activation_softmax.cpp
@@ -0,0 +1,102 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "activation_softmax.h"
+
+gctl::softmax::softmax() {}
+
+gctl::softmax::~softmax() {}
+
+void gctl::softmax::activate(const matrix<double> &z, matrix<double> &a)
+{
+    // a = activation(z) = softmax(z)
+    // Z = [z1, ..., zn], A = [a1, ..., an], n observations
+    a.resize(z.row_size(), z.col_size());
+
+    double max_z, row_s;
+    int i, j;
+#pragma omp parallel for private (i, j, max_z, row_s) schedule(guided) 
+    for (j = 0; j < a.col_size(); j++)
+    {
+        max_z = z[0][j];
+        for (i = 0; i < a.row_size(); i++)
+        {
+            max_z = max_z>z[i][j]?max_z:z[i][j];
+        }
+
+        row_s = 0.0;
+        for (i = 0; i < a.row_size(); i++)
+        {
+            row_s += std::exp(z[i][j] - max_z);
+        }
+
+        for (i = 0; i < a.row_size(); i++)
+        {
+            a[i][j] = std::exp(z[i][j] - max_z)/row_s;
+        }
+    }
+    return;
+}
+
+void gctl::softmax::apply_jacobian(const matrix<double> &z, 
+    const matrix<double> &a, const matrix<double> &f, matrix<double> &g)
+{
+    // Apply the Jacobian matrix J to a vector f
+    // J = d_a / d_z = diag(a) - a * a'
+    // g = J * f = a .* f - a * (a' * f) = a .* (f - a'f)
+    // Z = [z1, ..., zn], G = [g1, ..., gn], F = [f1, ..., fn]
+    // Note: When entering this function, Z and G may point to the same matrix
+    g.resize(a.row_size(), a.col_size());
+
+    double af;
+    int i, j;
+#pragma omp parallel for private (i, j, af) schedule(guided) 
+    for (j = 0; j < g.col_size(); j++)
+    {
+        af = 0.0;
+        for (i = 0; i < g.row_size(); i++)
+        {
+            af += a[i][j]*f[i][j];   
+        }
+
+        for (i = 0; i < g.row_size(); i++)
+        {
+            g[i][j] = a[i][j]*(f[i][j] - af);
+        }
+    }
+    return;
+}
+
+std::string gctl::softmax::activation_name() const
+{
+    return "SoftMax";
+}
+
+gctl::activation_type_e gctl::softmax::activation_type() const
+{
+    return SoftMax;
+}
\ No newline at end of file
diff --git a/lib/dnn/activation_softmax.h b/lib/dnn/activation_softmax.h
new file mode 100644
index 0000000..53399d3
--- /dev/null
+++ b/lib/dnn/activation_softmax.h
@@ -0,0 +1,49 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_ACTIVATION_SOFTMAX_H
+#define _GCTL_ACTIVATION_SOFTMAX_H
+
+#include "activation.h"
+
+namespace gctl
+{
+    class softmax : public dnn_activation
+    {
+    public:
+        softmax();
+        virtual ~softmax();
+
+        void activate(const matrix<double> &z, matrix<double> &a);
+        void apply_jacobian(const matrix<double> &z, const matrix<double> &a, 
+            const matrix<double> &f, matrix<double> &g);
+        std::string activation_name() const;
+        activation_type_e activation_type() const;
+    };
+}
+
+#endif // _GCTL_ACTIVATION_SOFTMAX_H
\ No newline at end of file
diff --git a/lib/dnn/activation_tanh.cpp b/lib/dnn/activation_tanh.cpp
new file mode 100644
index 0000000..5097a83
--- /dev/null
+++ b/lib/dnn/activation_tanh.cpp
@@ -0,0 +1,83 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "activation_tanh.h"
+
+gctl::tanh::tanh() {}
+
+gctl::tanh::~tanh() {}
+
+void gctl::tanh::activate(const matrix<double> &z, matrix<double> &a)
+{
+    // a = activation(z) = tanh(z)
+    // Z = [z1, ..., zn], A = [a1, ..., an], n observations
+    a.resize(z.row_size(), z.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < a.row_size(); i++)
+    {
+        for (j = 0; j < a.col_size(); j++)
+        {
+            a[i][j] = std::tanh(z[i][j]);   
+        }
+    }
+    return;
+}
+
+void gctl::tanh::apply_jacobian(const matrix<double> &z, 
+    const matrix<double> &a, const matrix<double> &f, matrix<double> &g)
+{
+    // Apply the Jacobian matrix J to a vector f
+    // tanh'(x) = 1 - tanh(x)^2
+    // J = d_a / d_z = diag(1 - a^2)
+    // g = J * f = (1 - a^2) .* f
+    // Z = [z1, ..., zn], G = [g1, ..., gn], F = [f1, ..., fn]
+    // Note: When entering this function, Z and G may point to the same matrix
+    g.resize(a.row_size(), a.col_size());
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided) 
+    for (i = 0; i < g.row_size(); i++)
+    {
+        for (j = 0; j < g.col_size(); j++)
+        {
+            g[i][j] = (1.0 - a[i][j]*a[i][j])*f[i][j];   
+        }
+    }
+    return;
+}
+
+std::string gctl::tanh::activation_name() const
+{
+    return "Tanh";
+}
+
+gctl::activation_type_e gctl::tanh::activation_type() const
+{
+    return Tanh;
+}
\ No newline at end of file
diff --git a/lib/dnn/activation_tanh.h b/lib/dnn/activation_tanh.h
new file mode 100644
index 0000000..d8c02ff
--- /dev/null
+++ b/lib/dnn/activation_tanh.h
@@ -0,0 +1,49 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_ACTIVATION_TANH_H
+#define _GCTL_ACTIVATION_TANH_H
+
+#include "activation.h"
+
+namespace gctl
+{
+    class tanh : public dnn_activation
+    {
+    public:
+        tanh();
+        virtual ~tanh();
+
+        void activate(const matrix<double> &z, matrix<double> &a);
+        void apply_jacobian(const matrix<double> &z, const matrix<double> &a, 
+            const matrix<double> &f, matrix<double> &g);
+        std::string activation_name() const;
+        activation_type_e activation_type() const;
+    };
+}
+
+#endif // _GCTL_ACTIVATION_TANH_H
\ No newline at end of file
diff --git a/lib/dnn/config.h b/lib/dnn/config.h
new file mode 100644
index 0000000..81867f0
--- /dev/null
+++ b/lib/dnn/config.h
@@ -0,0 +1 @@
+/* #undef GCTL_AI_OPENMP */
diff --git a/lib/dnn/dnn.cpp b/lib/dnn/dnn.cpp
new file mode 100644
index 0000000..6ea489d
--- /dev/null
+++ b/lib/dnn/dnn.cpp
@@ -0,0 +1,542 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "dnn.h"
+
+gctl::dnn::dnn(std::string info)
+{
+    hlayer_size_ = 0; 
+    hind_layers_.reserve(100);
+    output_layer_ = nullptr;
+    wd_st_ = 0; wd_size_ = 0;
+    batch_iter_ = 0;
+    h_ss = o_ss = t_ss = 0;
+    status_ = NoneSet;
+    info_ = info;
+    has_stats_ = false;
+}
+
+gctl::dnn::~dnn()
+{
+    for (int i = 0; i < hind_layers_.size(); i++)
+    {
+        if (hind_layers_[i] != nullptr)
+        {
+            delete hind_layers_[i];
+        }
+    }
+    
+    if (output_layer_ != nullptr) delete output_layer_;
+}
+
+void gctl::dnn::load_network(std::string filename)
+{
+    std::ifstream ifile;
+    open_infile(ifile, filename, ".gctl.dnn", std::ios::in|std::ios::binary);
+
+    int h_size, stats_int = 0;
+    hlayer_type_e l_type;
+    dnn_hlayer *hlayer_ptr = nullptr;
+
+    ifile.read((char*)&h_size, sizeof(int));
+    ifile.read((char*)&stats_int, sizeof(int));
+    if (stats_int) has_stats_ = true;
+
+    for (int i = 0; i < h_size; i++)
+    {
+        ifile.read((char*)&l_type, sizeof(hlayer_type_e));
+        if (l_type == FullyConnected) hlayer_ptr = new fully_connected;
+        if (l_type == MaxPooling) hlayer_ptr = new maxpooling;
+        if (l_type == AvgPooling) hlayer_ptr = new avgpooling;
+        if (l_type == Convolution) hlayer_ptr = new convolution;
+
+        hind_layers_.push_back(hlayer_ptr);
+        hlayer_size_++;
+    }
+
+    for (int i = 0; i < h_size; i++)
+    {
+        hind_layers_[i]->load_layer_setup(ifile);
+    }
+
+    olayer_type_e o_type;
+    ifile.read((char*)&o_type, sizeof(olayer_type_e));
+    add_output_layer(o_type);
+
+    ifile.read((char*)&wd_size_, sizeof(int));
+    weights_.resize(wd_size_);
+    ders_.resize(wd_size_);
+    
+    ifile.read((char*)weights_.get(), wd_size_*sizeof(double));
+    if (has_stats_)
+    {
+        weights_mean_.resize(wd_size_);
+        weights_std_.resize(wd_size_);
+        ifile.read((char*)weights_mean_.get(), wd_size_*sizeof(double));
+        ifile.read((char*)weights_std_.get(), wd_size_*sizeof(double));
+    }
+    ifile.close();
+
+    status_ = Trained;
+    return;
+}
+
+void gctl::dnn::save_network(std::string filename) const
+{
+    if (status_ != Trained)
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not trained.");
+    }
+
+    std::ofstream ofile;
+    open_outfile(ofile, filename, ".gctl.dnn", std::ios::out|std::ios::binary);
+
+    // save hind layers' setup
+    int stats_int = 0;
+    if (has_stats_) stats_int = 1;
+    ofile.write((char*)&hlayer_size_, sizeof(int));
+    ofile.write((char*)&stats_int, sizeof(int));
+
+    hlayer_type_e l_type;
+    for (int i = 0; i < hlayer_size_; i++)
+    {
+        l_type = hind_layers_[i]->get_layer_type();
+        ofile.write((char*)&l_type, sizeof(hlayer_type_e));
+    }
+
+    for (int i = 0; i < hlayer_size_; i++)
+    {
+        hind_layers_[i]->save_layer_setup(ofile);
+    }
+
+    // save output layer setup
+    olayer_type_e o_type = output_layer_->get_output_type();
+    ofile.write((char*)&o_type, sizeof(olayer_type_e));
+
+    // save hind layers' parameters
+    ofile.write((char*)&wd_size_, sizeof(int));
+    ofile.write((char*)weights_.get(), wd_size_*sizeof(double));
+
+    if (has_stats_)
+    {
+        ofile.write((char*)weights_mean_.get(), wd_size_*sizeof(double));
+        ofile.write((char*)weights_std_.get(),  wd_size_*sizeof(double));
+    }
+    ofile.close();
+    return;
+}
+
+void gctl::dnn::save_layer2text(int l_idx, std::string filename) const
+{
+    if (hind_layers_.empty())
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not setup properly.");
+    }
+
+    if (status_ != Trained)
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not trained.");
+    }
+
+    if (l_idx < 0 || l_idx >= hlayer_size_)
+    {
+        throw std::invalid_argument("[gctl::dnn] Invalid index of the hinden layers.");
+    }
+
+    std::ofstream ofile;
+    open_outfile(ofile, filename + "_weight", ".txt", std::ios::out);
+    hind_layers_[l_idx]->save_weights2text(weights_, ofile);
+    ofile.close();
+
+    if (has_stats_)
+    {
+        open_outfile(ofile, filename + "_weight_mean", ".txt", std::ios::out);
+        hind_layers_[l_idx]->save_weights2text(weights_mean_, ofile);
+        ofile.close();
+
+        open_outfile(ofile, filename + "_weight_std", ".txt", std::ios::out);
+        hind_layers_[l_idx]->save_weights2text(weights_std_, ofile);
+        ofile.close();
+    }
+    return;
+}
+
+void gctl::dnn::show_network()
+{
+    if (output_layer_ == nullptr || hind_layers_.empty())
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not setup properly.");
+    }
+
+    std::clog << "=============================\n";
+    std::clog << "GCTL's DNN Setup Panel\n";
+    std::clog << "-----------------------------\n";
+    std::clog << "Network Info: " << info_ << "\n";
+    std::clog << "Observation Number: " << obs_num_;
+    if (batch_size_ > 0) std::clog << ", Batch Size: " << batch_size_ << "\n";
+    else std::clog << ", Batch Size: Not set\n";
+    std::clog << "Output Layer: " << output_layer_->get_output_name() << "\n";
+    std::clog << "Hind Layer's Number: " << hlayer_size_ << "\n";
+
+    for (int i = 0; i < hlayer_size_; i++)
+    {
+        std::clog << "Layer-" << i+1 << ": " << hind_layers_[i]->layer_info() << "\n";
+    }
+    return;
+}
+
+void gctl::dnn::add_hind_layer(int in_s, int out_s, hlayer_type_e h_type, activation_type_e a_type)
+{
+    dnn_hlayer *hlayer_ptr = nullptr;
+    if (h_type == FullyConnected)
+    {
+        // 初始化新的隐藏层为全链接层
+        hlayer_ptr = new fully_connected(wd_st_, in_s, out_s, a_type);
+        hind_layers_.push_back(hlayer_ptr);
+
+        wd_size_ += out_s*(in_s+1);
+        wd_st_ = wd_size_;
+        hlayer_size_++;
+    }
+    else throw std::invalid_argument("[gctl::dnn] Invalid initiate parameters for current layer type.");
+
+    h_ss = 1;
+    if ((h_ss + o_ss + t_ss) == 3) status_ = AllSet;
+    else status_ = PartSet;
+    return;
+}
+
+void gctl::dnn::add_hind_layer(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, 
+    int stride_cols, hlayer_type_e h_type, pad_type_e p_type, activation_type_e acti_type)
+{
+    dnn_hlayer *hlayer_ptr = nullptr;
+    if (h_type == MaxPooling)
+    {
+        hlayer_ptr = new maxpooling(in_rows, in_cols, pool_rows, pool_cols, stride_rows, stride_cols, p_type, acti_type);
+        hind_layers_.push_back(hlayer_ptr);
+
+        hlayer_size_++;
+    }
+    else if (h_type == AvgPooling)
+    {
+        hlayer_ptr = new avgpooling(in_rows, in_cols, pool_rows, pool_cols, stride_rows, stride_cols, p_type, acti_type);
+        hind_layers_.push_back(hlayer_ptr);
+
+        hlayer_size_++;
+    }
+    else throw std::invalid_argument("[gctl::dnn] Invalid initiate parameters for current layer type.");
+
+    h_ss = 1;
+    if ((h_ss + o_ss + t_ss) == 3) status_ = AllSet;
+    else status_ = PartSet;
+    return;
+}
+
+void gctl::dnn::add_hind_layer(int channels, int in_rows, int in_cols, int filter_rows, int filter_cols, 
+    int stride_rows, int stride_cols, hlayer_type_e h_type, pad_type_e p_type, activation_type_e acti_type)
+{
+    dnn_hlayer *hlayer_ptr = nullptr;
+    if (h_type == Convolution)
+    {
+        hlayer_ptr = new convolution(wd_st_, channels, in_rows, in_cols, filter_rows, filter_cols, stride_rows, stride_cols, p_type, acti_type);
+        hind_layers_.push_back(hlayer_ptr);
+
+        wd_size_ += filter_rows*filter_cols*channels + 1;
+        wd_st_ = wd_size_;
+        hlayer_size_++;
+    }
+    else throw std::invalid_argument("[gctl::dnn] Invalid initiate parameters for current layer type.");
+
+    h_ss = 1;
+    if ((h_ss + o_ss + t_ss) == 3) status_ = AllSet;
+    else status_ = PartSet;
+    return;
+}
+
+void gctl::dnn::add_output_layer(olayer_type_e o_type)
+{
+    if (o_type == RegressionMSE) output_layer_ = new rmse;
+    else if (o_type == MultiClassEntropy) output_layer_ = new multi_entropy;
+    else if (o_type == BinaryClassEntropy) output_layer_ = new binary_entropy;
+    else throw std::invalid_argument("[gctl::dnn] Invalid output layer type.");
+
+    o_ss = 1;
+    if ((h_ss + o_ss + t_ss) == 3) status_ = AllSet;
+    else status_ = PartSet;
+    return;
+}
+
+void gctl::dnn::add_train_set(const matrix<double> &train_obs, const matrix<double> &train_tar, int batch_size)
+{
+    if (batch_size < 0)
+    {
+        throw std::invalid_argument("[gctl::dnn] Invalid batch size.");
+    }
+    batch_size_ = batch_size;
+
+    obs_num_ = train_obs.col_size();
+    if(obs_num_ != train_tar.col_size())
+    {
+        throw std::invalid_argument("[gctl::dnn] Observation sizes do not match.");
+    }
+
+    if (output_layer_ == nullptr)
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not setup properly.");
+    }
+    else output_layer_->check_target_data(train_tar);
+
+    if (batch_size_ == 0 || batch_size_ >= obs_num_)
+    {
+        obs_.resize(1);
+        targets_.resize(1);
+        obs_[0] = train_obs;
+        targets_[0] = train_tar;
+
+        batch_num_ = 1;
+
+        t_ss = 1;
+        if ((h_ss + o_ss + t_ss) == 3) status_ = AllSet;
+        else status_ = PartSet;
+        return;
+    }
+
+    batch_num_ = ceil(1.0*obs_num_/batch_size_);
+
+    obs_.resize(batch_num_);
+    targets_.resize(batch_num_);
+
+    for (size_t i = 0; i < batch_num_-1; i++)
+    {
+        obs_[i].resize(train_obs.row_size(), batch_size);
+        targets_[i].resize(train_tar.row_size(), batch_size);
+
+        for (size_t c = 0; c < batch_size; c++)
+        {
+            for (size_t r = 0; r < train_obs.row_size(); r++)
+            {
+                obs_[i][r][c] = train_obs[r][i*batch_size + c];
+            }
+
+            for (size_t r = 0; r < train_tar.row_size(); r++)
+            {
+                targets_[i][r][c] = train_tar[r][i*batch_size + c];
+            }
+        }
+    }
+
+    obs_[batch_num_-1].resize(train_obs.row_size(), obs_num_ - (batch_num_-1)*batch_size);
+    targets_[batch_num_-1].resize(train_tar.row_size(), obs_num_ - (batch_num_-1)*batch_size);
+
+    for (size_t c = 0; c < obs_num_ - (batch_num_-1)*batch_size; c++)
+    {
+        for (size_t r = 0; r < train_obs.row_size(); r++)
+        {
+            obs_[batch_num_-1][r][c] = train_obs[r][(batch_num_-1)*batch_size + c];
+        }
+
+        for (size_t r = 0; r < train_tar.row_size(); r++)
+        {
+            targets_[batch_num_-1][r][c] = train_tar[r][(batch_num_-1)*batch_size + c];
+        }
+    }
+
+    t_ss = 1;
+    if ((h_ss + o_ss + t_ss) == 3) status_ = AllSet;
+    else status_ = PartSet;
+    return;
+}
+
+void gctl::dnn::init_network(double mu, double sigma, unsigned int seed)
+{
+    if (status_ != AllSet)
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not setup properly.");
+    }
+    
+    if (hind_layers_.size() > 1)
+    {
+        for (int i = 1; i < hind_layers_.size(); i++)
+        {
+            if (hind_layers_[i]->in_size() != hind_layers_[i-1]->out_size())
+            {
+                throw std::invalid_argument("[gctl::dnn] Layer sizes do not match.");
+            }
+        }
+    }
+
+    weights_.resize(wd_size_);
+    ders_.resize(wd_size_);
+ 
+    gctl::random(weights_, mu, sigma, RdNormal, seed);
+    gctl::random(ders_, mu, sigma, RdNormal, seed);
+
+    status_ = Initialized;
+    return;
+}
+
+void gctl::dnn::train_network(const sgd_para &pa, sgd_solver_type solver)
+{
+    if (status_ != Initialized)
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not initialized.");
+    }
+
+    set_sgd_para(pa);
+    SGD_Minimize(weights_, solver);
+
+    status_ = Trained;
+    return;
+}
+
+void gctl::dnn::train_network(const lgd_para &pa)
+{
+    if (status_ != Initialized)
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not initialized.");
+    }
+
+    has_stats_ = true;
+
+    set_lgd_para(pa);
+    LGD_Minimize(weights_, weights_mean_, weights_std_);
+
+    status_ = Trained;
+    return;
+}
+
+void gctl::dnn::predict(const matrix<double> &pre_obs, matrix<double> &predicts)
+{
+    if (status_ != Trained)
+    {
+        throw std::invalid_argument("[gctl::dnn] The DNN is not trained.");
+    }
+
+    forward_propagation(pre_obs);
+    predicts = hind_layers_.back()->forward_propagation_data();
+    return;
+}
+
+void gctl::dnn::forward_propagation(const matrix<double> &input)
+{
+    // First layer
+    if (input.row_size() != hind_layers_[0]->in_size())
+    {
+        throw std::invalid_argument("[gctl::dnn] Input data have incorrect dimension.");
+    }
+
+    hind_layers_[0]->forward_propagation(weights_, input);
+    // The following layers
+    for (int i = 1; i < hlayer_size_; i++)
+    {
+        hind_layers_[i]->forward_propagation(weights_, hind_layers_[i-1]->forward_propagation_data());
+    }
+    return;
+}
+
+void gctl::dnn::backward_propagation(const matrix<double> &input, const matrix<double> &target)
+{
+    dnn_hlayer *first_layer = hind_layers_.front();
+    dnn_hlayer *last_layer = hind_layers_.back();
+
+    // Let output layer compute back-propagation data
+    output_layer_->evaluation(last_layer->forward_propagation_data(), target);
+
+    // If there is only one hidden layer, "prev_layer_data" will be the input data
+    if (hlayer_size_ == 1)
+    {
+        first_layer->backward_propagation(weights_, ders_, input, output_layer_->backward_propagation_data());
+        return;
+    }
+    
+    // Compute gradients for the last hidden layer
+    last_layer->backward_propagation(
+        weights_, ders_,
+        hind_layers_[hlayer_size_-2]->forward_propagation_data(), 
+        output_layer_->backward_propagation_data());
+
+    // Compute gradients for all the hidden layers except for the first one and the last one
+    for (int i = hlayer_size_-2; i > 0; i--)
+    {
+        hind_layers_[i]->backward_propagation(
+            weights_, ders_,
+            hind_layers_[i-1]->forward_propagation_data(),
+            hind_layers_[i+1]->backward_propagation_data());
+    }
+
+    // Compute gradients for the first layer
+    first_layer->backward_propagation(weights_, ders_, input, hind_layers_[1]->backward_propagation_data());
+    return;
+}
+
+double gctl::dnn::SGD_Evaluate(const array<double> &x, array<double> &g)
+{
+    // run the forward and backward process once
+    forward_propagation(obs_[batch_iter_%batch_num_]);
+    backward_propagation(obs_[batch_iter_%batch_num_], targets_[batch_iter_%batch_num_]);
+    batch_iter_++;
+    // get the network's derivatives
+    g = ders_;
+    return output_layer_->loss_value();
+}
+
+int gctl::dnn::SGD_Progress(double fx, const array<double> &x, const sgd_para &param, const int k)
+{
+    // 清屏
+	winsize win;
+	ioctl(0, TIOCGWINSZ, &win);
+	std::clog << "\033[2J" << "\033[" << win.ws_row << "A";
+
+    show_network();
+    sgd_solver::show_solver();
+    std::clog << GCTL_CLEARLINE << "\rF(x) = " << fx << ", Train-Times = " << k << "\n";
+    return 0;
+}
+
+double gctl::dnn::LGD_Evaluate(const array<double> &x, array<double> &g)
+{
+    // run the forward and backward process once
+    forward_propagation(obs_[batch_iter_%batch_num_]);
+    backward_propagation(obs_[batch_iter_%batch_num_], targets_[batch_iter_%batch_num_]);
+    batch_iter_++;
+    // get the network's derivatives
+    g = ders_;
+    return output_layer_->loss_value();
+}
+
+int gctl::dnn::LGD_Progress(const int curr_t, const double curr_fx, const double mean_fx, const double best_fx, const lgd_para &param)
+{
+    // 清屏
+	winsize win;
+	ioctl(0, TIOCGWINSZ, &win);
+	std::clog << "\033[2J" << "\033[" << win.ws_row << "A";
+
+    show_network();
+    lgd_solver::show_solver();
+    std::clog << "F(x) = " << curr_fx << ", Mean F(x) = " << mean_fx <<  ", Best F(x) = " << best_fx << ", Times = " << curr_t << "\n";
+    return 0;
+}
\ No newline at end of file
diff --git a/lib/dnn/dnn.h b/lib/dnn/dnn.h
new file mode 100644
index 0000000..69308d2
--- /dev/null
+++ b/lib/dnn/dnn.h
@@ -0,0 +1,138 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_DNN_H
+#define _GCTL_DNN_DNN_H
+
+#include "hlayer_fully_connected.h"
+#include "hlayer_maxpooling.h"
+#include "hlayer_avgpooling.h"
+#include "hlayer_convolution.h"
+
+#include "olayer_rmse.h"
+#include "olayer_binaryentropy.h"
+#include "olayer_multientropy.h"
+
+#include "gctl/optimization.h"
+
+namespace gctl
+{
+    enum dnn_status_e
+    {
+        NoneSet,
+        PartSet,
+        AllSet,
+        Initialized,
+        Trained,
+    };
+
+    class dnn : public sgd_solver, public lgd_solver
+    {
+    public:
+        dnn(std::string info);
+        virtual ~dnn();
+
+        /**
+         * @brief Load DNN network from file.
+         * 
+         * @param filename Name of the saved binary file.
+         */
+        void load_network(std::string filename);
+
+        /**
+         * @brief Save DNN network to file.
+         * 
+         * @param filename Name of the output binary file.
+         */
+        void save_network(std::string filename) const;
+
+        /**
+         * @brief Save DNN's layer to a text file.
+         * 
+         * @param l_idx    Layer index.
+         * @param filename Name of the output file.
+         */
+        void save_layer2text(int l_idx, std::string filename) const;
+
+        /**
+         * @brief Display DNN's setup on screen.
+         * 
+         */
+        void show_network();
+
+        /**
+         * @brief Add a fully connected hind layer
+         * 
+         * @param in_s Input parameter size.
+         * @param out_s Output parameter size.
+         * @param h_type Layer type. For now it must be FullyConnected.
+         * @param a_type Output type.
+         */
+        void add_hind_layer(int in_s, int out_s, hlayer_type_e h_type, activation_type_e a_type);
+        void add_hind_layer(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, 
+            int stride_cols, hlayer_type_e h_type, pad_type_e p_type, activation_type_e acti_type);
+        void add_hind_layer(int channels, int in_rows, int in_cols, int filter_rows, int filter_cols, 
+            int stride_rows, int stride_cols, hlayer_type_e h_type, pad_type_e p_type, activation_type_e acti_type);
+        void add_output_layer(olayer_type_e o_type);
+        void add_train_set(const matrix<double> &train_obs, const matrix<double> &train_tar, int batch_size = 0);
+        void init_network(double mu, double sigma, unsigned int seed = 0);
+        
+        void train_network(const sgd_para &pa, sgd_solver_type solver = ADAM);
+        void train_network(const lgd_para &pa);
+        void predict(const matrix<double> &pre_obs, matrix<double> &predicts);
+
+        double SGD_Evaluate(const array<double> &x, array<double> &g);
+        int SGD_Progress(double fx, const array<double> &x, const sgd_para &param, const int k);
+
+        double LGD_Evaluate(const array<double> &x, array<double> &g);
+        int LGD_Progress(const int curr_t, const double curr_fx, const double mean_fx, const double best_fx, const lgd_para &param);
+
+    private:
+        void forward_propagation(const matrix<double> &input);
+        void backward_propagation(const matrix<double> &input, const matrix<double> &target);
+
+    private:
+        int hlayer_size_;
+        std::vector<dnn_hlayer*> hind_layers_;
+        dnn_olayer *output_layer_;
+
+        bool has_stats_;
+        int wd_st_, wd_size_;
+        array<double> weights_, ders_;
+        array<double> weights_mean_, weights_std_;
+
+        int batch_iter_, batch_size_, batch_num_, obs_num_;
+        array<matrix<double> > obs_, targets_;
+
+        int h_ss, o_ss, t_ss;
+        dnn_status_e status_;
+
+        std::string info_;
+    };
+}
+
+#endif // _GCTL_DNN_DNN_H
\ No newline at end of file
diff --git a/lib/dnn/hlayer.cpp b/lib/dnn/hlayer.cpp
new file mode 100644
index 0000000..2dd9b5d
--- /dev/null
+++ b/lib/dnn/hlayer.cpp
@@ -0,0 +1,120 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "hlayer.h"
+
+gctl::dnn_hlayer::dnn_hlayer(){}
+
+gctl::dnn_hlayer::~dnn_hlayer()
+{
+    delete activator_;
+}
+
+int gctl::dnn_hlayer::in_size() const
+{
+    return w_is_;
+}
+
+int gctl::dnn_hlayer::out_size() const
+{
+    return w_outs_;
+}
+
+int gctl::dnn_hlayer::obs_size() const
+{
+    return o_is_;
+}
+
+const gctl::matrix<double> &gctl::dnn_hlayer::forward_propagation_data()
+{
+    return a_;
+}
+
+const gctl::matrix<double> &gctl::dnn_hlayer::backward_propagation_data()
+{
+    return der_in_;
+}
+
+const gctl::matrix<double> &gctl::dnn_hlayer::get_linear_term()
+{
+    return z_;
+}
+
+const gctl::matrix<double> &gctl::dnn_hlayer::get_linear_gradient()
+{
+    return der_z_;
+}
+
+gctl::activation_type_e gctl::dnn_hlayer::get_activation_type() const
+{
+    return activator_->activation_type();
+}
+
+std::string gctl::dnn_hlayer::get_activation_name() const
+{
+    return activator_->activation_name();
+}
+
+void gctl::dnn_hlayer::cal_valid_padding_idx(array<int> &idx, int i, int j, int i_rows, int i_cols, int p_rows, int p_cols, int s_rows, int s_cols)
+{
+    idx.resize(p_rows*p_cols);
+    idx.assign_all(-1);
+
+    int row_st = i*s_rows;
+    int col_st = j*s_cols;
+
+    for (size_t r = 0; r < p_rows; r++)
+    {
+        for (size_t c = 0; c < p_cols; c++)
+        {
+            idx[c + r*p_cols] = (row_st + r)*i_cols + col_st + c;
+        }
+    }
+    return;
+}
+
+void gctl::dnn_hlayer::cal_same_padding_idx(array<int> &idx, int i, int j, int i_rows, int i_cols, int p_rows, int p_cols, int s_rows, int s_cols, int u_pad, int l_pad)
+{
+    idx.resize(p_rows*p_cols);
+    idx.assign_all(-1);
+
+    int row_st = i*s_rows - u_pad;
+    int col_st = j*s_cols - l_pad;
+
+    for (size_t r = 0; r < p_rows; r++)
+    {
+        for (size_t c = 0; c < p_cols; c++)
+        {
+            if (row_st + r >= 0 && row_st + r < i_rows && 
+                col_st + c >= 0 && col_st + c < i_cols)
+            {
+                idx[c + r*p_cols] = (row_st + r)*i_cols + col_st + c;
+            }
+        }
+    }
+    return;
+}
\ No newline at end of file
diff --git a/lib/dnn/hlayer.h b/lib/dnn/hlayer.h
new file mode 100644
index 0000000..29edd8d
--- /dev/null
+++ b/lib/dnn/hlayer.h
@@ -0,0 +1,99 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_HLAYER_H
+#define _GCTL_DNN_HLAYER_H
+
+#include "activation_identity.h"
+#include "activation_mish.h"
+#include "activation_relu.h"
+#include "activation_prelu.h"
+#include "activation_sigmoid.h"
+#include "activation_softmax.h"
+#include "activation_tanh.h"
+
+namespace gctl
+{
+    enum hlayer_type_e
+    {
+        FullyConnected,
+        MaxPooling,
+        AvgPooling,
+        Convolution,
+    };
+
+    enum pad_type_e
+    {
+        Valid,
+        Same,
+    };
+
+    class dnn_hlayer
+    {
+    public:
+        dnn_hlayer();
+        virtual ~dnn_hlayer();
+
+        int in_size() const;
+        int out_size() const;
+        int obs_size() const;
+
+        const matrix<double> &forward_propagation_data();
+        const matrix<double> &backward_propagation_data();
+        const matrix<double> &get_linear_term();
+        const matrix<double> &get_linear_gradient();
+
+        activation_type_e get_activation_type() const;
+        std::string get_activation_name() const;
+
+        virtual void forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data) = 0;
+        virtual void backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+            const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data) = 0;
+
+        virtual hlayer_type_e get_layer_type() const = 0;
+        virtual std::string get_layer_name() const = 0;
+        virtual std::string layer_info() const = 0;
+        virtual void save_layer_setup(std::ofstream &os) const = 0;
+        virtual void load_layer_setup(std::ifstream &is) = 0;
+        virtual void save_weights2text(const array<double> &all_weights, std::ofstream &os) const = 0;
+
+    protected:
+        void cal_valid_padding_idx(array<int> &idx, int i, int j, int i_rows, int i_cols, int p_rows, int p_cols, int s_rows, int s_cols);
+        void cal_same_padding_idx(array<int> &idx, int i, int j, int i_rows, int i_cols, int p_rows, int p_cols, int s_rows, int s_cols, int u_pad, int l_pad);
+    
+    protected:
+        dnn_activation *activator_; // activation object
+        int w_is_, w_outs_; // weight start index, weight input size, weight output size
+        int o_is_; // observation input size
+        matrix<double> z_;           // Linear term, z = W' * in + b
+        matrix<double> a_;           // Output of this layer, a_ = act(z)
+        matrix<double> der_in_;      // Derivative of the input of this layer, which is also the output of previous layer
+        matrix<double> der_z_;       // Derivative of the linear term
+    };
+}
+
+#endif // _GCTL_DNN_HLAYER_H
\ No newline at end of file
diff --git a/lib/dnn/hlayer_avgpooling.cpp b/lib/dnn/hlayer_avgpooling.cpp
new file mode 100644
index 0000000..f7dfff6
--- /dev/null
+++ b/lib/dnn/hlayer_avgpooling.cpp
@@ -0,0 +1,220 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "hlayer_avgpooling.h"
+
+gctl::avgpooling::avgpooling(){}
+
+gctl::avgpooling::avgpooling(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, int stride_cols, 
+    pad_type_e pl_type, activation_type_e acti_type)
+{
+    init_avgpooling(in_rows, in_cols, pool_rows, pool_cols, stride_rows, stride_cols, pl_type, acti_type);
+}
+
+gctl::avgpooling::~avgpooling(){}
+
+void gctl::avgpooling::init_avgpooling(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, 
+    int stride_cols, pad_type_e pl_type, activation_type_e acti_type)
+{
+    i_rows_ = in_rows; i_cols_ = in_cols;
+    p_rows_ = pool_rows; p_cols_ = pool_cols;
+    s_rows_ = stride_rows; s_cols_ = stride_cols;
+    p_type_ = pl_type;
+
+    if (acti_type == Identity) activator_ = new identity;
+    else if (acti_type == Mish) activator_ = new mish;
+    else if (acti_type == ReLU) activator_ = new relu;
+    else if (acti_type == PReLU) activator_ = new prelu;
+    else if (acti_type == Sigmoid) activator_ = new sigmoid;
+    else if (acti_type == SoftMax) activator_ = new softmax;
+    else if (acti_type == Tanh) activator_ = new tanh;
+    else throw std::invalid_argument("[gctl::avgpooling] Invalid activation type.");
+
+    // 计算输出大小
+    o_rows_ = (in_rows + s_rows_ - p_rows_)/s_rows_;
+    o_cols_ = (in_cols + s_cols_ - p_cols_)/s_cols_;
+    if (pl_type == Same && (in_rows + s_rows_ - p_rows_)%s_rows_ != 0)
+    {
+        o_rows_++;
+        u_pad_ = ((in_rows + s_rows_ - p_rows_)%s_rows_)/2;
+    }
+    else u_pad_ = 0;
+
+    if (pl_type == Same && (in_cols + s_cols_ - p_cols_)%s_cols_ != 0)
+    {
+        o_cols_++;
+        l_pad_ = ((in_cols + s_cols_ - p_cols_)%s_cols_)/2;
+    }
+    else l_pad_ = 0;
+
+    w_is_ = i_rows_*i_cols_;
+    w_outs_ = o_rows_*o_cols_;
+    p_idx_.resize(p_rows_*p_cols_);
+    return;
+}
+
+void gctl::avgpooling::forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data)
+{
+    // z_: out_size x nobs
+    // a_: out_size x nobs
+    o_is_ = prev_layer_data.col_size();
+    // Forward linear terms
+    z_.resize(w_outs_, o_is_);
+    a_.resize(w_outs_, o_is_);
+
+    double p_val;
+    for (size_t i = 0; i < o_rows_; i++)
+    {
+        for (size_t j = 0; j < o_cols_; j++)
+        {
+            if (p_type_ == Valid) cal_valid_padding_idx(p_idx_, i, j, i_rows_, i_cols_, p_rows_, p_cols_, s_rows_, s_cols_);
+            if (p_type_ == Same) cal_same_padding_idx(p_idx_, i, j, i_rows_, i_cols_, p_rows_, p_cols_, s_rows_, s_cols_, u_pad_, l_pad_);
+
+            for (size_t o = 0; o < o_is_; o++)
+            {
+                p_val = 0.0;
+                for (size_t p = 0; p < p_idx_.size(); p++)
+                {
+                    if (p_idx_[p] != -1) p_val += prev_layer_data[p_idx_[p]][o];
+                }
+
+                z_[i*o_cols_ + j][o] = p_val/p_idx_.size();
+            }
+        }
+    }
+
+    // Apply activation function
+    activator_->activate(z_, a_);
+    return;
+}
+
+void gctl::avgpooling::backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+    const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data)
+{
+    der_z_.resize(w_outs_, o_is_);
+    der_in_.resize(w_is_, o_is_, 0.0);
+    // prev_layer_data: in_size x nobs
+    // next_layer_data: out_size x nobs
+    // After forward stage, m_z contains z = max_avgpooling(in)
+    // Now we need to calculate d(L) / d(z) = [d(a) / d(z)] * [d(L) / d(a)]
+    // d(L) / d(z) is computed in the next layer, contained in next_layer_data
+    // The Jacobian matrix J = d(a) / d(z) is determined by the activation function
+    // der_z_: out_size x nobs   
+    activator_->apply_jacobian(z_, a_, next_layer_data, der_z_);
+
+    for (size_t j = 0; j < o_is_; j++)
+    {
+        for (size_t r = 0; r < o_rows_; r++)
+        {
+            for (size_t c = 0; c < o_cols_; c++)
+            {
+                if (p_type_ == Valid) cal_valid_padding_idx(p_idx_, r, c, i_rows_, i_cols_, p_rows_, p_cols_, s_rows_, s_cols_);
+                if (p_type_ == Same) cal_same_padding_idx(p_idx_, r, c, i_rows_, i_cols_, p_rows_, p_cols_, s_rows_, s_cols_, u_pad_, l_pad_);
+
+                for (size_t p = 0; p < p_rows_*p_cols_; p++)
+                {
+                    if (p_idx_[p] != -1) der_in_[p_idx_[p]][j] += der_z_[r*o_cols_ + c][j]/(p_rows_*p_cols_);
+                }
+            }
+        }
+    }
+    return;
+}
+
+gctl::hlayer_type_e gctl::avgpooling::get_layer_type() const
+{
+    return AvgPooling;
+}
+
+std::string gctl::avgpooling::get_layer_name() const
+{
+    return "AvgPooling";
+}
+
+std::string gctl::avgpooling::layer_info() const
+{
+    std::string info = std::to_string(w_is_) + "x" + std::to_string(w_outs_)
+        + " (" + std::to_string(i_rows_) + "," + std::to_string(i_cols_) + ")x"
+        + "(" + std::to_string(o_rows_) + "," + std::to_string(o_cols_) + "), "
+        + std::to_string(p_rows_) + "x" + std::to_string(p_cols_)
+        + ", " + std::to_string(s_rows_) + "x" + std::to_string(s_cols_)
+        + ", AvgPooling ";
+
+    if (p_type_ == Same) info += "(Same), ";
+    if (p_type_ == Valid) info += "(Valid), ";
+
+    info += activator_->activation_name();
+    return info;
+}
+
+void gctl::avgpooling::save_layer_setup(std::ofstream &os) const
+{
+    hlayer_type_e h_type = get_layer_type();
+    activation_type_e a_type = get_activation_type();
+
+    os.write((char*)&i_rows_, sizeof(int));
+    os.write((char*)&i_cols_, sizeof(int));
+    os.write((char*)&p_rows_, sizeof(int));
+    os.write((char*)&p_cols_, sizeof(int));
+    os.write((char*)&s_rows_, sizeof(int));
+    os.write((char*)&s_cols_, sizeof(int));
+    os.write((char*)&h_type,  sizeof(hlayer_type_e));
+    os.write((char*)&p_type_, sizeof(pad_type_e));
+    os.write((char*)&a_type,  sizeof(activation_type_e));
+    return;
+}
+
+void gctl::avgpooling::load_layer_setup(std::ifstream &is)
+{
+    int in_rows, in_cols, pool_rows, pool_cols, stride_rows, stride_cols;
+    hlayer_type_e h_type;
+    pad_type_e p_type;
+    activation_type_e a_type;
+
+    is.read((char*)&in_rows, sizeof(int));
+    is.read((char*)&in_cols, sizeof(int));
+    is.read((char*)&pool_rows, sizeof(int));
+    is.read((char*)&pool_cols, sizeof(int));
+    is.read((char*)&stride_rows, sizeof(int));
+    is.read((char*)&stride_cols, sizeof(int));
+    is.read((char*)&h_type, sizeof(hlayer_type_e));
+    is.read((char*)&p_type, sizeof(pad_type_e));
+    is.read((char*)&a_type, sizeof(activation_type_e));
+
+    if (h_type != AvgPooling)
+    {
+        throw std::invalid_argument("[gctl::avgpooling] Invalid hind layer type.");
+    }
+
+    init_avgpooling(in_rows, in_cols, pool_rows, pool_cols, stride_rows, stride_cols, p_type, a_type);
+    return;
+}
+
+void gctl::avgpooling::save_weights2text(const array<double> &all_weights, std::ofstream &os) const
+{
+    return;
+}
\ No newline at end of file
diff --git a/lib/dnn/hlayer_avgpooling.h b/lib/dnn/hlayer_avgpooling.h
new file mode 100644
index 0000000..bc6d6e1
--- /dev/null
+++ b/lib/dnn/hlayer_avgpooling.h
@@ -0,0 +1,68 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_HLAYER_AVGPOOLING_H
+#define _GCTL_DNN_HLAYER_AVGPOOLING_H
+
+#include "hlayer.h"
+
+namespace gctl
+{
+    class avgpooling : public dnn_hlayer
+    {
+    public:
+        avgpooling();
+        avgpooling(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, 
+            int stride_cols, pad_type_e pl_type, activation_type_e acti_type);
+        virtual ~avgpooling();
+
+        void forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data);
+        void backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+            const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data);
+
+        hlayer_type_e get_layer_type() const;
+        std::string get_layer_name() const;
+        std::string layer_info() const;
+        void save_layer_setup(std::ofstream &os) const;
+        void load_layer_setup(std::ifstream &is);
+        void save_weights2text(const array<double> &all_weights, std::ofstream &os) const;
+
+    private:
+        void init_avgpooling(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, 
+            int stride_cols, pad_type_e pl_type, activation_type_e acti_type);
+    
+    private:
+        int i_rows_, i_cols_, p_rows_, p_cols_, s_rows_, s_cols_;
+        int o_rows_, o_cols_;
+        int l_pad_, u_pad_;
+
+        array<int> p_idx_;
+        pad_type_e p_type_;
+    };
+}
+
+#endif // _GCTL_DNN_HLAYER_AVGPOOLING_H
\ No newline at end of file
diff --git a/lib/dnn/hlayer_convolution.cpp b/lib/dnn/hlayer_convolution.cpp
new file mode 100644
index 0000000..870b8d9
--- /dev/null
+++ b/lib/dnn/hlayer_convolution.cpp
@@ -0,0 +1,303 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "hlayer_convolution.h"
+
+gctl::convolution::convolution(){}
+
+gctl::convolution::convolution(int p_st, int channels, int in_rows, int in_cols, int filter_rows, int filter_cols, int stride_rows, int stride_cols, 
+    pad_type_e pl_type, activation_type_e acti_type)
+{
+    init_convolution(p_st, channels, in_rows, in_cols, filter_rows, filter_cols, stride_rows, stride_cols, pl_type, acti_type);
+}
+
+gctl::convolution::~convolution(){}
+
+void gctl::convolution::init_convolution(int p_st, int channels, int in_rows, int in_cols, int filter_rows, int filter_cols, 
+    int stride_rows, int stride_cols, pad_type_e pl_type, activation_type_e acti_type)
+{
+    chls_ = channels;
+    i_rows_ = in_rows; i_cols_ = in_cols;
+    f_rows_ = filter_rows; f_cols_ = filter_cols;
+    s_rows_ = stride_rows; s_cols_ = stride_cols;
+    p_type_ = pl_type;
+
+    if (acti_type == Identity) activator_ = new identity;
+    else if (acti_type == Mish) activator_ = new mish;
+    else if (acti_type == ReLU) activator_ = new relu;
+    else if (acti_type == PReLU) activator_ = new prelu;
+    else if (acti_type == Sigmoid) activator_ = new sigmoid;
+    else if (acti_type == SoftMax) activator_ = new softmax;
+    else if (acti_type == Tanh) activator_ = new tanh;
+    else throw std::invalid_argument("[gctl::convolution] Invalid activation type.");
+
+    // 计算输出大小
+    o_rows_ = (in_rows + s_rows_ - f_rows_)/s_rows_;
+    o_cols_ = (in_cols + s_cols_ - f_cols_)/s_cols_;
+    if (pl_type == Same && (in_rows + s_rows_ - f_rows_)%s_rows_ != 0)
+    {
+        o_rows_++;
+        u_pad_ = ((in_rows + s_rows_ - f_rows_)%s_rows_)/2;
+    }
+    else u_pad_ = 0;
+
+    if (pl_type == Same && (in_cols + s_cols_ - f_cols_)%s_cols_ != 0)
+    {
+        o_cols_++;
+        l_pad_ = ((in_cols + s_cols_ - f_cols_)%s_cols_)/2;
+    }
+    else l_pad_ = 0;
+
+    w_is_ = i_rows_*i_cols_;
+    w_outs_ = o_rows_*o_cols_;
+    w_st_ = p_st;
+    b_st_ = p_st + chls_*f_rows_*f_cols_; // Size of filter's parameters is chls_*f_rows_*f_cols_
+    p_idx_.resize(f_rows_*f_cols_);
+    return;
+}
+
+void gctl::convolution::forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data)
+{
+    if (prev_layer_data.col_size()%chls_ != 0)
+    {
+        throw std::invalid_argument("[gctl::convolution] Incompatible observation size to the channel size.");
+    }
+    
+    o_is_ = prev_layer_data.col_size()/chls_;
+    // Forward linear terms
+    // z_: out_size x nobs
+    // a_: out_size x nobs
+    z_.resize(w_outs_, o_is_);
+    a_.resize(w_outs_, o_is_);
+
+    double p_val;
+    for (size_t i = 0; i < o_rows_; i++)
+    {
+        for (size_t j = 0; j < o_cols_; j++)
+        {
+            if (p_type_ == Valid) cal_valid_padding_idx(p_idx_, i, j, i_rows_, i_cols_, f_rows_, f_cols_, s_rows_, s_cols_);
+            if (p_type_ == Same) cal_same_padding_idx(p_idx_, i, j, i_rows_, i_cols_, f_rows_, f_cols_, s_rows_, s_cols_, u_pad_, l_pad_);
+
+            for (size_t o = 0; o < o_is_; o++)
+            {
+                p_val = 0.0;
+                for (size_t c = 0; c < chls_; c++)
+                {
+                    for (size_t p = 0; p < f_rows_*f_cols_; p++)
+                    {
+                        if (p_idx_[p] != -1) p_val += all_weights[w_st_ + p + c*f_rows_*f_cols_] * prev_layer_data[p_idx_[p]][chls_*o + c];
+                    }
+                }
+
+                z_[i*o_cols_ + j][o] = p_val + all_weights[b_st_]; // one filter has one bias value
+            }
+        }
+    }
+
+    // Apply activation function
+    activator_->activate(z_, a_);
+    return;
+}
+
+void gctl::convolution::backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+    const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data)
+{
+    der_z_.resize(w_outs_, o_is_);
+    der_in_.resize(w_is_, o_is_*chls_);
+    der_in_.assign_all(0.0);
+
+    all_ders_count_.resize(chls_*f_rows_*f_cols_);
+    all_ders_count_.assign_all(0);
+
+    der_in_count_.resize(w_is_, o_is_*chls_);
+    der_in_count_.assign_all(0);
+
+    
+    // After forward stage, m_z contains z = conv(in, w) + b
+    // Now we need to calculate d(L) / d(z) = [d(a) / d(z)] * [d(L) / d(a)]
+    // d(L) / d(a) is computed in the next layer, contained in next_layer_data
+    // The Jacobian matrix J = d(a) / d(z) is determined by the activation function
+    activator_->apply_jacobian(z_, a_, next_layer_data, der_z_);
+
+    // z_j = sum_i(conv(in_i, w_ij)) + b_j
+    //
+    // d(z_k) / d(w_ij) = 0, if k != j
+    // d(L) / d(w_ij) = [d(z_j) / d(w_ij)] * [d(L) / d(z_j)] = sum_i{ [d(z_j) / d(w_ij)] * [d(L) / d(z_j)] }
+    // = sum_i(conv(in_i, d(L) / d(z_j)))
+    //
+    // z_j is an image (matrix), b_j is a scalar
+    // d(z_j) / d(b_j) = a matrix of the same size of d(z_j) filled with 1
+    // d(L) / d(b_j) = (d(L) / d(z_j)).sum()
+    //
+    // d(z_j) / d(in_i) = conv_full_op(w_ij_rotate)
+    // d(L) / d(in_i) = sum_j((d(z_j) / d(in_i)) * (d(L) / d(z_j))) = sum_j(conv_full(d(L) / d(z_j), w_ij_rotate))
+    // Derivative for weights
+
+    for (size_t h = 0; h < chls_; h++)
+    {
+        for (size_t p = 0; p < f_rows_*f_cols_; p++)
+        {
+            all_ders[w_st_ + h*f_rows_*f_cols_ + p] = 0.0;
+        }
+    }
+
+    for (size_t j = 0; j < o_is_; j++)
+    {
+        for (size_t h = 0; h < chls_; h++)
+        {
+            for (size_t r = 0; r < o_rows_; r++)
+            {
+                for (size_t c = 0; c < o_cols_; c++)
+                {
+                    if (p_type_ == Valid) cal_valid_padding_idx(p_idx_, r, c, i_rows_, i_cols_, f_rows_, f_cols_, s_rows_, s_cols_);
+                    if (p_type_ == Same) cal_same_padding_idx(p_idx_, r, c, i_rows_, i_cols_, f_rows_, f_cols_, s_rows_, s_cols_, u_pad_, l_pad_);
+
+                    for (size_t p = 0; p < f_rows_*f_cols_; p++)
+                    {
+                        if (p_idx_[p] != -1)
+                        {   
+                            all_ders[w_st_ + h*f_rows_*f_cols_ + p] += prev_layer_data[p_idx_[p]][j*chls_ + h]*der_z_[r*o_cols_ + c][j];
+                            all_ders_count_[h*f_rows_*f_cols_ + p]++;
+
+                            der_in_[p_idx_[p]][j*chls_ + h] += all_weights[w_st_ + h*f_rows_*f_cols_ + p]*der_z_[r*o_cols_ + c][j];
+                            der_in_count_[p_idx_[p]][j*chls_ + h]++;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    for (size_t h = 0; h < chls_; h++)
+    {
+        for (size_t p = 0; p < f_rows_*f_cols_; p++)
+        {
+            all_ders[w_st_ + h*f_rows_*f_cols_ + p] /= all_ders_count_[h*f_rows_*f_cols_ + p];
+        }
+    }
+
+    for (size_t j = 0; j < o_is_; j++)
+    {
+        for (size_t h = 0; h < chls_; h++)
+        {
+            for (size_t i = 0; i < w_outs_; i++)
+            {
+                if (der_in_count_[i][j*chls_ + h] != 0) der_in_[i][j*chls_ + h] /= der_in_count_[i][j*chls_ + h];
+            }
+        }
+    }
+
+    // Derivative for bias, d(L) / d(b) = d(L) / d(z)
+    all_ders[b_st_] = 0.0;
+    for (size_t i = 0; i < w_outs_; i++)
+    {
+        for (size_t j = 0; j < o_is_; j++)
+        {
+            all_ders[b_st_] += der_z_[i][j];
+        }
+    }
+    all_ders[b_st_] /= (o_is_*w_outs_);
+    return;
+}
+
+gctl::hlayer_type_e gctl::convolution::get_layer_type() const
+{
+    return Convolution;
+}
+
+std::string gctl::convolution::get_layer_name() const
+{
+    return "Convolution";
+}
+
+std::string gctl::convolution::layer_info() const
+{
+    std::string info = std::to_string(w_is_) + "x" + std::to_string(w_outs_)
+        + " (" + std::to_string(i_rows_) + "," + std::to_string(i_cols_) + ")x"
+        + "(" + std::to_string(o_rows_) + "," + std::to_string(o_cols_) + "), "
+        + std::to_string(f_rows_) + "x" + std::to_string(f_cols_)
+        + ", " + std::to_string(s_rows_) + "x" + std::to_string(s_cols_)
+        + ", " + std::to_string(chls_) + " channel(s), Convolution ";
+
+    if (p_type_ == Same) info += "(Same), ";
+    if (p_type_ == Valid) info += "(Valid), ";
+
+    info += activator_->activation_name();
+    return info;
+}
+
+void gctl::convolution::save_layer_setup(std::ofstream &os) const
+{
+    hlayer_type_e h_type = get_layer_type();
+    activation_type_e a_type = get_activation_type();
+
+    os.write((char*)&w_st_, sizeof(int));
+    os.write((char*)&chls_, sizeof(int));
+    os.write((char*)&i_rows_, sizeof(int));
+    os.write((char*)&i_cols_, sizeof(int));
+    os.write((char*)&f_rows_, sizeof(int));
+    os.write((char*)&f_cols_, sizeof(int));
+    os.write((char*)&s_rows_, sizeof(int));
+    os.write((char*)&s_cols_, sizeof(int));
+    os.write((char*)&h_type,  sizeof(hlayer_type_e));
+    os.write((char*)&p_type_, sizeof(pad_type_e));
+    os.write((char*)&a_type,  sizeof(activation_type_e));
+    return;
+}
+
+void gctl::convolution::load_layer_setup(std::ifstream &is)
+{
+    int st, chls, in_rows, in_cols, filter_rows, filter_cols, stride_rows, stride_cols;
+    hlayer_type_e h_type;
+    pad_type_e p_type;
+    activation_type_e a_type;
+
+    is.read((char*)&st, sizeof(int));
+    is.read((char*)&chls, sizeof(int));
+    is.read((char*)&in_rows, sizeof(int));
+    is.read((char*)&in_cols, sizeof(int));
+    is.read((char*)&filter_rows, sizeof(int));
+    is.read((char*)&filter_cols, sizeof(int));
+    is.read((char*)&stride_rows, sizeof(int));
+    is.read((char*)&stride_cols, sizeof(int));
+    is.read((char*)&h_type, sizeof(hlayer_type_e));
+    is.read((char*)&p_type, sizeof(pad_type_e));
+    is.read((char*)&a_type, sizeof(activation_type_e));
+
+    if (h_type != Convolution)
+    {
+        throw std::invalid_argument("[gctl::convolution] Invalid hind layer type.");
+    }
+
+    init_convolution(st, chls, in_rows, in_cols, filter_rows, filter_cols, stride_rows, stride_cols, p_type, a_type);
+    return;
+}
+
+void gctl::convolution::save_weights2text(const array<double> &all_weights, std::ofstream &os) const
+{
+    return;
+}
\ No newline at end of file
diff --git a/lib/dnn/hlayer_convolution.h b/lib/dnn/hlayer_convolution.h
new file mode 100644
index 0000000..7469e3b
--- /dev/null
+++ b/lib/dnn/hlayer_convolution.h
@@ -0,0 +1,74 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_HLAYER_CONVOLUTION_H
+#define _GCTL_DNN_HLAYER_CONVOLUTION_H
+
+#include "hlayer.h"
+
+namespace gctl
+{
+    class convolution : public dnn_hlayer
+    {
+    public:
+        convolution();
+        convolution(int p_st, int channels, int in_rows, int in_cols, int filter_rows, int filter_cols, 
+            int stride_rows, int stride_cols, pad_type_e pl_type, activation_type_e acti_type);
+        virtual ~convolution();
+
+        void forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data);
+        void backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+            const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data);
+
+        hlayer_type_e get_layer_type() const;
+        std::string get_layer_name() const;
+        std::string layer_info() const;
+        void save_layer_setup(std::ofstream &os) const;
+        void load_layer_setup(std::ifstream &is);
+        void save_weights2text(const array<double> &all_weights, std::ofstream &os) const;
+
+    private:
+        void init_convolution(int p_st, int channels, int in_rows, int in_cols, int filter_rows, int filter_cols, int stride_rows, 
+            int stride_cols, pad_type_e pl_type, activation_type_e acti_type);
+    
+    private:
+        int i_rows_, i_cols_, f_rows_, f_cols_, s_rows_, s_cols_;
+        int o_rows_, o_cols_;
+        int l_pad_, u_pad_;
+        int chls_;
+        int w_st_; // weight start index
+        int b_st_; // bias start index
+
+        array<int>  all_ders_count_;
+        matrix<int> der_in_count_;
+
+        array<int> p_idx_;
+        pad_type_e p_type_;
+    };
+}
+
+#endif // _GCTL_DNN_HLAYER_CONVOLUTION_H
\ No newline at end of file
diff --git a/lib/dnn/hlayer_fully_connected.cpp b/lib/dnn/hlayer_fully_connected.cpp
new file mode 100644
index 0000000..44abe62
--- /dev/null
+++ b/lib/dnn/hlayer_fully_connected.cpp
@@ -0,0 +1,216 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "hlayer_fully_connected.h"
+
+gctl::fully_connected::fully_connected(){}
+
+gctl::fully_connected::fully_connected(int p_st, int p_ins, int p_outs, activation_type_e acti_type)
+{
+    init_fully_connected(p_st, p_ins, p_outs, acti_type);
+}
+
+gctl::fully_connected::~fully_connected(){}
+
+void gctl::fully_connected::init_fully_connected(int p_st, int p_ins, int p_outs, activation_type_e acti_type)
+{
+    w_st_ = p_st; w_is_ = p_ins; w_outs_ = p_outs;
+    b_st_ = p_st + p_ins*p_outs;
+
+    if (acti_type == Identity) activator_ = new identity;
+    else if (acti_type == Mish) activator_ = new mish;
+    else if (acti_type == ReLU) activator_ = new relu;
+    else if (acti_type == PReLU) activator_ = new prelu;
+    else if (acti_type == Sigmoid) activator_ = new sigmoid;
+    else if (acti_type == SoftMax) activator_ = new softmax;
+    else if (acti_type == Tanh) activator_ = new tanh;
+    else throw std::invalid_argument("[gctl::fully_connected] Invalid activation type.");
+    return;
+}
+
+void gctl::fully_connected::forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data)
+{
+    // z_: out_size x nobs
+    // a_: out_size x nobs
+    o_is_ = prev_layer_data.col_size();
+    // Forward linear terms
+    z_.resize(w_outs_, o_is_);
+    a_.resize(w_outs_, o_is_);
+
+    // Linear term z = W^T * in + b
+    int i, j, k;
+#pragma omp parallel for private (i, j, k) schedule(guided)
+    for (i = 0; i < w_outs_; i++)
+    {
+        for (j = 0; j < o_is_; j++)
+        {
+            z_[i][j] = 0.0;
+            for (k = 0; k < w_is_; k++)
+            {
+                z_[i][j] += all_weights[w_st_ + i + k*w_outs_]*prev_layer_data[k][j];
+            }
+        }
+    }
+
+//#pragma omp parallel for private (i, j) schedule(guided)
+    for (j = 0; j < o_is_; j++)
+    {
+        for (i = 0; i < w_outs_; i++)
+        {
+            z_[i][j] += all_weights[b_st_ + i];
+        }
+    }
+
+    // Apply activation function
+    activator_->activate(z_, a_);
+/*
+    for (j = 0; j < o_is_; j++)
+    {
+        for (i = 0; i < w_outs_; i++)
+        {
+            std::cout << a_[i][j] << "\n";
+        }
+    }
+
+    std::cout << "done\n";
+*/
+    return;
+}
+
+void gctl::fully_connected::backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+    const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data)
+{
+    der_z_.resize(w_outs_, o_is_);
+    der_in_.resize(w_is_, o_is_);
+    // prev_layer_data: in_size x nobs
+    // next_layer_data: out_size x nobs
+    // After forward stage, m_z contains z = W' * in + b
+    // Now we need to calculate d(L) / d(z) = [d(a) / d(z)] * [d(L) / d(a)]
+    // d(L) / d(a) is computed in the next layer, contained in next_layer_data
+    // The Jacobian matrix J = d(a) / d(z) is determined by the activation function
+    // der_z_: out_size x nobs   
+    activator_->apply_jacobian(z_, a_, next_layer_data, der_z_);
+    // Now dLz contains d(L) / d(z)
+    // Derivative for weights, d(L) / d(W) = [d(L) / d(z)] * in'
+    int i, j, k;
+#pragma omp parallel for private (i, j, k) schedule(guided)
+    for (i = 0; i < w_is_; i++)
+    {
+        for (j = 0; j < w_outs_; j++)
+        {
+            all_ders[w_st_ + i*w_outs_ + j] = 0;
+            for (k = 0; k < o_is_; k++)
+            {
+                all_ders[w_st_ + i*w_outs_ + j] += prev_layer_data[i][k]*der_z_[j][k];
+            }
+            all_ders[w_st_ + i*w_outs_ + j] /= o_is_;
+        }
+    }
+    // Derivative for bias, d(L) / d(b) = d(L) / d(z)
+#pragma omp parallel for private (i, j) schedule(guided)
+    for (i = 0; i < w_outs_; i++)
+    {
+        all_ders[b_st_ + i] = 0.0;
+        for (j = 0; j < o_is_; j++)
+        {
+            all_ders[b_st_ + i] += der_z_[i][j];
+        }
+        all_ders[b_st_ + i] /= o_is_;
+    }
+    // Compute d(L) / d_in = W * [d(L) / d(z)]
+    // der_in_: in_size x nobs
+#pragma omp parallel for private (i, j, k) schedule(guided)
+    for (i = 0; i < w_is_; i++)
+    {
+        for (j = 0; j < o_is_; j++)
+        {
+            der_in_[i][j] = 0.0;
+            for (k = 0; k < w_outs_; k++)
+            {
+                der_in_[i][j] += all_weights[w_st_ + i*w_outs_ + k]*der_z_[k][j];
+            }
+        }
+    }
+    return;
+}
+
+gctl::hlayer_type_e gctl::fully_connected::get_layer_type() const
+{
+    return FullyConnected;
+}
+
+std::string gctl::fully_connected::get_layer_name() const
+{
+    return "FullyConnected";
+}
+
+std::string gctl::fully_connected::layer_info() const
+{
+    std::string info = std::to_string(w_is_) + "x" + std::to_string(w_outs_) + ", FullyConnected, " + activator_->activation_name();
+    return info;
+}
+
+void gctl::fully_connected::save_layer_setup(std::ofstream &os) const
+{
+    hlayer_type_e l_type = get_layer_type();
+    activation_type_e a_type = get_activation_type();
+
+    os.write((char*)&w_st_, sizeof(int));
+    os.write((char*)&w_is_, sizeof(int));
+    os.write((char*)&w_outs_, sizeof(int));
+    os.write((char*)&l_type, sizeof(hlayer_type_e));
+    os.write((char*)&a_type, sizeof(activation_type_e));
+    return;
+}
+
+void gctl::fully_connected::load_layer_setup(std::ifstream &is)
+{
+    int st, iss, outs;
+    hlayer_type_e l_type;
+    activation_type_e a_type;
+
+    is.read((char*)&st, sizeof(int));
+    is.read((char*)&iss, sizeof(int));
+    is.read((char*)&outs, sizeof(int));
+    is.read((char*)&l_type, sizeof(hlayer_type_e));
+    is.read((char*)&a_type, sizeof(activation_type_e));
+
+    init_fully_connected(st, iss, outs, a_type);
+    return;
+}
+
+void gctl::fully_connected::save_weights2text(const array<double> &all_weights, std::ofstream &os) const
+{
+    for (int i = 0; i < w_outs_; i++)
+    {
+        for (int k = 0; k < w_is_; k++)
+        {
+            os << k << " " << i << " " << all_weights[w_st_ + i + k*w_outs_] << "\n";
+        }
+    }
+    return;
+}
\ No newline at end of file
diff --git a/lib/dnn/hlayer_fully_connected.h b/lib/dnn/hlayer_fully_connected.h
new file mode 100644
index 0000000..c7c52bf
--- /dev/null
+++ b/lib/dnn/hlayer_fully_connected.h
@@ -0,0 +1,62 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_HLAYER_FULLY_CONNECTED_H
+#define _GCTL_DNN_HLAYER_FULLY_CONNECTED_H
+
+#include "hlayer.h"
+
+namespace gctl
+{
+    class fully_connected : public dnn_hlayer
+    {
+    public:
+        fully_connected();
+        fully_connected(int p_st, int p_ins, int p_outs, activation_type_e acti_type);
+        virtual ~fully_connected();
+
+        void forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data);
+        void backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+            const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data);
+
+        hlayer_type_e get_layer_type() const;
+        std::string get_layer_name() const;
+        std::string layer_info() const;
+        void save_layer_setup(std::ofstream &os) const;
+        void load_layer_setup(std::ifstream &is);
+        void save_weights2text(const array<double> &all_weights, std::ofstream &os) const;
+    
+    private:
+        void init_fully_connected(int p_st, int p_ins, int p_outs, activation_type_e acti_type);
+    
+    private:
+        int w_st_; // weight start index
+        int b_st_; // bias start index
+    };
+}
+
+#endif // _GCTL_DNN_HLAYER_FULLY_CONNECTED_H
\ No newline at end of file
diff --git a/lib/dnn/hlayer_maxpooling.cpp b/lib/dnn/hlayer_maxpooling.cpp
new file mode 100644
index 0000000..3c42bdd
--- /dev/null
+++ b/lib/dnn/hlayer_maxpooling.cpp
@@ -0,0 +1,220 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "hlayer_maxpooling.h"
+
+gctl::maxpooling::maxpooling(){}
+
+gctl::maxpooling::maxpooling(int in_rows, int in_cols, int pool_rows, int pool_cols, 
+    int stride_rows, int stride_cols, pad_type_e pl_type, activation_type_e acti_type)
+{
+    init_maxpooling(in_rows, in_cols, pool_rows, pool_cols, stride_rows, stride_cols, pl_type, acti_type);
+}
+
+gctl::maxpooling::~maxpooling(){}
+
+void gctl::maxpooling::init_maxpooling(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, 
+    int stride_cols, pad_type_e pl_type, activation_type_e acti_type)
+{
+    i_rows_ = in_rows; i_cols_ = in_cols;
+    p_rows_ = pool_rows; p_cols_ = pool_cols;
+    s_rows_ = stride_rows; s_cols_ = stride_cols;
+    p_type_ = pl_type;
+
+    if (acti_type == Identity) activator_ = new identity;
+    else if (acti_type == Mish) activator_ = new mish;
+    else if (acti_type == ReLU) activator_ = new relu;
+    else if (acti_type == PReLU) activator_ = new prelu;
+    else if (acti_type == Sigmoid) activator_ = new sigmoid;
+    else if (acti_type == SoftMax) activator_ = new softmax;
+    else if (acti_type == Tanh) activator_ = new tanh;
+    else throw std::invalid_argument("[gctl::maxpooling] Invalid activation type.");
+
+    // 计算输出大小
+    o_rows_ = (in_rows + s_rows_ - p_rows_)/s_rows_;
+    o_cols_ = (in_cols + s_cols_ - p_cols_)/s_cols_;
+    if (pl_type == Same && (in_rows + s_rows_ - p_rows_)%s_rows_ != 0)
+    {
+        o_rows_++;
+        u_pad_ = ((in_rows + s_rows_ - p_rows_)%s_rows_)/2;
+    }
+    else u_pad_ = 0;
+
+    if (pl_type == Same && (in_cols + s_cols_ - p_cols_)%s_cols_ != 0)
+    {
+        o_cols_++;
+        l_pad_ = ((in_cols + s_cols_ - p_cols_)%s_cols_)/2;
+    }
+    else l_pad_ = 0;
+
+    w_is_ = i_rows_*i_cols_;
+    w_outs_ = o_rows_*o_cols_;
+    p_idx_.resize(p_rows_*p_cols_);
+    return;
+}
+
+void gctl::maxpooling::forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data)
+{
+    // z_: out_size x nobs
+    // a_: out_size x nobs
+    o_is_ = prev_layer_data.col_size();
+    // Forward linear terms
+    z_.resize(w_outs_, o_is_);
+    a_.resize(w_outs_, o_is_);
+    p_loc_.resize(w_outs_, o_is_);
+
+    double p_val;
+    for (size_t i = 0; i < o_rows_; i++)
+    {
+        for (size_t j = 0; j < o_cols_; j++)
+        {
+            if (p_type_ == Valid) cal_valid_padding_idx(p_idx_, i, j, i_rows_, i_cols_, p_rows_, p_cols_, s_rows_, s_cols_);
+            if (p_type_ == Same) cal_same_padding_idx(p_idx_, i, j, i_rows_, i_cols_, p_rows_, p_cols_, s_rows_, s_cols_, u_pad_, l_pad_);
+
+            for (size_t o = 0; o < o_is_; o++)
+            {
+                p_val = GCTL_BDL_MIN;
+                for (size_t p = 0; p < p_idx_.size(); p++)
+                {
+                    if (p_idx_[p] != -1 && prev_layer_data[p_idx_[p]][o] > p_val)
+                    {
+                        p_val = prev_layer_data[p_idx_[p]][o];
+                        p_loc_[i*o_cols_ + j][o] = p_idx_[p];
+                    }
+                }
+
+                z_[i*o_cols_ + j][o] = p_val;
+            }
+        }
+    }
+
+    // Apply activation function
+    activator_->activate(z_, a_);
+    return;
+}
+
+void gctl::maxpooling::backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+    const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data)
+{
+    der_z_.resize(w_outs_, o_is_);
+    der_in_.resize(w_is_, o_is_, 0.0);
+    // prev_layer_data: in_size x nobs
+    // next_layer_data: out_size x nobs
+    // After forward stage, m_z contains z = max_maxpooling(in)
+    // Now we need to calculate d(L) / d(z) = [d(a) / d(z)] * [d(L) / d(a)]
+    // d(L) / d(z) is computed in the next layer, contained in next_layer_data
+    // The Jacobian matrix J = d(a) / d(z) is determined by the activation function
+    // der_z_: out_size x nobs   
+    activator_->apply_jacobian(z_, a_, next_layer_data, der_z_);
+
+    // d(L) / d(in_i) = sum_j{ [d(z_j) / d(in_i)] * [d(L) / d(z_j)] }
+    // d(z_j) / d(in_i) = 1 if in_i is used to compute z_j and is the maximum
+    //                  = 0 otherwise
+    for (size_t j = 0; j < o_is_; j++)
+    {
+        for (size_t i = 0; i < w_outs_; i++)
+        {
+            der_in_[p_loc_[i][j]][j] += der_z_[i][j];
+        }
+    }
+    return;
+}
+
+gctl::hlayer_type_e gctl::maxpooling::get_layer_type() const
+{
+    return MaxPooling;
+}
+
+std::string gctl::maxpooling::get_layer_name() const
+{
+    return "MaxPooling";
+}
+
+std::string gctl::maxpooling::layer_info() const
+{
+    std::string info = std::to_string(w_is_) + "x" + std::to_string(w_outs_)
+        + " (" + std::to_string(i_rows_) + "," + std::to_string(i_cols_) + ")x"
+        + "(" + std::to_string(o_rows_) + "," + std::to_string(o_cols_) + "), "
+        + std::to_string(p_rows_) + "x" + std::to_string(p_cols_)
+        + ", " + std::to_string(s_rows_) + "x" + std::to_string(s_cols_)
+        + ", MaxPooling ";
+
+    if (p_type_ == Same) info += "(Same), ";
+    if (p_type_ == Valid) info += "(Valid), ";
+
+    info += activator_->activation_name();
+    return info;
+}
+
+void gctl::maxpooling::save_layer_setup(std::ofstream &os) const
+{
+    hlayer_type_e h_type = get_layer_type();
+    activation_type_e a_type = get_activation_type();
+
+    os.write((char*)&i_rows_, sizeof(int));
+    os.write((char*)&i_cols_, sizeof(int));
+    os.write((char*)&p_rows_, sizeof(int));
+    os.write((char*)&p_cols_, sizeof(int));
+    os.write((char*)&s_rows_, sizeof(int));
+    os.write((char*)&s_cols_, sizeof(int));
+    os.write((char*)&h_type,  sizeof(hlayer_type_e));
+    os.write((char*)&p_type_, sizeof(pad_type_e));
+    os.write((char*)&a_type,  sizeof(activation_type_e));
+    return;
+}
+
+void gctl::maxpooling::load_layer_setup(std::ifstream &is)
+{
+    int in_rows, in_cols, pool_rows, pool_cols, stride_rows, stride_cols;
+    hlayer_type_e h_type;
+    pad_type_e p_type;
+    activation_type_e a_type;
+
+    is.read((char*)&in_rows, sizeof(int));
+    is.read((char*)&in_cols, sizeof(int));
+    is.read((char*)&pool_rows, sizeof(int));
+    is.read((char*)&pool_cols, sizeof(int));
+    is.read((char*)&stride_rows, sizeof(int));
+    is.read((char*)&stride_cols, sizeof(int));
+    is.read((char*)&h_type, sizeof(hlayer_type_e));
+    is.read((char*)&p_type, sizeof(pad_type_e));
+    is.read((char*)&a_type, sizeof(activation_type_e));
+
+    if (h_type != MaxPooling)
+    {
+        throw std::invalid_argument("[gctl::maxpooling] Invalid hind layer type.");
+    }
+
+    init_maxpooling(in_rows, in_cols, pool_rows, pool_cols, stride_rows, stride_cols, p_type, a_type);
+    return;
+}
+
+void gctl::maxpooling::save_weights2text(const array<double> &all_weights, std::ofstream &os) const
+{
+    
+    return;
+}
\ No newline at end of file
diff --git a/lib/dnn/hlayer_maxpooling.h b/lib/dnn/hlayer_maxpooling.h
new file mode 100644
index 0000000..d0656df
--- /dev/null
+++ b/lib/dnn/hlayer_maxpooling.h
@@ -0,0 +1,69 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_HLAYER_MAXPOOLING_H
+#define _GCTL_DNN_HLAYER_MAXPOOLING_H
+
+#include "hlayer.h"
+
+namespace gctl
+{
+    class maxpooling : public dnn_hlayer
+    {
+    public:
+        maxpooling();
+        maxpooling(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, 
+            int stride_cols, pad_type_e pl_type, activation_type_e acti_type);
+        virtual ~maxpooling();
+
+        void forward_propagation(const array<double> &all_weights, const matrix<double> &prev_layer_data);
+        void backward_propagation(const array<double> &all_weights, const array<double> &all_ders, 
+            const matrix<double> &prev_layer_data, const matrix<double> &next_layer_data);
+
+        hlayer_type_e get_layer_type() const;
+        std::string get_layer_name() const;
+        std::string layer_info() const;
+        void save_layer_setup(std::ofstream &os) const;
+        void load_layer_setup(std::ifstream &is);
+        void save_weights2text(const array<double> &all_weights, std::ofstream &os) const;
+
+    private:
+        void init_maxpooling(int in_rows, int in_cols, int pool_rows, int pool_cols, int stride_rows, 
+            int stride_cols, pad_type_e pl_type, activation_type_e acti_type);
+    
+    private:
+        int i_rows_, i_cols_, p_rows_, p_cols_, s_rows_, s_cols_;
+        int o_rows_, o_cols_;
+        int l_pad_, u_pad_;
+
+        array<int> p_idx_;
+        matrix<int> p_loc_;
+        pad_type_e p_type_;
+    };
+}
+
+#endif // _GCTL_DNN_HLAYER_MAXPOOLING_H
\ No newline at end of file
diff --git a/lib/dnn/olayer.cpp b/lib/dnn/olayer.cpp
new file mode 100644
index 0000000..8afc9c1
--- /dev/null
+++ b/lib/dnn/olayer.cpp
@@ -0,0 +1,49 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "olayer.h"
+
+gctl::dnn_olayer::dnn_olayer() {}
+
+gctl::dnn_olayer::~dnn_olayer() {}
+
+const gctl::matrix<double> &gctl::dnn_olayer::backward_propagation_data()
+{
+    return der_in_;
+}
+
+void gctl::dnn_olayer::check_target_data(const matrix<double> &target) {}
+
+void gctl::dnn_olayer::check_target_data(const array<int> &target)
+{
+    throw std::invalid_argument("[gctl::dnn_olayer] This output type cannot take class labels as target data.");
+}
+
+void gctl::dnn_olayer::evaluation(const matrix<double> &prev_layer_data, const array<int> &target)
+{
+    throw std::invalid_argument("[gctl::dnn_olayer] This output type cannot take class labels as target data.");
+}
\ No newline at end of file
diff --git a/lib/dnn/olayer.h b/lib/dnn/olayer.h
new file mode 100644
index 0000000..8b1bcdc
--- /dev/null
+++ b/lib/dnn/olayer.h
@@ -0,0 +1,89 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_OLAYER_H
+#define _GCTL_DNN_OLAYER_H
+
+#include "gctl/core.h"
+#include "gctl/algorithm.h"
+
+namespace gctl
+{
+    enum olayer_type_e
+    {
+        RegressionMSE,
+        MultiClassEntropy,
+        BinaryClassEntropy,
+    };
+
+    class dnn_olayer
+    {
+    public:
+        dnn_olayer();
+        virtual ~dnn_olayer();
+
+        // The derivative of the input of this layer, which is also the derivative
+        // of the output of previous layer
+        const matrix<double> &backward_propagation_data();
+
+        // Check the format of target data, e.g. in classification problems the
+        // target data should be binary (either 0 or 1)
+        virtual void check_target_data(const matrix<double> &target);
+
+        // Another type of target data where each element is a class label
+        // This version may not be sensible for regression tasks, so by default
+        // we raise an exception
+        virtual void check_target_data(const array<int> &target);
+
+        // A combination of the forward stage and the back-propagation stage for the output layer
+        // The computed derivative of the input should be stored in this layer, and can be retrieved by
+        // the backward_propagation_data() function
+        virtual void evaluation(const matrix<double> &prev_layer_data, const matrix<double> &target) = 0;
+
+        // Another type of target data where each element is a class label
+        // This version may not be sensible for regression tasks, so by default
+        // we raise an exception
+        virtual void evaluation(const matrix<double> &prev_layer_data, const array<int> &target);
+
+        // Return the loss function value after the evaluation
+        // This function can be assumed to be called after evaluate(), so that it can make use of the
+        // intermediate result to save some computation
+        virtual double loss_value() const = 0;
+
+        // Return the output layer name. It is used to export the NN model.
+        virtual std::string get_output_name() const = 0;
+        
+        // Return the output layer type. It is used to export the NN model.
+        virtual olayer_type_e get_output_type() const = 0;
+    
+    protected:
+        matrix<double> der_in_; // Derivative of the input of this layer
+        // Note that input of this layer is also the output of previous layer
+    };
+}
+
+#endif // _GCTL_DNN_OLAYER_H
\ No newline at end of file
diff --git a/lib/dnn/olayer_binaryentropy.cpp b/lib/dnn/olayer_binaryentropy.cpp
new file mode 100644
index 0000000..ca6b20e
--- /dev/null
+++ b/lib/dnn/olayer_binaryentropy.cpp
@@ -0,0 +1,160 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "olayer_binaryentropy.h"
+
+gctl::binary_entropy::binary_entropy() {}
+
+gctl::binary_entropy::~binary_entropy() {}
+
+void gctl::binary_entropy::check_target_data(const matrix<double> &target)
+{
+    int nobs = target.col_size();
+    int ncls = target.row_size();
+
+    // Each element should be either 0 or 1
+    for (size_t i = 0; i < nobs; i++)
+    {
+        for (size_t j = 0; j < ncls; j++)
+        {
+            if (fabs(target[j][i] - 1.0) > 1e-6 && fabs(target[j][i]) > 1e-6)
+            {
+                throw std::invalid_argument("[gctl::binary_entropy] Target data should only contain zero or one.");
+            }
+        }
+    }
+    return;
+}
+
+void gctl::binary_entropy::check_target_data(const array<int> &target)
+{
+    int nobs = target.size();
+    for (size_t i = 0; i < nobs; i++)
+    {
+        if (target[i] != 0 && target[i] != 1)
+        {
+            throw std::invalid_argument("[gctl::binary_entropy] Target data should only contain zero or one.");
+        }
+    }
+    return;
+}
+
+void gctl::binary_entropy::evaluation(const matrix<double> &prev_layer_data, const matrix<double> &target)
+{
+    int nobs = prev_layer_data.col_size();
+    int ncls = prev_layer_data.row_size();
+
+    if (target.col_size() != nobs || target.row_size() != ncls)
+    {
+        throw std::invalid_argument("[gctl::binary_entropy] Target data have incorrect dimension.");
+    }
+    
+    // Compute the derivative of the input of this layer
+    // L = -y * log(phat) - (1 - y) * log(1 - phat)
+    // in = phat
+    // d（L） / d（in） = -y / phat + (1 - y) / (1 - phat), y is either 0 or 1
+    der_in_.resize(ncls, nobs);
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided)
+    for (i = 0; i < ncls; i++)
+    {
+        for (j = 0; j < nobs; j++)
+        {
+            der_in_[i][j] = -1.0*target[i][j]/prev_layer_data[i][j] + (1.0 - target[i][j])/(1.0 - prev_layer_data[i][j]);
+        }
+    }
+    return;
+}
+
+void gctl::binary_entropy::evaluation(const matrix<double> &prev_layer_data, const array<int> &target)
+{
+    // target is a vector of class labels that take values from [0, 1, ..., nclass - 1]
+    // The i-th element of target is the class label for observation i
+    int nobs = prev_layer_data.col_size();
+    int ncls = prev_layer_data.row_size();
+
+    if (ncls != 1)
+    {
+        throw std::invalid_argument("[gctl::binary_entropy] Only one response variable is allowed when class labels are used as target data.");
+    }
+
+    if (target.size() != nobs)
+    {
+        throw std::invalid_argument("[gctl::binary_entropy] Target data have incorrect dimension.");
+    }
+
+    // Compute the derivative of the input of this layer
+    // L = -log(phat[y])
+    // in = phat
+    // d(L) / d(in) = [0, 0, ..., -1/phat[y], 0, ..., 0]
+    der_in_.resize(ncls, nobs);
+    der_in_.assign_all(0.0);
+
+    int j;
+#pragma omp parallel for private (j) schedule(guided)
+    for (j = 0; j < nobs; j++)
+    {
+        der_in_[0][j] = -1.0*target[j]/prev_layer_data[0][j] + (1.0 - target[j])/(1.0 - prev_layer_data[0][j]);
+    }
+    return;
+}
+
+double gctl::binary_entropy::loss_value() const
+{
+    // L = -y * log(phat) - (1 - y) * log(1 - phat)
+    // y = 0 => L = -log(1 - phat)
+    // y = 1 => L = -log(phat)
+    // m_din contains 1/(1 - phat) if y = 0, and -1/phat if y = 1, so
+    // L = log(abs(m_din)).sum()
+    double res = 0.0;
+
+    int nobs = der_in_.col_size();
+    int ncls = der_in_.row_size();
+
+    int i, j;
+//#pragma omp parallel for private (i, j) schedule(guided)
+    for (i = 0; i < ncls; i++)
+    {
+        for (j = 0; j < nobs; j++)
+        {
+            res += std::log(fabs(der_in_[i][j]));
+        }
+    }
+
+    return res/der_in_.col_size();
+}
+
+std::string gctl::binary_entropy::get_output_name() const
+{
+    return "BinaryClassEntropy";
+}
+
+gctl::olayer_type_e gctl::binary_entropy::get_output_type() const
+{
+    return BinaryClassEntropy;
+}
\ No newline at end of file
diff --git a/lib/dnn/olayer_binaryentropy.h b/lib/dnn/olayer_binaryentropy.h
new file mode 100644
index 0000000..78d69d2
--- /dev/null
+++ b/lib/dnn/olayer_binaryentropy.h
@@ -0,0 +1,73 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_OLAYER_BINARYENTROPY_H
+#define _GCTL_DNN_OLAYER_BINARYENTROPY_H
+
+#include "olayer.h"
+
+namespace gctl
+{
+    class binary_entropy : public dnn_olayer
+    {
+    public:
+        binary_entropy();
+        virtual ~binary_entropy();
+
+        // Check the format of target data, e.g. in classification problems the
+        // target data should be binary (either 0 or 1)
+        void check_target_data(const matrix<double> &target);
+
+        // Another type of target data where each element is a class label
+        // This version may not be sensible for regression tasks, so by default
+        // we raise an exception
+        void check_target_data(const array<int> &target);
+
+        // A combination of the forward stage and the back-propagation stage for the output layer
+        // The computed derivative of the input should be stored in this layer, and can be retrieved by
+        // the backward_propagation_data() function
+        void evaluation(const matrix<double> &prev_layer_data, const matrix<double> &target);
+
+        // Another type of target data where each element is a class label
+        // This version may not be sensible for regression tasks, so by default
+        // we raise an exception
+        void evaluation(const matrix<double> &prev_layer_data, const array<int> &target);
+
+        // Return the loss function value after the evaluation
+        // This function can be assumed to be called after evaluate(), so that it can make use of the
+        // intermediate result to save some computation
+        double loss_value() const;
+
+        // Return the output layer name. It is used to export the NN model.
+        std::string get_output_name() const;
+
+        // Return the output layer type. It is used to export the NN model.
+        olayer_type_e get_output_type() const;
+    };
+}
+
+#endif // _GCTL_DNN_OLAYER_BINARYENTROPY_H
\ No newline at end of file
diff --git a/lib/dnn/olayer_multientropy.cpp b/lib/dnn/olayer_multientropy.cpp
new file mode 100644
index 0000000..e5e5d0e
--- /dev/null
+++ b/lib/dnn/olayer_multientropy.cpp
@@ -0,0 +1,163 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "olayer_multientropy.h"
+
+gctl::multi_entropy::multi_entropy() {}
+
+gctl::multi_entropy::~multi_entropy() {}
+
+void gctl::multi_entropy::check_target_data(const matrix<double> &target)
+{
+    int nobs = target.col_size();
+    int ncls = target.row_size();
+
+    int one;
+    for (size_t i = 0; i < nobs; i++)
+    {
+        one = 0;
+        for (size_t j = 0; j < ncls; j++)
+        {
+            if (fabs(target[j][i] - 1.0) < 1e-6)
+            {
+                one++;
+                continue;
+            }
+            
+            if (fabs(target[j][i]) > 1e-6)
+            {
+                throw std::invalid_argument("[gctl::multi_entropy] Target data should only contain zero or one.");
+            }
+        }
+        
+        if (one != 1)
+        {
+            throw std::invalid_argument("[gctl::multi_entropy] Each column of target data should only contain one \"1\".");
+        }
+    }
+    return;
+}
+
+void gctl::multi_entropy::check_target_data(const array<int> &target)
+{
+    int nobs = target.size();
+    for (size_t i = 0; i < nobs; i++)
+    {
+        if (target[i] < 0)
+        {
+            throw std::invalid_argument("[gctl::multi_entropy] Target data must be non-negative.");
+        }
+    }
+    return;
+}
+
+void gctl::multi_entropy::evaluation(const matrix<double> &prev_layer_data, const matrix<double> &target)
+{
+    int nobs = prev_layer_data.col_size();
+    int ncls = prev_layer_data.row_size();
+
+    if (target.col_size() != nobs || target.row_size() != ncls)
+    {
+        throw std::invalid_argument("[gctl::multi_entropy] Target data have incorrect dimension.");
+    }
+
+    // Compute the derivative of the input of this layer
+    // L = -sum(log(phat) * y)
+    // in = phat
+    // d(L) / d(in) = -y / phat
+    der_in_.resize(ncls, nobs, 0.0);
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided)
+    for (i = 0; i < ncls; i++)
+    {
+        for (j = 0; j < nobs; j++)
+        {
+            der_in_[i][j] = -1.0*target[i][j]/prev_layer_data[i][j];
+        }
+    }
+
+    // L = -sum(log(phat) * y)
+    // in = phat
+    // d(L) / d(in) = -y / phat
+    // m_din contains 0 if y = 0, and -1/phat if y = 1
+    loss_ = 0.0;
+//#pragma omp parallel for private (i, j) schedule(guided)
+    for (i = 0; i < ncls; i++)
+    {
+        for (j = 0; j < nobs; j++)
+        {
+            loss_ -= std::log(prev_layer_data[i][j])*target[i][j];
+        }
+    }
+
+    loss_ /= der_in_.col_size();
+    return;
+}
+
+void gctl::multi_entropy::evaluation(const matrix<double> &prev_layer_data, const array<int> &target)
+{
+    // target is a vector of class labels that take values from [0, 1, ..., nclass - 1]
+    // The i-th element of target is the class label for observation i
+    int nobs = prev_layer_data.col_size();
+    int ncls = prev_layer_data.row_size();
+
+    if (target.size() != nobs)
+    {
+        throw std::invalid_argument("[gctl::multi_entropy] Target data have incorrect dimension.");
+    }
+
+    // Compute the derivative of the input of this layer
+    // L = -log(phat[y])
+    // in = phat
+    // d(L) / d(in) = [0, 0, ..., -1/phat[y], 0, ..., 0]
+    der_in_.resize(ncls, nobs);
+    der_in_.assign_all(0.0);
+
+    int i;
+#pragma omp parallel for private (i) schedule(guided)
+    for (i = 0; i < nobs; i++)
+    {
+        der_in_[target[i]][i] = -1.0/prev_layer_data[target[i]][i];
+    }
+    return;
+}
+
+double gctl::multi_entropy::loss_value() const
+{
+    return loss_;
+}
+
+std::string gctl::multi_entropy::get_output_name() const
+{
+    return "MultiClassEntropy";
+}
+
+gctl::olayer_type_e gctl::multi_entropy::get_output_type() const
+{
+    return MultiClassEntropy;
+}
\ No newline at end of file
diff --git a/lib/dnn/olayer_multientropy.h b/lib/dnn/olayer_multientropy.h
new file mode 100644
index 0000000..cf96504
--- /dev/null
+++ b/lib/dnn/olayer_multientropy.h
@@ -0,0 +1,76 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_OLAYER_MULTIENTROPY_H
+#define _GCTL_DNN_OLAYER_MULTIENTROPY_H
+
+#include "olayer.h"
+
+namespace gctl
+{
+    class multi_entropy : public dnn_olayer
+    {
+    public:
+        multi_entropy();
+        virtual ~multi_entropy();
+
+        // Check the format of target data, e.g. in classification problems the
+        // target data should be binary (either 0 or 1)
+        void check_target_data(const matrix<double> &target);
+
+        // Another type of target data where each element is a class label
+        // This version may not be sensible for regression tasks, so by default
+        // we raise an exception
+        void check_target_data(const array<int> &target);
+
+        // A combination of the forward stage and the back-propagation stage for the output layer
+        // The computed derivative of the input should be stored in this layer, and can be retrieved by
+        // the backward_propagation_data() function
+        void evaluation(const matrix<double> &prev_layer_data, const matrix<double> &target);
+
+        // Another type of target data where each element is a class label
+        // This version may not be sensible for regression tasks, so by default
+        // we raise an exception
+        void evaluation(const matrix<double> &prev_layer_data, const array<int> &target);
+
+        // Return the loss function value after the evaluation
+        // This function can be assumed to be called after evaluate(), so that it can make use of the
+        // intermediate result to save some computation
+        double loss_value() const;
+
+        // Return the output layer name. It is used to export the NN model.
+        std::string get_output_name() const;
+
+        // Return the output layer type. It is used to export the NN model.
+        olayer_type_e get_output_type() const;
+    
+    private:
+        double loss_;
+    };
+}
+
+#endif // _GCTL_DNN_OLAYER_MULTIENTROPY_H
\ No newline at end of file
diff --git a/lib/dnn/olayer_rmse.cpp b/lib/dnn/olayer_rmse.cpp
new file mode 100644
index 0000000..6e84bfb
--- /dev/null
+++ b/lib/dnn/olayer_rmse.cpp
@@ -0,0 +1,77 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#include "olayer_rmse.h"
+
+gctl::rmse::rmse() {}
+
+gctl::rmse::~rmse() {}
+
+void gctl::rmse::evaluation(const matrix<double> &prev_layer_data, const matrix<double> &target)
+{
+    int nobs = prev_layer_data.col_size();
+    int nvar = prev_layer_data.row_size();
+
+    if (target.col_size() != nobs || target.row_size() != nvar)
+    {
+        throw std::invalid_argument("[gctl::rmse] Target data have incorrect dimension.");
+    }
+    
+    // Compute the derivative of the input of this layer
+    // L = 0.5 * ||yhat - y||^2
+    // in = yhat
+    // d(L) / d(in) = yhat - y
+    der_in_.resize(nvar, nobs);
+
+    int i, j;
+#pragma omp parallel for private (i, j) schedule(guided)
+    for (i = 0; i < nvar; i++)
+    {
+        for (j = 0; j < nobs; j++)
+        {
+            der_in_[i][j] = prev_layer_data[i][j] - target[i][j];
+        }
+    }
+    return;
+}
+
+double gctl::rmse::loss_value() const
+{
+    // L = 0.5 * ||yhat - y||^2
+    double n = der_in_.norm(L2);
+    return 0.5*power2(n)/der_in_.col_size();
+}
+
+std::string gctl::rmse::get_output_name() const
+{
+    return "RegressionMSE";
+}
+
+gctl::olayer_type_e gctl::rmse::get_output_type() const
+{
+    return RegressionMSE;
+}
\ No newline at end of file
diff --git a/lib/dnn/olayer_rmse.h b/lib/dnn/olayer_rmse.h
new file mode 100644
index 0000000..6770663
--- /dev/null
+++ b/lib/dnn/olayer_rmse.h
@@ -0,0 +1,59 @@
+/********************************************************
+ *  ██████╗  ██████╗████████╗██╗
+ * ██╔════╝ ██╔════╝╚══██╔══╝██║
+ * ██║  ███╗██║        ██║   ██║
+ * ██║   ██║██║        ██║   ██║
+ * ╚██████╔╝╚██████╗   ██║   ███████╗
+ *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
+ * Geophysical Computational Tools & Library (GCTL)
+ *
+ * Copyright (c) 2022  Yi Zhang (yizhang-geo@zju.edu.cn)
+ *
+ * GCTL is distributed under a dual licensing scheme. You can redistribute 
+ * it and/or modify it under the terms of the GNU Lesser General Public 
+ * License as published by the Free Software Foundation, either version 2 
+ * of the License, or (at your option) any later version. You should have 
+ * received a copy of the GNU Lesser General Public License along with this 
+ * program. If not, see <http://www.gnu.org/licenses/>.
+ * 
+ * If the terms and conditions of the LGPL v.2. would prevent you from using 
+ * the GCTL, please consider the option to obtain a commercial license for a 
+ * fee. These licenses are offered by the GCTL's original author. As a rule, 
+ * licenses are provided "as-is", unlimited in time for a one time fee. Please 
+ * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
+ * to include some description of your company and the realm of its activities. 
+ * Also add information on how to contact you by electronic and paper mail.
+ ******************************************************/
+
+#ifndef _GCTL_DNN_OLAYER_RMSE_H
+#define _GCTL_DNN_OLAYER_RMSE_H
+
+#include "olayer.h"
+
+namespace gctl
+{
+    class rmse : public dnn_olayer
+    {
+    public:
+        rmse();
+        virtual ~rmse();
+
+        // A combination of the forward stage and the back-propagation stage for the output layer
+        // The computed derivative of the input should be stored in this layer, and can be retrieved by
+        // the backward_propagation_data() function
+        void evaluation(const matrix<double> &prev_layer_data, const matrix<double> &target);
+
+        // Return the loss function value after the evaluation
+        // This function can be assumed to be called after evaluate(), so that it can make use of the
+        // intermediate result to save some computation
+        double loss_value() const;
+
+        // Return the output layer name. It is used to export the NN model.
+        std::string get_output_name() const;
+
+        // Return the output layer type. It is used to export the NN model.
+        olayer_type_e get_output_type() const;
+    };
+}
+
+#endif // _GCTL_DNN_OLAYER_RMSE_H
\ No newline at end of file