From 70217103f9f83961dd269972dcf69422fb8c5a5f Mon Sep 17 00:00:00 2001 From: guido Date: Fri, 30 Jun 2000 20:54:43 +0000 Subject: [PATCH] experimental LAPACK support git-svn-id: https://svn.dealii.org/trunk@3115 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/source/full_matrix.double.cc | 70 ++++++++++++++++++++++++ 1 file changed, 70 insertions(+) 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 +#include #define TYPEMAT double @@ -61,3 +62,72 @@ template double FullMatrix::least_squares(Vector&, Vector::residual(Vector&, const Vector&, const Vector&) 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::invert (const FullMatrix &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 -- 2.39.5