Eigen  3.3.3
SparseLU_Memory.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends