householder.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_householder_h
#define dealii_householder_h


#include <deal.II/base/config.h>

#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/vector_memory.h>

#include <cmath>
#include <vector>

DEAL_II_NAMESPACE_OPEN


// forward declarations
template <typename number>
class Vector;


template <typename number>
class Householder
{
public:
  using size_type = types::global_dof_index;

  Householder() = default;

  template <typename number2>
  Householder(const FullMatrix<number2> &A);

  template <typename number2>
  void
  initialize(const FullMatrix<number2> &A);

  template <typename number2>
  double
  least_squares(Vector<number2> &dst, const Vector<number2> &src) const;

  template <typename number2>
  double
  least_squares(BlockVector<number2> &      dst,
                const BlockVector<number2> &src) const;

  template <class VectorType>
  void
  vmult(VectorType &dst, const VectorType &src) const;

  template <class VectorType>
  void
  Tvmult(VectorType &dst, const VectorType &src) const;


private:
  std::vector<number> diagonal;

  FullMatrix<double> storage;
};

#ifndef DOXYGEN
/*-------------------------Inline functions -------------------------------*/

// QR-transformation cf. Stoer 1 4.8.2 (p. 191)

template <typename number>
template <typename number2>
void
Householder<number>::initialize(const FullMatrix<number2> &M)
{
  const size_type m = M.n_rows(), n = M.n_cols();
  storage.reinit(m, n);
  storage.fill(M);
  Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
  diagonal.resize(m);

  // m > n, src.n() = m
  Assert(storage.n_cols() <= storage.n_rows(),
         ExcDimensionMismatch(storage.n_cols(), storage.n_rows()));

  for (size_type j = 0; j < n; ++j)
    {
      number2   sigma = 0;
      size_type i;
      // sigma = ||v||^2
      for (i = j; i < m; ++i)
        sigma += storage(i, j) * storage(i, j);
      // We are ready if the column is
      // empty. Are we?
      if (std::fabs(sigma) < 1.e-15)
        return;

      number2 s = (storage(j, j) < 0) ? std::sqrt(sigma) : -std::sqrt(sigma);
      //
      number2 beta = std::sqrt(1. / (sigma - s * storage(j, j)));

      // Make column j the Householder
      // vector, store first entry in
      // diagonal
      diagonal[j]   = beta * (storage(j, j) - s);
      storage(j, j) = s;

      for (i = j + 1; i < m; ++i)
        storage(i, j) *= beta;


      // For all subsequent columns do
      // the Householder reflection
      for (size_type k = j + 1; k < n; ++k)
        {
          number2 sum = diagonal[j] * storage(j, k);
          for (i = j + 1; i < m; ++i)
            sum += storage(i, j) * storage(i, k);

          storage(j, k) -= sum * this->diagonal[j];
          for (i = j + 1; i < m; ++i)
            storage(i, k) -= sum * storage(i, j);
        }
    }
}



template <typename number>
template <typename number2>
Householder<number>::Householder(const FullMatrix<number2> &M)
{
  initialize(M);
}



template <typename number>
template <typename number2>
double
Householder<number>::least_squares(Vector<number2> &      dst,
                                   const Vector<number2> &src) const
{
  Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
  AssertDimension(dst.size(), storage.n());
  AssertDimension(src.size(), storage.m());

  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;
  // m > n, m = src.n, n = dst.n

  // Multiply Q_n ... Q_2 Q_1 src
  // Where Q_i = I - v_i v_i^T
  for (size_type j = 0; j < n; ++j)
    {
      // sum = v_i^T dst
      number2 sum = diagonal[j] * (*aux)(j);
      for (size_type i = j + 1; i < m; ++i)
        sum += static_cast<number2>(storage(i, j)) * (*aux)(i);
      // dst -= v * sum
      (*aux)(j) -= sum * diagonal[j];
      for (size_type i = j + 1; i < m; ++i)
        (*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);
  AssertIsFinite(sum);

  // Compute solution
  storage.backward(dst, *aux);

  return std::sqrt(sum);
}



template <typename number>
template <typename number2>
double
Householder<number>::least_squares(BlockVector<number2> &      dst,
                                   const BlockVector<number2> &src) const
{
  Assert(!storage.empty(), typename FullMatrix<number2>::ExcEmptyMatrix());
  AssertDimension(dst.size(), storage.n());
  AssertDimension(src.size(), storage.m());

  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;
  // m > n, m = src.n, n = dst.n

  // Multiply Q_n ... Q_2 Q_1 src
  // Where Q_i = I-v_i v_i^T
  for (size_type j = 0; j < n; ++j)
    {
      // sum = v_i^T dst
      number2 sum = diagonal[j] * (*aux)(j);
      for (size_type i = j + 1; i < m; ++i)
        sum += storage(i, j) * (*aux)(i);
      // dst -= v * sum
      (*aux)(j) -= sum * diagonal[j];
      for (size_type i = j + 1; i < m; ++i)
        (*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);
  AssertIsFinite(sum);

  // backward works for
  // Vectors only, so copy
  // them before
  Vector<number2> v_dst, v_aux;
  v_dst = dst;
  v_aux = *aux;
  // Compute solution
  storage.backward(v_dst, v_aux);
  // copy the result back
  // to the BlockVector
  dst = v_dst;

  return std::sqrt(sum);
}


template <typename number>
template <class VectorType>
void
Householder<number>::vmult(VectorType &dst, const VectorType &src) const
{
  least_squares(dst, src);
}


template <typename number>
template <class VectorType>
void
Householder<number>::Tvmult(VectorType &, const VectorType &) const
{
  Assert(false, ExcNotImplemented());
}



#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif