Anasazi Version of the Day
Tsqr_Random_NormalGenerator.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef __TSQR_Random_NormalGenerator_hpp
00030 #define __TSQR_Random_NormalGenerator_hpp
00031 
00032 #include <Tsqr_Lapack.hpp>
00033 #include <algorithm>
00034 #include <vector>
00035 
00038 
00039 namespace TSQR {
00040   namespace Random {
00041 
00042     template< class Ordinal, class Scalar >
00043     class NormalGenerator {
00044     private:
00045       static const int defaultBufferLength = 100;
00046 
00047     public:
00048       typedef Ordinal ordinal_type;
00049       typedef Scalar scalar_type;
00050 
00063       NormalGenerator (const std::vector<int>& iseed,
00064            const int buffer_length = defaultBufferLength) :
00065   iseed_ (4),
00066   buffer_ (buffer_length),
00067   buffer_length_ (buffer_length),
00068   cur_pos_ (0)
00069       {
00070   std::copy (iseed.begin(), iseed.end(), iseed_.begin());
00071   fill_buffer ();
00072       }
00073 
00074 
00084       NormalGenerator (const int buffer_length = defaultBufferLength) :
00085   iseed_ (4),
00086   buffer_ (buffer_length),
00087   buffer_length_ (buffer_length),
00088   cur_pos_ (0)
00089       {
00090   iseed_[0] = 0;
00091   iseed_[1] = 0;
00092   iseed_[2] = 0;
00093   iseed_[3] = 1;
00094   fill_buffer ();
00095       }
00096       
00100       Scalar operator() () { return next(); }
00101 
00104       void 
00105       getSeed (std::vector<int>& iseed) const
00106       {
00107   std::copy (iseed_.begin(), iseed_.end(), iseed.begin());
00108       }
00109 
00110     private:
00111       std::vector< int > iseed_;
00112       std::vector< Scalar > buffer_;
00113       int buffer_length_, cur_pos_;
00114 
00115       void
00116       fill_buffer () 
00117       {
00118   LAPACK< int, Scalar > lapack;
00119 
00120   // LAPACK's _LARNV routine defines this "enum" (just an
00121   // integer, because it's Fortran) that lets users choose from
00122   // one of three different pseudorandom distributions:
00123   // uniform(0,1), uniform(-1,1), and normal(0,1).
00124   enum distribution_type { 
00125     uniform_0_1 = 1, 
00126     uniform_m1_1 = 2, 
00127     normal_0_1 = 3 
00128   };
00129   lapack.LARNV (normal_0_1, &iseed_[0], buffer_length_, &buffer_[0]);
00130       }
00131 
00132       Scalar 
00133       next () 
00134       { 
00135   // Greater-than impossible, but we check for robustness' sake.
00136   if (cur_pos_ >= buffer_length_) 
00137     {
00138       fill_buffer ();
00139       cur_pos_ = 0;
00140     }
00141   return buffer_[cur_pos_++];
00142       }
00143     };
00144   } // namespace Random
00145 } // namespace TSQR
00146 
00147 
00148 #endif // __TSQR_Random_NormalGenerator_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends