add svd by lapack, just for debugging

[[Imported from SVN: r1602]]

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);
+        }
+
+    }
+}

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