1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Dsir Nuentsa-Wakam <desire.nuentsa_wakam (at) inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 /* 11 12 * NOTE: This file is the modified version of [s,d,c,z]memory.c files in SuperLU 13 14 * -- SuperLU routine (version 3.1) -- 15 * Univ. of California Berkeley, Xerox Palo Alto Research Center, 16 * and Lawrence Berkeley National Lab. 17 * August 1, 2008 18 * 19 * Copyright (c) 1994 by Xerox Corporation. All rights reserved. 20 * 21 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 22 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 23 * 24 * Permission is hereby granted to use or copy this program for any 25 * purpose, provided the above notices are retained on all copies. 26 * Permission to modify the code and to distribute modified code is 27 * granted, provided the above notices are retained, and a notice that 28 * the code was modified is included with the above copyright notice. 29 */ 30 31 #ifndef EIGEN_SPARSELU_MEMORY 32 #define EIGEN_SPARSELU_MEMORY 33 34 namespace Eigen { 35 namespace internal { 36 37 enum { LUNoMarker = 3 }; 38 enum {emptyIdxLU = -1}; 39 inline Index LUnumTempV(Index& m, Index& w, Index& t, Index& b) 40 { 41 return (std::max)(m, (t+b)*w); 42 } 43 44 template< typename Scalar> 45 inline Index LUTempSpace(Index&m, Index& w) 46 { 47 return (2*w + 4 + LUNoMarker) * m * sizeof(Index) + (w + 1) * m * sizeof(Scalar); 48 } 49 50 51 52 53 /** 54 * Expand the existing storage to accomodate more fill-ins 55 * \param vec Valid pointer to the vector to allocate or expand 56 * \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 57 * \param[in] nbElts Current number of elements in the factors 58 * \param keep_prev 1: use length and do not expand the vector; 0: compute new_len and expand 59 * \param[in,out] num_expansions Number of times the memory has been expanded 60 */ 61 template <typename Scalar, typename StorageIndex> 62 template <typename VectorType> 63 Index SparseLUImpl<Scalar,StorageIndex>::expand(VectorType& vec, Index& length, Index nbElts, Index keep_prev, Index& num_expansions) 64 { 65 66 float alpha = 1.5; // Ratio of the memory increase 67 Index new_len; // New size of the allocated memory 68 69 if(num_expansions == 0 || keep_prev) 70 new_len = length ; // First time allocate requested 71 else 72 new_len = (std::max)(length+1,Index(alpha * length)); 73 74 VectorType old_vec; // Temporary vector to hold the previous values 75 if (nbElts > 0 ) 76 old_vec = vec.segment(0,nbElts); 77 78 //Allocate or expand the current vector 79 #ifdef EIGEN_EXCEPTIONS 80 try 81 #endif 82 { 83 vec.resize(new_len); 84 } 85 #ifdef EIGEN_EXCEPTIONS 86 catch(std::bad_alloc& ) 87 #else 88 if(!vec.size()) 89 #endif 90 { 91 if (!num_expansions) 92 { 93 // First time to allocate from LUMemInit() 94 // Let LUMemInit() deals with it. 95 return -1; 96 } 97 if (keep_prev) 98 { 99 // In this case, the memory length should not not be reduced 100 return new_len; 101 } 102 else 103 { 104 // Reduce the size and increase again 105 Index tries = 0; // Number of attempts 106 do 107 { 108 alpha = (alpha + 1)/2; 109 new_len = (std::max)(length+1,Index(alpha * length)); 110 #ifdef EIGEN_EXCEPTIONS 111 try 112 #endif 113 { 114 vec.resize(new_len); 115 } 116 #ifdef EIGEN_EXCEPTIONS 117 catch(std::bad_alloc& ) 118 #else 119 if (!vec.size()) 120 #endif 121 { 122 tries += 1; 123 if ( tries > 10) return new_len; 124 } 125 } while (!vec.size()); 126 } 127 } 128 //Copy the previous values to the newly allocated space 129 if (nbElts > 0) 130 vec.segment(0, nbElts) = old_vec; 131 132 133 length = new_len; 134 if(num_expansions) ++num_expansions; 135 return 0; 136 } 137 138 /** 139 * \brief Allocate various working space for the numerical factorization phase. 140 * \param m number of rows of the input matrix 141 * \param n number of columns 142 * \param annz number of initial nonzeros in the matrix 143 * \param lwork if lwork=-1, this routine returns an estimated size of the required memory 144 * \param glu persistent data to facilitate multiple factors : will be deleted later ?? 145 * \param fillratio estimated ratio of fill in the factors 146 * \param panel_size Size of a panel 147 * \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 148 * \note Unlike SuperLU, this routine does not support successive factorization with the same pattern and the same row permutation 149 */ 150 template <typename Scalar, typename StorageIndex> 151 Index SparseLUImpl<Scalar,StorageIndex>::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size, GlobalLU_t& glu) 152 { 153 Index& num_expansions = glu.num_expansions; //No memory expansions so far 154 num_expansions = 0; 155 glu.nzumax = glu.nzlumax = (std::min)(fillratio * (annz+1) / n, m) * n; // estimated number of nonzeros in U 156 glu.nzlmax = (std::max)(Index(4), fillratio) * (annz+1) / 4; // estimated nnz in L factor 157 // Return the estimated size to the user if necessary 158 Index tempSpace; 159 tempSpace = (2*panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar); 160 if (lwork == emptyIdxLU) 161 { 162 Index estimated_size; 163 estimated_size = (5 * n + 5) * sizeof(Index) + tempSpace 164 + (glu.nzlmax + glu.nzumax) * sizeof(Index) + (glu.nzlumax+glu.nzumax) * sizeof(Scalar) + n; 165 return estimated_size; 166 } 167 168 // Setup the required space 169 170 // First allocate Integer pointers for L\U factors 171 glu.xsup.resize(n+1); 172 glu.supno.resize(n+1); 173 glu.xlsub.resize(n+1); 174 glu.xlusup.resize(n+1); 175 glu.xusub.resize(n+1); 176 177 // Reserve memory for L/U factors 178 do 179 { 180 if( (expand<ScalarVector>(glu.lusup, glu.nzlumax, 0, 0, num_expansions)<0) 181 || (expand<ScalarVector>(glu.ucol, glu.nzumax, 0, 0, num_expansions)<0) 182 || (expand<IndexVector> (glu.lsub, glu.nzlmax, 0, 0, num_expansions)<0) 183 || (expand<IndexVector> (glu.usub, glu.nzumax, 0, 1, num_expansions)<0) ) 184 { 185 //Reduce the estimated size and retry 186 glu.nzlumax /= 2; 187 glu.nzumax /= 2; 188 glu.nzlmax /= 2; 189 if (glu.nzlumax < annz ) return glu.nzlumax; 190 } 191 } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size()); 192 193 ++num_expansions; 194 return 0; 195 196 } // end LuMemInit 197 198 /** 199 * \brief Expand the existing storage 200 * \param vec vector to expand 201 * \param[in,out] maxlen On input, previous size of vec (Number of elements to copy ). on output, new size 202 * \param nbElts current number of elements in the vector. 203 * \param memtype Type of the element to expand 204 * \param num_expansions Number of expansions 205 * \return 0 on success, > 0 size of the memory allocated so far 206 */ 207 template <typename Scalar, typename StorageIndex> 208 template <typename VectorType> 209 Index SparseLUImpl<Scalar,StorageIndex>::memXpand(VectorType& vec, Index& maxlen, Index nbElts, MemType memtype, Index& num_expansions) 210 { 211 Index failed_size; 212 if (memtype == USUB) 213 failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 1, num_expansions); 214 else 215 failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 0, num_expansions); 216 217 if (failed_size) 218 return failed_size; 219 220 return 0 ; 221 } 222 223 } // end namespace internal 224 225 } // end namespace Eigen 226 #endif // EIGEN_SPARSELU_MEMORY 227
