]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make Householder class compile with complex numbers
authorMartin Kronbichler <martin.kronbichler@rub.de>
Tue, 20 May 2025 12:38:52 +0000 (14:38 +0200)
committerMartin Kronbichler <martin.kronbichler@rub.de>
Tue, 20 May 2025 12:38:52 +0000 (14:38 +0200)
include/deal.II/lac/householder.h

index c0a64e3b70ba2d86b32e075c1dc9ce716eeebef1..0f5051b36c51835fbb367b054ff8e6c3fb7178b0 100644 (file)
@@ -115,7 +115,7 @@ public:
    * return.
    */
   template <typename number2>
-  double
+  number2
   least_squares(Vector<number2> &dst, const Vector<number2> &src) const;
 
   /**
@@ -153,7 +153,7 @@ private:
   /**
    * Storage that is internally used for the Householder transformation.
    */
-  FullMatrix<double> storage;
+  FullMatrix<number> storage;
 };
 
 /** @} */
@@ -190,7 +190,11 @@ Householder<number>::initialize(const FullMatrix<number2> &M)
       if (std::fabs(sigma) < 1.e-15)
         return;
 
-      number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
+      number2 s;
+      if constexpr (numbers::NumberTraits<number2>::is_complex)
+        s = storage(j, j).real() < 0 ? std::sqrt(sigma) : -std::sqrt(sigma);
+      else
+        s = storage(j, j) < 0 ? std::sqrt(sigma) : -std::sqrt(sigma);
       //
       number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));
 
@@ -232,7 +236,7 @@ Householder<number>::Householder(const FullMatrix<number2> &M)
 
 template <typename number>
 template <typename number2>
-double
+number2
 Householder<number>::least_squares(Vector<number2>       &dst,
                                    const Vector<number2> &src) const
 {
@@ -242,10 +246,7 @@ Householder<number>::least_squares(Vector<number2>       &dst,
 
   const size_type m = storage.m(), n = storage.n();
 
-  GrowingVectorMemory<Vector<number2>>            mem;
-  typename VectorMemory<Vector<number2>>::Pointer aux(mem);
-  aux->reinit(src, true);
-  *aux = src;
+  Vector<number2> aux(src);
   // m > n, m = src.n, n = dst.n
 
   // Multiply Q_n ... Q_2 Q_1 src
@@ -253,22 +254,22 @@ Householder<number>::least_squares(Vector<number2>       &dst,
   for (size_type j = 0; j < n; ++j)
     {
       // sum = v_i^T dst
-      number2 sum = diagonal[j] * (*aux)(j);
+      number2 sum = diagonal[j] * aux(j);
       for (size_type i = j + 1; i < m; ++i)
-        sum += static_cast<number2>(storage(i, j)) * (*aux)(i);
+        sum += static_cast<number2>(storage(i, j)) * aux(i);
       // dst -= v * sum
-      (*aux)(j) -= sum * diagonal[j];
+      aux(j) -= sum * diagonal[j];
       for (size_type i = j + 1; i < m; ++i)
-        (*aux)(i) -= sum * static_cast<number2>(storage(i, j));
+        aux(i) -= sum * static_cast<number2>(storage(i, j));
     }
   // Compute norm of residual
   number2 sum = 0.;
   for (size_type i = n; i < m; ++i)
-    sum += (*aux)(i) * (*aux)(i);
+    sum += aux(i) * aux(i);
   AssertIsFinite(sum);
 
   // Compute solution
-  storage.backward(dst, *aux);
+  storage.backward(dst, aux);
 
   return std::sqrt(sum);
 }
@@ -287,10 +288,9 @@ Householder<number>::least_squares(BlockVector<number2>       &dst,
 
   const size_type m = storage.m(), n = storage.n();
 
-  GrowingVectorMemory<BlockVector<number2>>            mem;
-  typename VectorMemory<BlockVector<number2>>::Pointer aux(mem);
-  aux->reinit(src, true);
-  *aux = src;
+  BlockVector<number2> aux;
+  aux.reinit(src, true);
+  aux = src;
   // m > n, m = src.n, n = dst.n
 
   // Multiply Q_n ... Q_2 Q_1 src
@@ -298,30 +298,27 @@ Householder<number>::least_squares(BlockVector<number2>       &dst,
   for (size_type j = 0; j < n; ++j)
     {
       // sum = v_i^T dst
-      number2 sum = diagonal[j] * (*aux)(j);
+      number2 sum = diagonal[j] * aux(j);
       for (size_type i = j + 1; i < m; ++i)
-        sum += storage(i, j) * (*aux)(i);
+        sum += storage(i, j) * aux(i);
       // dst -= v * sum
-      (*aux)(j) -= sum * diagonal[j];
+      aux(j) -= sum * diagonal[j];
       for (size_type i = j + 1; i < m; ++i)
-        (*aux)(i) -= sum * storage(i, j);
+        aux(i) -= sum * storage(i, j);
     }
   // Compute norm of residual
   number2 sum = 0.;
   for (size_type i = n; i < m; ++i)
-    sum += (*aux)(i) * (*aux)(i);
+    sum += *aux(i) * aux(i);
   AssertIsFinite(sum);
 
-  // backward works for
-  // Vectors only, so copy
-  // them before
+  // backward works for Vectors only, so copy them before
   Vector<number2> v_dst, v_aux;
   v_dst = dst;
-  v_aux = *aux;
+  v_aux = aux;
   // Compute solution
   storage.backward(v_dst, v_aux);
-  // copy the result back
-  // to the BlockVector
+  // copy the result back to the BlockVector
   dst = v_dst;
 
   return std::sqrt(sum);

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.