![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 /* 00011 00012 * NOTE: This file is the modified version of [s,d,c,z]memory.c files in SuperLU 00013 00014 * -- SuperLU routine (version 3.1) -- 00015 * Univ. of California Berkeley, Xerox Palo Alto Research Center, 00016 * and Lawrence Berkeley National Lab. 00017 * August 1, 2008 00018 * 00019 * Copyright (c) 1994 by Xerox Corporation. All rights reserved. 00020 * 00021 * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY 00022 * EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK. 00023 * 00024 * Permission is hereby granted to use or copy this program for any 00025 * purpose, provided the above notices are retained on all copies. 00026 * Permission to modify the code and to distribute modified code is 00027 * granted, provided the above notices are retained, and a notice that 00028 * the code was modified is included with the above copyright notice. 00029 */ 00030 00031 #ifndef EIGEN_SPARSELU_MEMORY 00032 #define EIGEN_SPARSELU_MEMORY 00033 00034 namespace Eigen { 00035 namespace internal { 00036 00037 enum { LUNoMarker = 3 }; 00038 enum {emptyIdxLU = -1}; 00039 inline Index LUnumTempV(Index& m, Index& w, Index& t, Index& b) 00040 { 00041 return (std::max)(m, (t+b)*w); 00042 } 00043 00044 template< typename Scalar> 00045 inline Index LUTempSpace(Index&m, Index& w) 00046 { 00047 return (2*w + 4 + LUNoMarker) * m * sizeof(Index) + (w + 1) * m * sizeof(Scalar); 00048 } 00049 00050 00051 00052 00061 template <typename Scalar, typename StorageIndex> 00062 template <typename VectorType> 00063 Index SparseLUImpl<Scalar,StorageIndex>::expand(VectorType& vec, Index& length, Index nbElts, Index keep_prev, Index& num_expansions) 00064 { 00065 00066 float alpha = 1.5; // Ratio of the memory increase 00067 Index new_len; // New size of the allocated memory 00068 00069 if(num_expansions == 0 || keep_prev) 00070 new_len = length ; // First time allocate requested 00071 else 00072 new_len = (std::max)(length+1,Index(alpha * length)); 00073 00074 VectorType old_vec; // Temporary vector to hold the previous values 00075 if (nbElts > 0 ) 00076 old_vec = vec.segment(0,nbElts); 00077 00078 //Allocate or expand the current vector 00079 #ifdef EIGEN_EXCEPTIONS 00080 try 00081 #endif 00082 { 00083 vec.resize(new_len); 00084 } 00085 #ifdef EIGEN_EXCEPTIONS 00086 catch(std::bad_alloc& ) 00087 #else 00088 if(!vec.size()) 00089 #endif 00090 { 00091 if (!num_expansions) 00092 { 00093 // First time to allocate from LUMemInit() 00094 // Let LUMemInit() deals with it. 00095 return -1; 00096 } 00097 if (keep_prev) 00098 { 00099 // In this case, the memory length should not not be reduced 00100 return new_len; 00101 } 00102 else 00103 { 00104 // Reduce the size and increase again 00105 Index tries = 0; // Number of attempts 00106 do 00107 { 00108 alpha = (alpha + 1)/2; 00109 new_len = (std::max)(length+1,Index(alpha * length)); 00110 #ifdef EIGEN_EXCEPTIONS 00111 try 00112 #endif 00113 { 00114 vec.resize(new_len); 00115 } 00116 #ifdef EIGEN_EXCEPTIONS 00117 catch(std::bad_alloc& ) 00118 #else 00119 if (!vec.size()) 00120 #endif 00121 { 00122 tries += 1; 00123 if ( tries > 10) return new_len; 00124 } 00125 } while (!vec.size()); 00126 } 00127 } 00128 //Copy the previous values to the newly allocated space 00129 if (nbElts > 0) 00130 vec.segment(0, nbElts) = old_vec; 00131 00132 00133 length = new_len; 00134 if(num_expansions) ++num_expansions; 00135 return 0; 00136 } 00137 00150 template <typename Scalar, typename StorageIndex> 00151 Index SparseLUImpl<Scalar,StorageIndex>::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size, GlobalLU_t& glu) 00152 { 00153 Index& num_expansions = glu.num_expansions; //No memory expansions so far 00154 num_expansions = 0; 00155 glu.nzumax = glu.nzlumax = (std::min)(fillratio * (annz+1) / n, m) * n; // estimated number of nonzeros in U 00156 glu.nzlmax = (std::max)(Index(4), fillratio) * (annz+1) / 4; // estimated nnz in L factor 00157 // Return the estimated size to the user if necessary 00158 Index tempSpace; 00159 tempSpace = (2*panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar); 00160 if (lwork == emptyIdxLU) 00161 { 00162 Index estimated_size; 00163 estimated_size = (5 * n + 5) * sizeof(Index) + tempSpace 00164 + (glu.nzlmax + glu.nzumax) * sizeof(Index) + (glu.nzlumax+glu.nzumax) * sizeof(Scalar) + n; 00165 return estimated_size; 00166 } 00167 00168 // Setup the required space 00169 00170 // First allocate Integer pointers for L\U factors 00171 glu.xsup.resize(n+1); 00172 glu.supno.resize(n+1); 00173 glu.xlsub.resize(n+1); 00174 glu.xlusup.resize(n+1); 00175 glu.xusub.resize(n+1); 00176 00177 // Reserve memory for L/U factors 00178 do 00179 { 00180 if( (expand<ScalarVector>(glu.lusup, glu.nzlumax, 0, 0, num_expansions)<0) 00181 || (expand<ScalarVector>(glu.ucol, glu.nzumax, 0, 0, num_expansions)<0) 00182 || (expand<IndexVector> (glu.lsub, glu.nzlmax, 0, 0, num_expansions)<0) 00183 || (expand<IndexVector> (glu.usub, glu.nzumax, 0, 1, num_expansions)<0) ) 00184 { 00185 //Reduce the estimated size and retry 00186 glu.nzlumax /= 2; 00187 glu.nzumax /= 2; 00188 glu.nzlmax /= 2; 00189 if (glu.nzlumax < annz ) return glu.nzlumax; 00190 } 00191 } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size()); 00192 00193 ++num_expansions; 00194 return 0; 00195 00196 } // end LuMemInit 00197 00207 template <typename Scalar, typename StorageIndex> 00208 template <typename VectorType> 00209 Index SparseLUImpl<Scalar,StorageIndex>::memXpand(VectorType& vec, Index& maxlen, Index nbElts, MemType memtype, Index& num_expansions) 00210 { 00211 Index failed_size; 00212 if (memtype == USUB) 00213 failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 1, num_expansions); 00214 else 00215 failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 0, num_expansions); 00216 00217 if (failed_size) 00218 return failed_size; 00219 00220 return 0 ; 00221 } 00222 00223 } // end namespace internal 00224 00225 } // end namespace Eigen 00226 #endif // EIGEN_SPARSELU_MEMORY