From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Fri, 30 Jun 2000 20:54:43 +0000 (+0000)
Subject: experimental LAPACK support
X-Git-Tag: v8.0.0~20308
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3c9d01a10e910d99e7263f451387f8e55965ee18;p=dealii.git

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 f69230e76e..1cc3f7b6cc 100644
--- 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