]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make LapackFullMatrix also work with std::complex data types. 5395/head
authorWolfgang Bangerth <bangerth@colostate.edu>
Tue, 31 Oct 2017 21:48:10 +0000 (15:48 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Sat, 4 Nov 2017 02:00:16 +0000 (20:00 -0600)
This is slightly complicated by the fact that the LAPACK functions sometimes have
a different signature for complex data types. This patch does not address this
-- these functions will simply fail if the signature is wrong. However, at least
we can compile things with this patch.

include/deal.II/lac/lapack_full_matrix.h
source/lac/lapack_full_matrix.cc

index 8228c234262a69f6c14a7be96de197c6ad8cd08f..7d75ba796b3f6d4353019f2e2ad2e1a10cb8ac55 100644 (file)
@@ -117,10 +117,11 @@ public:
 
   /**
    * This operator assigns a scalar to a matrix. To avoid confusion with
-   * constructors, zero is the only value allowed for <tt>d</tt>
+   * constructors, zero (when cast to the @p number type) is the only
+   * value allowed for <tt>d</tt>
    */
   LAPACKFullMatrix<number> &
-  operator = (const double d);
+  operator = (const number d);
 
   /**
    * This operator multiplies all entries by a fixed factor.
index c6918d2a3bdbe567b026b36f4b9f462eed375ffc..8d00a95018872b65f97b4f89cfa77bdefae4d8a3 100644 (file)
@@ -123,10 +123,10 @@ LAPACKFullMatrix<number>::operator = (const SparseMatrix<number2> &M)
 
 template <typename number>
 LAPACKFullMatrix<number> &
-LAPACKFullMatrix<number>::operator = (const double d)
+LAPACKFullMatrix<number>::operator = (const number d)
 {
   (void)d;
-  Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
+  Assert (d==number(0), ExcScalarAssignmentOnlyForZeroValue());
 
   if (this->n_elements() != 0)
     this->reset_values();
@@ -718,7 +718,7 @@ LAPACKFullMatrix<number>::compute_svd()
   AssertThrow (info==0, LAPACKSupport::ExcErrorCode("gesdd", info));
   // Resize the work array. Add one to the size computed by LAPACK to be on
   // the safe side.
-  lwork = static_cast<int>(work[0] + 1);
+  lwork = static_cast<int>(std::fabs(work[0]) + 1);
 
   work.resize(lwork);
   // Do the actual SVD.
@@ -743,11 +743,11 @@ LAPACKFullMatrix<number>::compute_inverse_svd(const double threshold)
 
   Assert (state==LAPACKSupport::svd, ExcState(state));
 
-  const double lim = wr[0]*threshold;
+  const double lim = std::fabs(wr[0])*threshold;
   for (size_type i=0; i<wr.size(); ++i)
     {
-      if (wr[i] > lim)
-        wr[i] = 1./wr[i];
+      if (std::fabs(wr[i]) > lim)
+        wr[i] = number(1.)/wr[i];
       else
         wr[i] = 0.;
     }
@@ -911,7 +911,7 @@ LAPACKFullMatrix<number>::compute_eigenvalues(const bool right,
   Assert (info == 0, ExcInternalError());
   // Allocate working array according to suggestion (same strategy as was
   // noted in compute_svd).
-  lwork = static_cast<int>(work[0] + 1);
+  lwork = static_cast<int>(std::fabs(work[0]) + 1);
 
   // resize workspace array
   work.resize((size_type ) lwork);
@@ -984,7 +984,7 @@ LAPACKFullMatrix<number>::compute_eigenvalues_symmetric(const number        lowe
   Assert (info == 0, ExcInternalError());
   // Allocate working array according to suggestion (same strategy as was noted in
   // compute_svd).
-  lwork = static_cast<int>(work[0] + 1);
+  lwork = static_cast<int>(std::fabs(work[0]) + 1);
   work.resize(static_cast<size_type> (lwork));
 
   // Finally compute the eigenvalues.
@@ -1075,7 +1075,7 @@ LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
   Assert (info == 0, ExcInternalError());
   // Allocate working array according to suggestion (same strategy as was
   // noted in compute_svd).
-  lwork = static_cast<int>(work[0] + 1);
+  lwork = static_cast<int>(std::fabs(work[0]) + 1);
 
   // resize workspace arrays
   work.resize(static_cast<size_type> (lwork));
@@ -1157,7 +1157,7 @@ LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric (
   Assert (info == 0, ExcInternalError());
   // Allocate working array according to suggestion (same strategy as was
   // noted in compute_svd).
-  lwork = static_cast<int>(work[0] + 1);
+  lwork = static_cast<int>(std::fabs(work[0]) + 1);
 
   // resize workspace array
   work.resize(static_cast<size_type>(lwork));

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.