src/SparseLU/SparseLU_Memory.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 /*
0011
0012  * NOTE: This file is the modified version of [s,d,c,z]memory.c files in SuperLU
0013
0014  * -- SuperLU routine (version 3.1) --
0015  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
0016  * and Lawrence Berkeley National Lab.
0017  * August 1, 2008
0018  *
0019  * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
0020  *
0021  * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
0022  * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
0023  *
0024  * Permission is hereby granted to use or copy this program for any
0025  * purpose, provided the above notices are retained on all copies.
0026  * Permission to modify the code and to distribute modified code is
0027  * granted, provided the above notices are retained, and a notice that
0028  * the code was modified is included with the above copyright notice.
0029  */
0030
0031 #ifndef EIGEN_SPARSELU_MEMORY
0032 #define EIGEN_SPARSELU_MEMORY
0033
0034 namespace RivetEigen {
0035 namespace internal {
0036
0037 enum { LUNoMarker = 3 };
0038 enum {emptyIdxLU = -1};
0039 inline Index LUnumTempV(Index& m, Index& w, Index& t, Index& b)
0040 {
0041   return (std::max)(m, (t+b)*w);
0042 }
0043
0044 template< typename Scalar>
0045 inline Index LUTempSpace(Index&m, Index& w)
0046 {
0047   return (2*w + 4 + LUNoMarker) * m * sizeof(Index) + (w + 1) * m * sizeof(Scalar);
0048 }
0049
0050
0051
0052
0053 /**
0054   * Expand the existing storage to accommodate more fill-ins
0055   * \param vec Valid pointer to the vector to allocate or expand
0056   * \param[in,out] length  At input, contain the current length of the vector that is to be increased. At output, length of the newly allocated vector
0057   * \param[in] nbElts Current number of elements in the factors
0058   * \param keep_prev  1: use length  and do not expand the vector; 0: compute new_len and expand
0059   * \param[in,out] num_expansions Number of times the memory has been expanded
0060   */
0061 template <typename Scalar, typename StorageIndex>
0062 template <typename VectorType>
0063 Index  SparseLUImpl<Scalar,StorageIndex>::expand(VectorType& vec, Index& length, Index nbElts, Index keep_prev, Index& num_expansions)
0064 {
0065
0066   float alpha = 1.5; // Ratio of the memory increase
0067   Index new_len; // New size of the allocated memory
0068
0069   if(num_expansions == 0 || keep_prev)
0070     new_len = length ; // First time allocate requested
0071   else
0072     new_len = (std::max)(length+1,Index(alpha * length));
0073
0074   VectorType old_vec; // Temporary vector to hold the previous values
0075   if (nbElts > 0 )
0076     old_vec = vec.segment(0,nbElts);
0077
0078   //Allocate or expand the current vector
0079 #ifdef EIGEN_EXCEPTIONS
0080   try
0081 #endif
0082   {
0083     vec.resize(new_len);
0084   }
0085 #ifdef EIGEN_EXCEPTIONS
0086   catch(std::bad_alloc& )
0087 #else
0088   if(!vec.size())
0089 #endif
0090   {
0091     if (!num_expansions)
0092     {
0093       // First time to allocate from LUMemInit()
0094       // Let LUMemInit() deals with it.
0095       return -1;
0096     }
0097     if (keep_prev)
0098     {
0099       // In this case, the memory length should not not be reduced
0100       return new_len;
0101     }
0102     else
0103     {
0104       // Reduce the size and increase again
0105       Index tries = 0; // Number of attempts
0106       do
0107       {
0108         alpha = (alpha + 1)/2;
0109         new_len = (std::max)(length+1,Index(alpha * length));
0110 #ifdef EIGEN_EXCEPTIONS
0111         try
0112 #endif
0113         {
0114           vec.resize(new_len);
0115         }
0116 #ifdef EIGEN_EXCEPTIONS
0117         catch(std::bad_alloc& )
0118 #else
0119         if (!vec.size())
0120 #endif
0121         {
0122           tries += 1;
0123           if ( tries > 10) return new_len;
0124         }
0125       } while (!vec.size());
0126     }
0127   }
0128   //Copy the previous values to the newly allocated space
0129   if (nbElts > 0)
0130     vec.segment(0, nbElts) = old_vec;
0131
0132
0133   length  = new_len;
0134   if(num_expansions) ++num_expansions;
0135   return 0;
0136 }
0137
0138 /**
0139  * \brief  Allocate various working space for the numerical factorization phase.
0140  * \param m number of rows of the input matrix
0141  * \param n number of columns
0142  * \param annz number of initial nonzeros in the matrix
0143  * \param lwork  if lwork=-1, this routine returns an estimated size of the required memory
0144  * \param glu persistent data to facilitate multiple factors : will be deleted later ??
0145  * \param fillratio estimated ratio of fill in the factors
0146  * \param panel_size Size of a panel
0147  * \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated memory when allocation failed, and 0 on success
0148  * \note Unlike SuperLU, this routine does not support successive factorization with the same pattern and the same row permutation
0149  */
0150 template <typename Scalar, typename StorageIndex>
0151 Index SparseLUImpl<Scalar,StorageIndex>::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size,  GlobalLU_t& glu)
0152 {
0153   Index& num_expansions = glu.num_expansions; //No memory expansions so far
0154   num_expansions = 0;
0155   glu.nzumax = glu.nzlumax = (std::min)(fillratio * (annz+1) / n, m) * n; // estimated number of nonzeros in U
0156   glu.nzlmax = (std::max)(Index(4), fillratio) * (annz+1) / 4; // estimated  nnz in L factor
0157   // Return the estimated size to the user if necessary
0158   Index tempSpace;
0159   tempSpace = (2*panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar);
0160   if (lwork == emptyIdxLU)
0161   {
0162     Index estimated_size;
0163     estimated_size = (5 * n + 5) * sizeof(Index)  + tempSpace
0164                     + (glu.nzlmax + glu.nzumax) * sizeof(Index) + (glu.nzlumax+glu.nzumax) *  sizeof(Scalar) + n;
0165     return estimated_size;
0166   }
0167
0168   // Setup the required space
0169
0170   // First allocate Integer pointers for L\U factors
0171   glu.xsup.resize(n+1);
0172   glu.supno.resize(n+1);
0173   glu.xlsub.resize(n+1);
0174   glu.xlusup.resize(n+1);
0175   glu.xusub.resize(n+1);
0176
0177   // Reserve memory for L/U factors
0178   do
0179   {
0180     if(     (expand<ScalarVector>(glu.lusup, glu.nzlumax, 0, 0, num_expansions)<0)
0181         ||  (expand<ScalarVector>(glu.ucol,  glu.nzumax,  0, 0, num_expansions)<0)
0182         ||  (expand<IndexVector> (glu.lsub,  glu.nzlmax,  0, 0, num_expansions)<0)
0183         ||  (expand<IndexVector> (glu.usub,  glu.nzumax,  0, 1, num_expansions)<0) )
0184     {
0185       //Reduce the estimated size and retry
0186       glu.nzlumax /= 2;
0187       glu.nzumax /= 2;
0188       glu.nzlmax /= 2;
0189       if (glu.nzlumax < annz ) return glu.nzlumax;
0190     }
0191   } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size());
0192
0193   ++num_expansions;
0194   return 0;
0195
0196 } // end LuMemInit
0197
0198 /**
0199  * \brief Expand the existing storage
0200  * \param vec vector to expand
0201  * \param[in,out] maxlen On input, previous size of vec (Number of elements to copy ). on output, new size
0202  * \param nbElts current number of elements in the vector.
0203  * \param memtype Type of the element to expand
0204  * \param num_expansions Number of expansions
0205  * \return 0 on success, > 0 size of the memory allocated so far
0206  */
0207 template <typename Scalar, typename StorageIndex>
0208 template <typename VectorType>
0209 Index SparseLUImpl<Scalar,StorageIndex>::memXpand(VectorType& vec, Index& maxlen, Index nbElts, MemType memtype, Index& num_expansions)
0210 {
0211   Index failed_size;
0212   if (memtype == USUB)
0213      failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 1, num_expansions);
0214   else
0215     failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 0, num_expansions);
0216
0217   if (failed_size)
0218     return failed_size;
0219
0220   return 0 ;
0221 }
0222
0223 } // end namespace internal
0224
0225 } // end namespace RivetEigen
0226 #endif // EIGEN_SPARSELU_MEMORY