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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_householder_h
#define dealii_householder_h


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

#include <vector>

DEAL_II_NAMESPACE_OPEN


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


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

  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