const Vector<TYPEVEC>&,
const Vector<TYPERES>&) const;
-// Experimental code
-
-#ifdef HAVE_LIBLAPACK_AND_REALLY_USE_IT_HERE
-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)
-{
- double* val = const_cast<double*> (data());
-
- Assert (val != 0, ExcEmptyMatrix());
-
- unsigned int dim_range = n_cols ();
- unsigned int dim_image = n_rows ();
-
- Assert (dim_range == dim_image, ExcNotQuadratic());
- Assert (dim_range == M.n_cols(),
- ExcDimensionMismatch(dim_range,M.n_cols()));
- Assert (dim_image == M.n_rows(),
- ExcDimensionMismatch(dim_image,M.n_rows()));
-
- 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];
-
- double* matrix = new double[dim_range*dim_range];
- std::copy (M.data(), M.data()+dim_image*dim_range, matrix);
-
-
- int erg = dgelss_ (&dim_range, &dim_range, &dim_range,
- matrix, &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 << std::endl;
- if (rank<(int)dim_range)
- deallog << "rank deficiency " << rank << std::endl;
- delete[] work;
- delete[] s;
- delete[] matrix;
-}
-
-#endif