diff --git a/src/cotmatrix.cpp b/src/cotmatrix.cpp
index c888e5c..fa1b77f 100644
--- a/src/cotmatrix.cpp
+++ b/src/cotmatrix.cpp
@@ -1,4 +1,13 @@
 #include "cotmatrix.h"
+#include <Eigen/Sparse>
+#include <math.h>
+#include <ppl.h>
+#include <concurrent_vector.h>
+inline void vertex_index(const int i, int& il, int& ir)
+{
+  il = i - 1 < 0 ? 2 : i - 1;
+  ir = i + 1 > 2 ? 0 : i + 1;
+}
 
 void cotmatrix(
   const Eigen::MatrixXd & l,
@@ -6,5 +15,45 @@ void cotmatrix(
   Eigen::SparseMatrix<double> & L)
 {
   // Add your code here
-}
+  typedef Eigen::Triplet<double> tuple;
+  Concurrency::concurrent_vector<tuple> tuple_list;
+  const int n = F.rows();
+  tuple_list.reserve(12 * n);
+  Concurrency::parallel_for(size_t(0),size_t(n), [&l,&F,n,&tuple_list](const int m)
+  {
+    int il, ir;
+    int v[3];
+    double e[3];
+    double s = 0;
+    for (int i = 0; i < 3; i++)
+    {
+      v[i] = F(m, i);
+      e[i] = l(m, i);
+      s += 0.5*e[i];
+    }
+    double A = sqrt(s*(s - e[0])*(s - e[1])*(s - e[2]));
+    double d = 2 * A / (e[0] * e[1] * e[2]);
+
+    double cosine[3];
+    double sine[3];
+    double cot[3];
+
+    for (int i = 0; i < 3; i++)
+    {
+      vertex_index(i, il, ir);
+      cosine[i] = (e[il] * e[il] + e[ir] * e[ir] - e[i] * e[i]) / (2 * e[il] * e[ir]);
+      sine[i] = e[i] * d;
+      cot[i] = cosine[i] / sine[i];
+    }
 
+    for (int i = 0; i < 3; i++)
+    {
+      vertex_index(i, il, ir);
+      tuple_list.push_back(tuple(v[il], v[ir], 0.5*cot[i]));
+      tuple_list.push_back(tuple(v[ir], v[il], 0.5*cot[i]));
+      tuple_list.push_back(tuple(v[il], v[il], -0.5*cot[i]));
+      tuple_list.push_back(tuple(v[ir], v[ir], -0.5*cot[i]));
+    } 
+  }, Concurrency::auto_partitioner());
+  L.setFromTriplets(tuple_list.begin(), tuple_list.end());
+}
\ No newline at end of file
diff --git a/src/massmatrix.cpp b/src/massmatrix.cpp
index 577e90f..01745e8 100644
--- a/src/massmatrix.cpp
+++ b/src/massmatrix.cpp
@@ -1,10 +1,22 @@
 #include "massmatrix.h"
+#include <ppl.h>
 
 void massmatrix(
   const Eigen::MatrixXd & l,
   const Eigen::MatrixXi & F,
   Eigen::DiagonalMatrix<double,Eigen::Dynamic> & M)
 {
-  // Add your code here
-}
-
+  const int n = F.maxCoeff();
+  const int m = F.rows();
+  Concurrency::parallel_for(size_t(0),size_t(n),[m,&l,&F,&M](const int i)
+  {
+    double sum = 0;
+    for (int j = 0; j < m; j++)
+      if (i == F(j, 0) || i == F(j, 1) || i == F(j, 2))
+      {
+        double s = 0.5*(l(j, 0) + l(j, 1) + l(j, 2));
+        sum += sqrt(s*(s - l(j, 0))*(s - l(j, 1))*(s - l(j, 2)));
+      }
+    M.diagonal()(i) = sum / 3;
+  }, Concurrency::auto_partitioner());
+}
\ No newline at end of file
diff --git a/src/smooth.cpp b/src/smooth.cpp
index 6f72147..c3f44a6 100644
--- a/src/smooth.cpp
+++ b/src/smooth.cpp
@@ -1,4 +1,11 @@
 #include "smooth.h"
+#include "cotmatrix.h"
+#include "massmatrix.h"
+#include <ppl.h>
+#include <iostream>
+#include <Eigen/Sparse>
+#include <concurrent_vector.h>
+#include <igl\edge_lengths.h>
 
 void smooth(
     const Eigen::MatrixXd & V,
@@ -8,5 +15,32 @@ void smooth(
     Eigen::MatrixXd & U)
 {
   // Replace with your code
-  U = G;
-}
+  int face_num = F.rows();
+  int vertex_num = F.maxCoeff() + 1;
+  Eigen::MatrixXd l(face_num, 3); l.setZero();
+  igl::edge_lengths(V, F, l);
+
+  Eigen::DiagonalMatrix<double, Eigen::Dynamic> M(vertex_num);
+  Eigen::SparseMatrix<double>SM(vertex_num, vertex_num); SM.setZero();
+  massmatrix(l, F, M);
+
+  auto diag_to_sparse = [&M, &SM](Eigen::DiagonalMatrix<double, Eigen::Dynamic> M,
+    Eigen::SparseMatrix<double>&SM)
+  {
+    typedef Eigen::Triplet<double> tuple;
+    Concurrency::concurrent_vector<tuple> tuple_list;
+    const int n = M.rows();
+    Concurrency::parallel_for(size_t(0), size_t(n), [&M, &SM, &tuple_list](const int m)
+    {
+      tuple_list.push_back(tuple(m, m, M.diagonal().array()[m]));
+    }, Concurrency::auto_partitioner());
+    SM.setFromTriplets(tuple_list.begin(), tuple_list.end());
+  };
+
+  diag_to_sparse(M, SM);
+  Eigen::SparseMatrix<double>L(vertex_num, vertex_num); L.setZero();
+  cotmatrix(l, F, L);
+
+  Eigen::SimplicialCholesky<Eigen::SparseMatrix<double>> solver(SM - lambda*L);
+  U = solver.solve(M*G);
+}
\ No newline at end of file