sparse_direct.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2019 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.
//
// ---------------------------------------------------------------------

#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/thread_management.h>

#include <deal.II/lac/block_sparse_matrix.h>
#include <deal.II/lac/sparse_direct.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>

#include <cerrno>
#include <iostream>
#include <list>
#include <typeinfo>
#include <vector>


DEAL_II_NAMESPACE_OPEN


// include UMFPACK file.
#ifdef DEAL_II_WITH_UMFPACK
#  include <umfpack.h>
#endif



SparseDirectUMFPACK::~SparseDirectUMFPACK()
{
  clear();
}


void
SparseDirectUMFPACK::initialize(const SparsityPattern &)
{}


#ifdef DEAL_II_WITH_UMFPACK

SparseDirectUMFPACK::SparseDirectUMFPACK()
  : _m(0)
  , _n(0)
  , symbolic_decomposition(nullptr)
  , numeric_decomposition(nullptr)
  , control(UMFPACK_CONTROL)
{
  umfpack_dl_defaults(control.data());
}



void
SparseDirectUMFPACK::clear()
{
  // delete objects that haven't been deleted yet
  if (symbolic_decomposition != nullptr)
    {
      umfpack_dl_free_symbolic(&symbolic_decomposition);
      symbolic_decomposition = nullptr;
    }

  if (numeric_decomposition != nullptr)
    {
      umfpack_dl_free_numeric(&numeric_decomposition);
      numeric_decomposition = nullptr;
    }

  {
    std::vector<long int> tmp;
    tmp.swap(Ap);
  }

  {
    std::vector<long int> tmp;
    tmp.swap(Ai);
  }

  {
    std::vector<double> tmp;
    tmp.swap(Ax);
  }

  umfpack_dl_defaults(control.data());
}



template <typename number>
void
SparseDirectUMFPACK::sort_arrays(const SparseMatrix<number> &matrix)
{
  // do the copying around of entries so that the diagonal entry is in the
  // right place. note that this is easy to detect: since all entries apart
  // from the diagonal entry are sorted, we know that the diagonal entry is
  // in the wrong place if and only if its column index is larger than the
  // column index of the second entry in a row
  //
  // ignore rows with only one or no entry
  for (size_type row = 0; row < matrix.m(); ++row)
    {
      // we may have to move some elements that are left of the diagonal
      // but presently after the diagonal entry to the left, whereas the
      // diagonal entry has to move to the right. we could first figure out
      // where to move everything to, but for simplicity we just make a
      // series of swaps instead (this is kind of a single run of
      // bubble-sort, which gives us the desired result since the array is
      // already "almost" sorted)
      //
      // in the first loop, the condition in the while-header also checks
      // that the row has at least two entries and that the diagonal entry
      // is really in the wrong place
      long int cursor = Ap[row];
      while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
        {
          std::swap(Ai[cursor], Ai[cursor + 1]);
          std::swap(Ax[cursor], Ax[cursor + 1]);
          ++cursor;
        }
    }
}



template <typename number>
void
SparseDirectUMFPACK::sort_arrays(const SparseMatrixEZ<number> &matrix)
{
  // same thing for SparseMatrixEZ
  for (size_type row = 0; row < matrix.m(); ++row)
    {
      long int cursor = Ap[row];
      while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] > Ai[cursor + 1]))
        {
          std::swap(Ai[cursor], Ai[cursor + 1]);
          std::swap(Ax[cursor], Ax[cursor + 1]);
          ++cursor;
        }
    }
}



template <typename number>
void
SparseDirectUMFPACK::sort_arrays(const BlockSparseMatrix<number> &matrix)
{
  // the case for block matrices is a bit more difficult, since all we know
  // is that *within each block*, the diagonal of that block may come
  // first. however, that means that there may be as many entries per row
  // in the wrong place as there are block columns. we can do the same
  // thing as above, but we have to do it multiple times
  for (size_type row = 0; row < matrix.m(); ++row)
    {
      long int cursor = Ap[row];
      for (size_type block = 0; block < matrix.n_block_cols(); ++block)
        {
          // find the next out-of-order element
          while ((cursor < Ap[row + 1] - 1) && (Ai[cursor] < Ai[cursor + 1]))
            ++cursor;

          // if there is none, then just go on
          if (cursor == Ap[row + 1] - 1)
            break;

          // otherwise swap this entry with successive ones as long as
          // necessary
          long int element = cursor;
          while ((element < Ap[row + 1] - 1) && (Ai[element] > Ai[element + 1]))
            {
              std::swap(Ai[element], Ai[element + 1]);
              std::swap(Ax[element], Ax[element + 1]);
              ++element;
            }
        }
    }
}



template <class Matrix>
void
SparseDirectUMFPACK::factorize(const Matrix &matrix)
{
  Assert(matrix.m() == matrix.n(), ExcNotQuadratic())

    clear();

  _m = matrix.m();
  _n = matrix.n();

  const size_type N = matrix.m();

  // copy over the data from the matrix to the data structures UMFPACK
  // wants. note two things: first, UMFPACK wants compressed column storage
  // whereas we always do compressed row storage; we work around this by,
  // rather than shuffling things around, copy over the data we have, but
  // then call the umfpack_dl_solve function with the UMFPACK_At argument,
  // meaning that we want to solve for the transpose system
  //
  // second: the data we have in the sparse matrices is "almost" right
  // already; UMFPACK wants the entries in each row (i.e. really: column)
  // to be sorted in ascending order. we almost have that, except that we
  // usually store the diagonal first in each row to allow for some
  // optimizations. thus, we have to resort things a little bit, but only
  // within each row
  //
  // final note: if the matrix has entries in the sparsity pattern that are
  // actually occupied by entries that have a zero numerical value, then we
  // keep them anyway. people are supposed to provide accurate sparsity
  // patterns.
  Ap.resize(N + 1);
  Ai.resize(matrix.n_nonzero_elements());
  Ax.resize(matrix.n_nonzero_elements());

  // first fill row lengths array
  Ap[0] = 0;
  for (size_type row = 1; row <= N; ++row)
    Ap[row] = Ap[row - 1] + matrix.get_row_length(row - 1);
  Assert(static_cast<size_type>(Ap.back()) == Ai.size(), ExcInternalError());

  // then copy over matrix elements. note that for sparse matrices,
  // iterators are sorted so that they traverse each row from start to end
  // before moving on to the next row. however, this isn't true for block
  // matrices, so we have to do a bit of book keeping
  {
    // have an array that for each row points to the first entry not yet
    // written to
    std::vector<long int> row_pointers = Ap;

    // loop over the elements of the matrix row by row, as suggested in the
    // documentation of the sparse matrix iterator class
    for (size_type row = 0; row < matrix.m(); ++row)
      {
        for (typename Matrix::const_iterator p = matrix.begin(row);
             p != matrix.end(row);
             ++p)
          {
            // write entry into the first free one for this row
            Ai[row_pointers[row]] = p->column();
            Ax[row_pointers[row]] = p->value();

            // then move pointer ahead
            ++row_pointers[row];
          }
      }

    // at the end, we should have written all rows completely
    for (size_type i = 0; i < Ap.size() - 1; ++i)
      Assert(row_pointers[i] == Ap[i + 1], ExcInternalError());
  }

  // make sure that the elements in each row are sorted. we have to be more
  // careful for block sparse matrices, so ship this task out to a
  // different function
  sort_arrays(matrix);

  int status;
  status = umfpack_dl_symbolic(N,
                               N,
                               Ap.data(),
                               Ai.data(),
                               Ax.data(),
                               &symbolic_decomposition,
                               control.data(),
                               nullptr);
  AssertThrow(status == UMFPACK_OK,
              ExcUMFPACKError("umfpack_dl_symbolic", status));

  status = umfpack_dl_numeric(Ap.data(),
                              Ai.data(),
                              Ax.data(),
                              symbolic_decomposition,
                              &numeric_decomposition,
                              control.data(),
                              nullptr);
  AssertThrow(status == UMFPACK_OK,
              ExcUMFPACKError("umfpack_dl_numeric", status));

  umfpack_dl_free_symbolic(&symbolic_decomposition);
}



void
SparseDirectUMFPACK::solve(Vector<double> &rhs_and_solution,
                           bool            transpose /*=false*/) const
{
  // make sure that some kind of factorize() call has happened before
  Assert(Ap.size() != 0, ExcNotInitialized());
  Assert(Ai.size() != 0, ExcNotInitialized());
  Assert(Ai.size() == Ax.size(), ExcNotInitialized());

  Vector<double> rhs(rhs_and_solution.size());
  rhs = rhs_and_solution;

  // solve the system. note that since UMFPACK wants compressed column
  // storage instead of the compressed row storage format we use in
  // deal.II's SparsityPattern classes, we solve for UMFPACK's A^T instead

  // Conversely, if we solve for the transpose, we have to use UMFPACK_A
  // instead.
  const int status = umfpack_dl_solve(transpose ? UMFPACK_A : UMFPACK_At,
                                      Ap.data(),
                                      Ai.data(),
                                      Ax.data(),
                                      rhs_and_solution.begin(),
                                      rhs.begin(),
                                      numeric_decomposition,
                                      control.data(),
                                      nullptr);
  AssertThrow(status == UMFPACK_OK,
              ExcUMFPACKError("umfpack_dl_solve", status));
}


void
SparseDirectUMFPACK::solve(BlockVector<double> &rhs_and_solution,
                           bool                 transpose /*=false*/) const
{
  // the UMFPACK functions want a contiguous array of elements, so
  // there is no way around copying data around. thus, just copy the
  // data into a regular vector and back
  Vector<double> tmp(rhs_and_solution.size());
  tmp = rhs_and_solution;
  solve(tmp, transpose);
  rhs_and_solution = tmp;
}



template <class Matrix>
void
SparseDirectUMFPACK::solve(const Matrix &  matrix,
                           Vector<double> &rhs_and_solution,
                           bool            transpose /*=false*/)
{
  factorize(matrix);
  solve(rhs_and_solution, transpose);
}


template <class Matrix>
void
SparseDirectUMFPACK::solve(const Matrix &       matrix,
                           BlockVector<double> &rhs_and_solution,
                           bool                 transpose /*=false*/)
{
  factorize(matrix);
  solve(rhs_and_solution, transpose);
}


#else


SparseDirectUMFPACK::SparseDirectUMFPACK()
  : _m(0)
  , _n(0)
  , symbolic_decomposition(nullptr)
  , numeric_decomposition(nullptr)
  , control(0)
{}


void
SparseDirectUMFPACK::clear()
{}


template <class Matrix>
void
SparseDirectUMFPACK::factorize(const Matrix &)
{
  AssertThrow(
    false,
    ExcMessage(
      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
}


void
SparseDirectUMFPACK::solve(Vector<double> &, bool) const
{
  AssertThrow(
    false,
    ExcMessage(
      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
}



void
SparseDirectUMFPACK::solve(BlockVector<double> &, bool) const
{
  AssertThrow(
    false,
    ExcMessage(
      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
}


template <class Matrix>
void
SparseDirectUMFPACK::solve(const Matrix &, Vector<double> &, bool)
{
  AssertThrow(
    false,
    ExcMessage(
      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
}



template <class Matrix>
void
SparseDirectUMFPACK::solve(const Matrix &, BlockVector<double> &, bool)
{
  AssertThrow(
    false,
    ExcMessage(
      "To call this function you need UMFPACK, but you configured deal.II without passing the necessary switch to 'cmake'. Please consult the installation instructions in doc/readme.html."));
}

#endif


template <class Matrix>
void
SparseDirectUMFPACK::initialize(const Matrix &M, const AdditionalData)
{
  this->factorize(M);
}


void
SparseDirectUMFPACK::vmult(Vector<double> &dst, const Vector<double> &src) const
{
  dst = src;
  this->solve(dst);
}



void
SparseDirectUMFPACK::vmult(BlockVector<double> &      dst,
                           const BlockVector<double> &src) const
{
  dst = src;
  this->solve(dst);
}


void
SparseDirectUMFPACK::Tvmult(Vector<double> &      dst,
                            const Vector<double> &src) const
{
  dst = src;
  this->solve(dst, /*transpose=*/true);
}



void
SparseDirectUMFPACK::Tvmult(BlockVector<double> &      dst,
                            const BlockVector<double> &src) const
{
  dst = src;
  this->solve(dst, /*transpose=*/true);
}

SparseDirectUMFPACK::size_type
SparseDirectUMFPACK::m() const
{
  Assert(_m != 0, ExcNotInitialized());
  return _m;
}

SparseDirectUMFPACK::size_type
SparseDirectUMFPACK::n() const
{
  Assert(_n != 0, ExcNotInitialized());
  return _n;
}


// explicit instantiations for SparseMatrixUMFPACK
#define InstantiateUMFPACK(MatrixType)                              \
  template void SparseDirectUMFPACK::factorize(const MatrixType &); \
  template void SparseDirectUMFPACK::solve(const MatrixType &,      \
                                           Vector<double> &,        \
                                           bool);                   \
  template void SparseDirectUMFPACK::solve(const MatrixType &,      \
                                           BlockVector<double> &,   \
                                           bool);                   \
  template void SparseDirectUMFPACK::initialize(const MatrixType &, \
                                                const AdditionalData)

InstantiateUMFPACK(SparseMatrix<double>);
InstantiateUMFPACK(SparseMatrix<float>);
InstantiateUMFPACK(SparseMatrixEZ<double>);
InstantiateUMFPACK(SparseMatrixEZ<float>);
InstantiateUMFPACK(BlockSparseMatrix<double>);
InstantiateUMFPACK(BlockSparseMatrix<float>);


DEAL_II_NAMESPACE_CLOSE