SparseLU_Memory.h
Go to the documentation of this file.
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 = 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   try 
00081   {
00082     vec.resize(new_len); 
00083   }
00084   catch(std::bad_alloc& )
00085   {
00086     if ( !num_expansions )
00087     {
00088       // First time to allocate from LUMemInit()
00089       throw; // Pass the exception to LUMemInit() which has a try... catch block
00090     }
00091     if (keep_prev)
00092     {
00093       // In this case, the memory length should not not be reduced
00094       return new_len;
00095     }
00096     else 
00097     {
00098       // Reduce the size and increase again 
00099       Index tries = 0; // Number of attempts
00100       do 
00101       {
00102         alpha = (alpha + 1)/2;
00103         new_len = Index(alpha * length);
00104         try
00105         {
00106           vec.resize(new_len); 
00107         }
00108         catch(std::bad_alloc& )
00109         {
00110           tries += 1; 
00111           if ( tries > 10) return new_len; 
00112         }
00113       } while (!vec.size());
00114     }
00115   }
00116   //Copy the previous values to the newly allocated space 
00117   if (nbElts > 0)
00118     vec.segment(0, nbElts) = old_vec;   
00119    
00120   
00121   length  = new_len;
00122   if(num_expansions) ++num_expansions;
00123   return 0; 
00124 }
00125 
00138 template <typename Scalar, typename Index>
00139 Index SparseLUImpl<Scalar,Index>::memInit(Index m, Index n, Index annz, Index lwork, Index fillratio, Index panel_size,  GlobalLU_t& glu)
00140 {
00141   Index& num_expansions = glu.num_expansions; //No memory expansions so far
00142   num_expansions = 0; 
00143   glu.nzumax = glu.nzlumax = (std::max)(fillratio * annz, m*n); // estimated number of nonzeros in U 
00144   glu.nzlmax = (std::max)(Index(4), fillratio) * annz / 4; // estimated  nnz in L factor
00145 
00146   // Return the estimated size to the user if necessary
00147   Index tempSpace;
00148   tempSpace = (2*panel_size + 4 + LUNoMarker) * m * sizeof(Index) + (panel_size + 1) * m * sizeof(Scalar);
00149   if (lwork == emptyIdxLU) 
00150   {
00151     Index estimated_size;
00152     estimated_size = (5 * n + 5) * sizeof(Index)  + tempSpace
00153                     + (glu.nzlmax + glu.nzumax) * sizeof(Index) + (glu.nzlumax+glu.nzumax) *  sizeof(Scalar) + n; 
00154     return estimated_size;
00155   }
00156   
00157   // Setup the required space 
00158   
00159   // First allocate Integer pointers for L\U factors
00160   glu.xsup.resize(n+1);
00161   glu.supno.resize(n+1);
00162   glu.xlsub.resize(n+1);
00163   glu.xlusup.resize(n+1);
00164   glu.xusub.resize(n+1);
00165 
00166   // Reserve memory for L/U factors
00167   do 
00168   {
00169     try
00170     {
00171       expand<ScalarVector>(glu.lusup, glu.nzlumax, 0, 0, num_expansions); 
00172       expand<ScalarVector>(glu.ucol,glu.nzumax, 0, 0, num_expansions); 
00173       expand<IndexVector>(glu.lsub,glu.nzlmax, 0, 0, num_expansions); 
00174       expand<IndexVector>(glu.usub,glu.nzumax, 0, 1, num_expansions); 
00175     }
00176     catch(std::bad_alloc& )
00177     {
00178       //Reduce the estimated size and retry
00179       glu.nzlumax /= 2;
00180       glu.nzumax /= 2;
00181       glu.nzlmax /= 2;
00182       if (glu.nzlumax < annz ) return glu.nzlumax; 
00183     }
00184     
00185   } while (!glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size());
00186 
00187   
00188   
00189   ++num_expansions;
00190   return 0;
00191   
00192 } // end LuMemInit
00193 
00203 template <typename Scalar, typename Index>
00204 template <typename VectorType>
00205 Index SparseLUImpl<Scalar,Index>::memXpand(VectorType& vec, Index& maxlen, Index nbElts, MemType memtype, Index& num_expansions)
00206 {
00207   Index failed_size; 
00208   if (memtype == USUB)
00209      failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 1, num_expansions);
00210   else
00211     failed_size = this->expand<VectorType>(vec, maxlen, nbElts, 0, num_expansions);
00212 
00213   if (failed_size)
00214     return failed_size; 
00215   
00216   return 0 ;  
00217 }
00218 
00219 } // end namespace internal
00220 
00221 } // end namespace Eigen
00222 #endif // EIGEN_SPARSELU_MEMORY


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:39:23