SparseLU_Memory.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
/*
 
 * NOTE: This file is the modified version of [s,d,c,z]memory.c files in SuperLU
 
 * -- SuperLU routine (version 3.1) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * August 1, 2008
 *
 * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
 *
 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
 * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
 *
 * Permission is hereby granted to use or copy this program for any
 * purpose, provided the above notices are retained on all copies.
 * Permission to modify the code and to distribute modified code is
 * granted, provided the above notices are retained, and a notice that
 * the code was modified is included with the above copyright notice.
 */
 
#ifndef EIGEN_SPARSELU_MEMORY
#define EIGEN_SPARSELU_MEMORY
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
namespace internal {
 
enum { LUNoMarker = 3 };
enum { emptyIdxLU = -1 };
inline Index LUnumTempV(Index& m, Index& w, Index& t, Index& b) { return (std::max)(m, (t + b) * w); }
 
template <typename Scalar>
inline Index LUTempSpace(Index& m, Index& w) {
  return (2 * w + 4 + LUNoMarker) * m * sizeof(Index) + (w + 1) * m * sizeof(Scalar);
}

template <typename Scalar, typename StorageIndex>
template <typename VectorType>
Index SparseLUImpl<Scalar, StorageIndex>::expand(VectorType& vec, Index& length, Index nbElts, Index keep_prev,
                                                 Index& num_expansions) {
  float alpha = 1.5;  // Ratio of the memory increase
  Index new_len;      // New size of the allocated memory
 
  if (num_expansions == 0 || keep_prev)
    new_len = length;  // First time allocate requested
  else
    new_len = (std::max)(length + 1, Index(alpha * length));
 
  VectorType old_vec;  // Temporary vector to hold the previous values
  if (nbElts > 0) old_vec = vec.segment(0, nbElts);
 
    // Allocate or expand the current vector
#ifdef EIGEN_EXCEPTIONS
  try
#endif
  {
    vec.resize(new_len);
  }
#ifdef EIGEN_EXCEPTIONS
  catch (std::bad_alloc&)
#else
  if (!vec.size())
#endif
  {
    if (!num_expansions) {
      // First time to allocate from LUMemInit()
      // Let LUMemInit() deals with it.
      return -1;
    }
    if (keep_prev) {
      // In this case, the memory length should not not be reduced
      return new_len;
    } else {
      // Reduce the size and increase again
      Index tries = 0;  // Number of attempts
      do {
        alpha = (alpha + 1) / 2;
        new_len = (std::max)(length + 1, Index(alpha * length));
#ifdef EIGEN_EXCEPTIONS
        try
#endif
        {
          vec.resize(new_len);
        }
#ifdef EIGEN_EXCEPTIONS
        catch (std::bad_alloc&)
#else
        if (!vec.size())
#endif
        {
          tries += 1;
          if (tries > 10) return new_len;
        }
      } while (!vec.size());
    }
  }
  // Copy the previous values to the newly allocated space
  if (nbElts > 0) vec.segment(0, nbElts) = old_vec;
 
  length = new_len;
  if (num_expansions) ++num_expansions;
  return 0;
}

template <typename Scalar, typename StorageIndex>
Index SparseLUImpl<Scalar, StorageIndex>::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio,
                                                  Index panel_size, GlobalLU_t& glu) {
  Index& num_expansions = glu.num_expansions;  // No memory expansions so far
  num_expansions = 0;
  glu.nzumax = glu.nzlumax = (std::min)(fillratio * (annz + 1) / n, m) * n;  // estimated number of nonzeros in U
  glu.nzlmax = (std::max)(Index(4), fillratio) * (annz + 1) / 4;             // estimated  nnz in L factor
  // Return the estimated size to the user if necessary
  Index tempSpace;
  tempSpace = (2 * panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar);
  if (lwork == emptyIdxLU) {
    Index estimated_size;
    estimated_size = (5 * n + 5) * sizeof(Index) + tempSpace + (glu.nzlmax + glu.nzumax) * sizeof(Index) +
                     (glu.nzlumax + glu.nzumax) * sizeof(Scalar) + n;
    return estimated_size;
  }
 
  // Setup the required space
 
  // First allocate Integer pointers for L\U factors
  glu.xsup.resize(n + 1);
  glu.supno.resize(n + 1);
  glu.xlsub.resize(n + 1);
  glu.xlusup.resize(n + 1);
  glu.xusub.resize(n + 1);
 
  // Reserve memory for L/U factors
  do {
    if ((expand<ScalarVector>(glu.lusup, glu.nzlumax, 0, 0, num_expansions) < 0) ||
        (expand<ScalarVector>(glu.ucol, glu.nzumax, 0, 0, num_expansions) < 0) ||
        (expand<IndexVector>(glu.lsub, glu.nzlmax, 0, 0, num_expansions) < 0) ||
        (expand<IndexVector>(glu.usub, glu.nzumax, 0, 1, num_expansions) < 0)) {
      // Reduce the estimated size and retry
      glu.nzlumax /= 2;
      glu.nzumax /= 2;
      glu.nzlmax /= 2;
      if (glu.nzlumax < annz) return glu.nzlumax;
    }
  } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size());
 
  ++num_expansions;
  return 0;
 
}  // end LuMemInit

template <typename Scalar, typename StorageIndex>
template <typename VectorType>
Index SparseLUImpl<Scalar, StorageIndex>::memXpand(VectorType& vec, Index& maxlen, Index nbElts, MemType memtype,
                                                   Index& num_expansions) {
  Index failed_size;
  if (memtype == USUB)
    failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 1, num_expansions);
  else
    failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 0, num_expansions);
 
  if (failed_size) return failed_size;
 
  return 0;
}
 
}  // end namespace internal
 
}  // end namespace Eigen
#endif  // EIGEN_SPARSELU_MEMORY