$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 template<typename Index> 00040 inline Index LUnumTempV(Index& m, Index& w, Index& t, Index& b) 00041 { 00042 return (std::max)(m, (t+b)*w); 00043 } 00044 00045 template< typename Scalar, typename Index> 00046 inline Index LUTempSpace(Index&m, Index& w) 00047 { 00048 return (2*w + 4 + LUNoMarker) * m * sizeof(Index) + (w + 1) * m * sizeof(Scalar); 00049 } 00050 00051 00052 00053 00062 template <typename Scalar, typename Index> 00063 template <typename VectorType> 00064 Index SparseLUImpl<Scalar,Index>::expand(VectorType& vec, Index& length, Index nbElts, Index keep_prev, Index& num_expansions) 00065 { 00066 00067 float alpha = 1.5; // Ratio of the memory increase 00068 Index new_len; // New size of the allocated memory 00069 00070 if(num_expansions == 0 || keep_prev) 00071 new_len = length ; // First time allocate requested 00072 else 00073 new_len = (std::max)(length+1,Index(alpha * length)); 00074 00075 VectorType old_vec; // Temporary vector to hold the previous values 00076 if (nbElts > 0 ) 00077 old_vec = vec.segment(0,nbElts); 00078 00079 //Allocate or expand the current vector 00080 #ifdef EIGEN_EXCEPTIONS 00081 try 00082 #endif 00083 { 00084 vec.resize(new_len); 00085 } 00086 #ifdef EIGEN_EXCEPTIONS 00087 catch(std::bad_alloc& ) 00088 #else 00089 if(!vec.size()) 00090 #endif 00091 { 00092 if (!num_expansions) 00093 { 00094 // First time to allocate from LUMemInit() 00095 // Let LUMemInit() deals with it. 00096 return -1; 00097 } 00098 if (keep_prev) 00099 { 00100 // In this case, the memory length should not not be reduced 00101 return new_len; 00102 } 00103 else 00104 { 00105 // Reduce the size and increase again 00106 Index tries = 0; // Number of attempts 00107 do 00108 { 00109 alpha = (alpha + 1)/2; 00110 new_len = (std::max)(length+1,Index(alpha * length)); 00111 #ifdef EIGEN_EXCEPTIONS 00112 try 00113 #endif 00114 { 00115 vec.resize(new_len); 00116 } 00117 #ifdef EIGEN_EXCEPTIONS 00118 catch(std::bad_alloc& ) 00119 #else 00120 if (!vec.size()) 00121 #endif 00122 { 00123 tries += 1; 00124 if ( tries > 10) return new_len; 00125 } 00126 } while (!vec.size()); 00127 } 00128 } 00129 //Copy the previous values to the newly allocated space 00130 if (nbElts > 0) 00131 vec.segment(0, nbElts) = old_vec; 00132 00133 00134 length = new_len; 00135 if(num_expansions) ++num_expansions; 00136 return 0; 00137 } 00138 00151 template <typename Scalar, typename Index> 00152 Index SparseLUImpl<Scalar,Index>::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size, GlobalLU_t& glu) 00153 { 00154 Index& num_expansions = glu.num_expansions; //No memory expansions so far 00155 num_expansions = 0; 00156 glu.nzumax = glu.nzlumax = (std::min)(fillratio * annz / n, m) * n; // estimated number of nonzeros in U 00157 glu.nzlmax = (std::max)(Index(4), fillratio) * annz / 4; // estimated nnz in L factor 00158 // Return the estimated size to the user if necessary 00159 Index tempSpace; 00160 tempSpace = (2*panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar); 00161 if (lwork == emptyIdxLU) 00162 { 00163 Index estimated_size; 00164 estimated_size = (5 * n + 5) * sizeof(Index) + tempSpace 00165 + (glu.nzlmax + glu.nzumax) * sizeof(Index) + (glu.nzlumax+glu.nzumax) * sizeof(Scalar) + n; 00166 return estimated_size; 00167 } 00168 00169 // Setup the required space 00170 00171 // First allocate Integer pointers for L\U factors 00172 glu.xsup.resize(n+1); 00173 glu.supno.resize(n+1); 00174 glu.xlsub.resize(n+1); 00175 glu.xlusup.resize(n+1); 00176 glu.xusub.resize(n+1); 00177 00178 // Reserve memory for L/U factors 00179 do 00180 { 00181 if( (expand<ScalarVector>(glu.lusup, glu.nzlumax, 0, 0, num_expansions)<0) 00182 || (expand<ScalarVector>(glu.ucol, glu.nzumax, 0, 0, num_expansions)<0) 00183 || (expand<IndexVector> (glu.lsub, glu.nzlmax, 0, 0, num_expansions)<0) 00184 || (expand<IndexVector> (glu.usub, glu.nzumax, 0, 1, num_expansions)<0) ) 00185 { 00186 //Reduce the estimated size and retry 00187 glu.nzlumax /= 2; 00188 glu.nzumax /= 2; 00189 glu.nzlmax /= 2; 00190 if (glu.nzlumax < annz ) return glu.nzlumax; 00191 } 00192 } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size()); 00193 00194 ++num_expansions; 00195 return 0; 00196 00197 } // end LuMemInit 00198 00208 template <typename Scalar, typename Index> 00209 template <typename VectorType> 00210 Index SparseLUImpl<Scalar,Index>::memXpand(VectorType& vec, Index& maxlen, Index nbElts, MemType memtype, Index& num_expansions) 00211 { 00212 Index failed_size; 00213 if (memtype == USUB) 00214 failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 1, num_expansions); 00215 else 00216 failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 0, num_expansions); 00217 00218 if (failed_size) 00219 return failed_size; 00220 00221 return 0 ; 00222 } 00223 00224 } // end namespace internal 00225 00226 } // end namespace Eigen 00227 #endif // EIGEN_SPARSELU_MEMORY