#include <deal.II/base/table.h>
#include <deal.II/lac/lapack_support.h>
#include <deal.II/lac/vector_memory.h>
+#include <deal.II/base/thread_management.h>
#include <memory>
#include <vector>
*/
mutable std::vector<number> work;
+ /**
+ * Integer working array used for some LAPACK functions.
+ */
+ mutable std::vector<int> iwork;
+
/**
* The vector storing the permutations applied for pivoting in the LU-
* factorization.
* <i>USV<sup>T</sup></i>.
*/
std::shared_ptr<LAPACKFullMatrix<number> > svd_vt;
+
+ /**
+ * Thread mutex.
+ */
+ mutable Threads::Mutex mutex;
};
const Vector<number> &v,
const bool adding) const
{
+ Threads::Mutex::ScopedLock lock (mutex);
+
const int mm = this->n_rows();
const int nn = this->n_cols();
const number alpha = 1.;
const Vector<number> &v,
const bool adding) const
{
+ Threads::Mutex::ScopedLock lock (mutex);
+
const int mm = this->n_rows();
const int nn = this->n_cols();
const number alpha = 1.;
template <typename number>
number LAPACKFullMatrix<number>::norm(const char type) const
{
+ Threads::Mutex::ScopedLock lock (mutex);
+
Assert (state == LAPACKSupport::matrix ||
state == LAPACKSupport::inverse_matrix,
ExcMessage("norms can be called in matrix state only."));
const int N = this->n_cols();
const int M = this->n_rows();
const number *values = &this->values[0];
- std::vector<number> w;
if (property == symmetric)
{
const int lda = std::max(1,N);
const int lwork = (type == 'I' || type == 'O') ?
std::max(1,N) :
0;
- w.resize(lwork);
- return lansy (&type, &LAPACKSupport::L, &N, values, &lda, &w[0]);
+ work.resize(lwork);
+ return lansy (&type, &LAPACKSupport::L, &N, values, &lda, &work[0]);
}
else
{
const int lwork = (type == 'I') ?
std::max(1,M) :
0;
- w.resize(lwork);
- return lange (&type, &M, &N, values, &lda, &w[0]);
+ work.resize(lwork);
+ return lange (&type, &M, &N, values, &lda, &work[0]);
}
}
number
LAPACKFullMatrix<number>::reciprocal_condition_number(const number a_norm) const
{
+ Threads::Mutex::ScopedLock lock (mutex);
Assert(state == cholesky, ExcState(state));
number rcond = 0.;
const number *values = &this->values[0];
int info = 0;
const int lda = std::max(1,N);
- std::vector<number> w(3*N);
- std::vector<int> iw(N);
+ work.resize(3*N);
+ iwork.resize(N);
// use the same uplo as in Cholesky
pocon (&LAPACKSupport::L, &N, values, &lda,
&a_norm, &rcond,
- &w[0], &iw[0], &info);
+ &work[0], &iwork[0], &info);
Assert(info >= 0, ExcInternalError());
const char *const uplo(&U);
const char *const range(&V);
const int *const dummy(&one);
- std::vector<int> iwork(static_cast<size_type> (5*nn));
+ iwork.resize(static_cast<size_type> (5*nn));
std::vector<int> ifail(static_cast<size_type> (nn));