From f14f34f5af7859e315ecdca2c7d1bc34a9f93e76 Mon Sep 17 00:00:00 2001 From: prill Date: Tue, 16 May 2006 07:48:47 +0000 Subject: [PATCH] Fix for LAPACK without workspace query mode. git-svn-id: https://svn.dealii.org/trunk@13116 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/source/lapack_full_matrix.cc | 33 +++++++++++++++++++++--- 1 file changed, 30 insertions(+), 3 deletions(-) diff --git a/deal.II/lac/source/lapack_full_matrix.cc b/deal.II/lac/source/lapack_full_matrix.cc index 092d912642..4af23fbaaa 100644 --- a/deal.II/lac/source/lapack_full_matrix.cc +++ b/deal.II/lac/source/lapack_full_matrix.cc @@ -161,9 +161,31 @@ LAPACKFullMatrix::compute_eigenvalues() wi.resize(nn); number* values = const_cast (this->data()); - int info; - int lwork = -1; + int info = 0; + int lwork = 1; + // Optimal workspace query: + + // The LAPACK routine DGEEV requires + // a sufficient large workspace variable, + // minimum requirement is + // work.size>=4*nn. + // However, to improve performance, a + // somewhat larger workspace may be needed. + + // SOME implementations of the LAPACK routine + // provide a workspace query call, + // info:=0, lwork:=-1 + // which returns an optimal value for the + // size of the workspace array + // (the PETSc 2.3.0 implementation does NOT + // provide this functionality!). + + // define the HAVE_LIBLAPACK_NOQUERYMODE flag to + // disable the workspace query. +#ifndef HAVE_LIBLAPACK_NOQUERYMODE + lwork = -1; work.resize(1); + geev(&N, &N, &nn, values, &nn, &wr[0], &wi[0], 0, &one, 0, &one, @@ -174,11 +196,15 @@ LAPACKFullMatrix::compute_eigenvalues() // everything else would not be // acceptable. Assert (info == 0, ExcInternalError()); - // Allocate working array according // to suggestion. lwork = (int) (work[0]+.1); +#else + lwork = 4*nn; // no query mode +#endif + // resize workspace array work.resize((unsigned int) lwork); + // Finally compute the eigenvalues. geev(&N, &N, &nn, values, &nn, &wr[0], &wi[0], @@ -187,6 +213,7 @@ LAPACKFullMatrix::compute_eigenvalues() // Negative return value implies a // wrong argument. This should be // internal. + Assert (info >=0, ExcInternalError()); //TODO:[GK] What if the QR method fails? if (info != 0) -- 2.39.5