diff --git a/linearsolver.cc b/linearsolver.cc
index d66153a4c895c2fe8a3589ada15c9f4c2400daae..d8dac993b5ed4cdbcf9e467f550b1b9533267f9f 100644
--- a/linearsolver.cc
+++ b/linearsolver.cc
@@ -33,3 +33,30 @@ void linearSolver(const FieldMatrix<double,6,12>& A,
     for (int i=0; i<M; i++)
         x[i] = X(i);
 }
+
+// Compute the svd using lapack++
+void lapackSVD(const FieldMatrix<double,3,3>& A,
+               FieldMatrix<double,3,3>& U,
+               FieldVector<double,3>& sigma,
+               FieldMatrix<double,3,3>& VT)
+{
+    LaGenMatDouble lpA(3,3), lpU(3,3), lpVT(3,3);
+    LaVectorDouble lpSigma(3);
+
+    for (int i=0; i<3; i++)
+        for (int j=0; j<3; j++)
+            lpA(i,j) = A[i][j];
+
+    LaSVD_IP(lpA, lpSigma, lpU, lpVT);
+
+    for (int i=0; i<3; i++) {
+
+        sigma[i] = lpSigma(i);
+
+        for (int j=0; j<3; j++) {
+            U[i][j]  = lpU(i,j);
+            VT[i][j] = lpVT(i,j);
+        }
+
+    }
+}
diff --git a/linearsolver.hh b/linearsolver.hh
index e42607ac9e22bdd5b64c698decfd5290df7a3a7b..b66376e1105d70e35dd9583246a2914282feb05c 100644
--- a/linearsolver.hh
+++ b/linearsolver.hh
@@ -7,4 +7,9 @@ void linearSolver(const Dune::FieldMatrix<double,6,12>& A,
                   Dune::FieldVector<double,12>& x,
                   const Dune::FieldVector<double,6>& b);
 
+void lapackSVD(const Dune::FieldMatrix<double,3,3>& A,
+               Dune::FieldMatrix<double,3,3>& U,
+               Dune::FieldVector<double,3>& sigma,
+              Dune::FieldMatrix<double,3,3>& VT);
+
 #endif