experimental LAPACK support

git-svn-id: https://svn.dealii.org/trunk@3115 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/lac/source/full_matrix.double.cc b/deal.II/lac/source/full_matrix.double.cc
index f69230e76ec16ddb45449283d9b07485c759ec1b..1cc3f7b6cc50f2d24eb17304f6df0a431aad4926 100644 (file)
--- a/deal.II/lac/source/full_matrix.double.cc
+++ b/deal.II/lac/source/full_matrix.double.cc
@@ -13,6 +13,7 @@
 
 
 #include <lac/full_matrix.templates.h>
+#include <base/logstream.h>
 
 #define TYPEMAT double
 
@@ -61,3 +62,72 @@ template double FullMatrix<TYPEMAT>::least_squares(Vector<TYPEVEC>&, Vector<TYPE
 
 template double FullMatrix<TYPEMAT>::residual(Vector<TYPEVEC>&, const Vector<TYPEVEC>&, const Vector<TYPERES>&) const;
 
+// Experimental code
+
+#ifdef HAVE_LIBLAPACK
+extern "C" int dgels_ (const char* trans,
+                      const unsigned int* M, const unsigned int* N,
+                      const unsigned int* NRHS,
+                      double* A, const unsigned int* LDA,
+                      double* B, const unsigned int* LDB,
+                      double* WORK, const unsigned int* LWORK,
+                      int* INFO);
+extern "C" int dgelss_ (const unsigned int* M, const unsigned int* N,
+                       const unsigned int* NRHS,
+                       double* A, const unsigned int* LDA,
+                       double* B, const unsigned int* LDB,
+                       double* S, const double* RCOND,
+                       int* RANK,
+                       double* WORK, const unsigned int* LWORK,
+                       int* INFO);
+
+
+template<>
+void
+FullMatrix<double>::invert (const FullMatrix<double> &M)
+{
+  Assert (val != 0, ExcEmptyMatrix());
+  
+  Assert (dim_range == dim_image, ExcNotQuadratic());
+  Assert (dim_range == M.dim_range,
+          ExcDimensionMismatch(dim_range,M.dim_range));
+  Assert (dim_image == M.dim_image,
+         ExcDimensionMismatch(dim_image,M.dim_image));
+
+  clear();
+  diagadd(1.);
+
+  const unsigned int lwork = 10*dim_range*dim_range;
+  double* work = new double[lwork];
+  int info;
+  
+  const char* trans = "N";
+  int rank;
+  const double rcond=-1.;
+  double* s = new double[dim_range];
+  
+  int erg = dgelss_ (&dim_range, &dim_range, &dim_range,
+                    M.val, &dim_range,
+                    val, &dim_range,
+                    s, &rcond, &rank,
+                    work, &lwork,
+                    &info);
+//    int erg = dgels_ (trans, &dim_range, &dim_range, &dim_range,
+//                 M.val, &dim_range,
+//                 val, &dim_range,
+//                 work, &lwork,
+//                 &info);
+
+  double condition = s[0]/s[dim_range-1];
+  
+  if (info!=0)
+    deallog << "inverting error " << info << ' ' << erg << endl;
+  if (rank<(int)dim_range)
+    deallog << "rank deficiency " << rank << endl;
+  if (condition > 1.e6)
+    deallog << "Condition " <<  condition << endl;
+  delete[] work;
+  delete[] s;
+}
+
+#endif